Gauß-Verfahren - Matrix lösen

[Bjoerk]

Japp. Selbst mein Blödgauss aus den 80gern schafft in 1 sec eine Matrix 500 x 500. Ich denke mal TE würde sowas womöglich schon reichen. Etwas umformuliert, erste Tests, die Quelle hab ich oben angegeben.

Delphi-Quellcode:

function FloatToFracStr(const Value: Extended): string; // Quelle "delfiphan"
const
  Eps = 1E-12; // Fehlertoleranz
var
  P, LastP, Q, LastQ, TempP, TempQ, U, A, D: Extended;
  Numerator, Denominator: int64;
begin
  Numerator := 0;
  Denominator := 0;
  // Initialisierung
  A := 1;
  P := 1;
  Q := 0;
  LastP := 0;
  LastQ := 1;
  U := Value;
  // Abbruchkriterien
  while (CompareValue(U, 0) <> 0) and (CompareValue(Value + A, Value) <> 0) and (CompareValue(A, Eps) >= 0) do
  begin
    // Einen ganzzahligen Anteil abspalten
    D := Round(U);
    U := U - D;
    // Update von P und Q: Kettenbruch nachführen. Es gilt: P / Q ~= Value
    TempP := P * D + LastP;
    TempQ := Q * D + LastQ;
    LastP := P;
    LastQ := Q;
    P := TempP;
    Q := TempQ;
    // Approximationsfehler
    A := 0.25 * Abs(P / Q - Value);
    // Bruch umkehren
    if U <> 0 then U := 1 / U;
  end;
  // Vor Integerkonversion auf Bereich überprüfen
  if (P > High(int64)) or (Q > High(int64)) or (P < Low(int64)) or (P < Low(int64)) then
    raise EIntOverflow.Create('FloatToFrac: int64 overflow.');
  // Vorzeichen von Nenner zum Zähler
  if Q < 0 then
    Numerator := -Trunc(P)
  else
    Numerator := Trunc(P);
  Denominator := Abs(Trunc(Q));
  if Denominator = 1 then
    Result := IntToStr(Numerator)
  else
    Result := Format('%d/%d', [Numerator, Denominator]);
end;

[Danny92]

Zitat von jfheins:

Zum Thema: Du musst auf jeden Fall mit dem String-gerechne aufhören Vll. als erste Maßnahme mal Int64 für Zähler/Nenner hernehmen, oder wenn es WIRKLICH rational bleiben muss, dann auf Byte-Arrays für Zähle/Nenner gehen. Oder halt etwas Gehirnschmalz in eine eigene Div/Mod Funktion für Extended stecken und doch Gleitkommazahlen verwenden.

Natürlich ist es möglich, eigene Div und Modulo-Funktionen mit extended zu schreiben, wie:

Code:

function AModB(A, B: extended): integer;
begin
  result := Trunc(A-Trunc(A/B)*B);
end;

Dieser ist zwar mit bis auf extended ausdehnbare Zahlen anwendbar, aber auch nicht zu gebrauchen, da er mir bei Zahlen A und B jenseits von integer leider falsche Ergebnisse liefert, z. B. ist A=1000000000733131 MOD B=429596729935 = 328410174386, jedoch liefert mir diese Funktion in diesem Falle einen falschen Rest von result=1992659890.

Ich hätte außerdem doch von Anfang an mit extended gerechnet und mir die aufwändigen Strings gespart, wenn ich da die Ungenauigkeit nicht hätte, denn spätestens ab der 19-20ten Nachkommastelle ist bei extended Schluss.
Habe ich eine 40x40-Dreiecksmatrix, deren 40 Einträge in der Hauptdiagonalen alle 1/40 sind, ist die Determinante D=(1/40)^40 bzw. 1/40^40, und da sind wir schon weit unterhalb von 20 Nachkommastellen. Extended würde mir dann also Null liefern, was ja aber nicht richtig ist, und 1/40 ist eindeutig nicht Nichts sondern immerhin 0,025!

Ich schätze ich bin einfach auf einen ganzzahligen Datentypen angewiesen, insofern hoffe ich, dass ich mit BigInt mehr erreichen kann.
Wie gesagt sind nur Brüche exakt, Gleitkommazahlen nützen mir nichts.
Und mit extended zu rechnen, um diese dann mittels FloatToFrac (Eps = 1E-12) in Brüche umzuformen nützt mir ebenfalls nichts, denn aus falschen Dezimalzahlen entstehen falsche Brüche!
Brüche in Dezimalzahlen umwandeln ist problemlos, sofern es die Genauigkeit zulässt - umgekehrt aber nicht.
Je größer die Matrizen, umso wichtiger scheint die Rechengenauigkeit zu werden.
Relativ Schnell, (scheinbar) beliebig groß, aber vor allem eindeutig UND exakt - das ist das Problem.

Integer (2^31-1) ist zu klein - und extended ist ungenau. There is a need of an data type with a range of extended - paired with a precision of integer.

In diesem Sinne sagt Fritz zu seinen Freund Karl: "Ich kann ganz schnell multiplizieren. Nenne mir zwei Zahlen!" Da fragt Karl ihn: "Was ist 17*13?" Da antwortet Fritz: "200." Nach kurzem Überlegen meint Karl: "Das ist aber falsch!" "Na und?", meint Fritz, "aber schnell."

Und ja, Richtigkeit sollte eine höhere Priorität haben als Schnelligkeit.

[Bjoerk]

Sagt Karl, keine Ahnung was rauskommt, Ergebnis kommt übermorgen..

Ich wäre dir übrigens dankbar wenn du meine Posts nicht mit LOL kommentiert. Ein Eps von 1E-12 ist durchaus üblich und für den Rest der Welt in den allermeisten Fällen ausreichend..

[Namenloser]

Zitat von Danny92:

Habe ich eine 40x40-Dreiecksmatrix, deren 40 Einträge in der Hauptdiagonalen alle 1/40 sind, ist die Determinante D=(1/40)^40 bzw. 1/40^40, und da sind wir schon weit unterhalb von 20 Nachkommastellen. Extended würde mir dann also Null liefern, was ja aber nicht richtig ist,

Nein, würde es nicht. Die kleinste Zahl, die Extended darstellen kann ist 1.9*10^−4951 (laut Wikipedia). 1/40^40 ist gerade mal 8.2718061*10^-65 (laut Google). Da bist du also noch weit von den Grenzen entfernt. Bei dir liegt offenbar ein Missverständnis vor, wie Gleitkommazahlen funktionieren. Lies dich dazu schlau.

Wie wäre es außerdem, wenn du mal verrätst, wofür du das ganze brauchst? Als legitimen Anwendungsfall kann ich mir höchstens noch irgendwas im Bereich Kryptographie vorstellen. Aber für mich sieht es momentan eher so aus, als ob es komplett von hinten durch die Brust ins Auge ist.

[Danny92]

Zitat von Namenloser:

Zitat von Danny92:

Habe ich eine 40x40-Dreiecksmatrix, deren 40 Einträge in der Hauptdiagonalen alle 1/40 sind, ist die Determinante D=(1/40)^40 bzw. 1/40^40, und da sind wir schon weit unterhalb von 20 Nachkommastellen. Extended würde mir dann also Null liefern, was ja aber nicht richtig ist,

Nein, würde es nicht. Die kleinste Zahl, die Extended darstellen kann ist 1.9*10^−4951 (laut Wikipedia). 1/40^40 ist gerade mal 8.2718061*10^-65 (laut Google). Da bist du also noch weit von den Grenzen entfernt. Bei dir liegt offenbar ein Missverständnis vor, wie Gleitkommazahlen funktionieren. Lies dich dazu schlau.

Wie wäre es außerdem, wenn du mal verrätst, wofür du das ganze brauchst? Als legitimen Anwendungsfall kann ich mir höchstens irgendwas im Bereich Kryptographie vorstellen. Aber ich tippe eher darauf, dass es komplett von hinten durch die Brust ins Auge ist.

Danke für den Hinweis meines Missverständnisses. Aber ich habe mich wohl nur falsch ausgedrückt. Der Kern meiner Aussage sollte darin liegen, dass extended eben nur bis auf maximal 19 Nachkommastellen genau arbeitet (so viel weiß sogar ich).
Sollte also die Determinante D also einmal D=1/40^40 groß sein, wäre diese durch den Datentypen extended nur unzureichend erfasst, denn alle folgenden Nachkommastellen würden durch Null ersetzt (ja das ist so). Würde man dann mit diesen nicht exakten Werten noch weiterrechnen, würden sich diese Ungenauigkeiten (je größer die Matrix ist) noch potenzieren - Stichwort Numerische Stabilität.

Die Antwort zu deiner letzten Frage steht im Titel dieser Diskussion: Matrizen lösen. Aber man könnte auch etliche weitere Anwendungen nennen, Verschlüsselung/Primzahlen wäre nur eine von vielen.
Eine Genauigkeit von 19 Nachkommastellen mag für Mittelschulmathematik wie 1+1 reichen, aber schon dieses Beispiel zeigt schon, dass man sehr schnell an den informationstechnischen Grenzen gelangt wenn man einen großen ganzzahligen Datentypen benötigt - zumindest vorerst in Delphi. Da ich viele mathematische Problemstellungen in Programme implementiere, bin ich auf dieses Problem des zu kleinen Integers und zu ungenauen Extended schon öfter gestoßen.

[Namenloser]

Zitat von Danny92:

Sollte also die Determinante D also einmal D=1/40^40 groß sein, wäre diese durch den Datentypen extended nur unzureichend erfasst, denn alle folgenden Nachkommastellen würden durch Null ersetzt (ja das ist so).

Nein, das ist eben nicht so. Es gibt immer 19 signifikante Stellen. 1/40 wird mit der gleichen relativen Präzision repräsentiert wie 1/40^100. Es ist egal in welcher Größenordnung die Zahl ist. Wenn man 1/40^100 naiv berechnet als 1/40*1/40*..., verliert man zwar u.U. durch die Berechnung etwas Genauigkeit, aber das macht kaum etwas aus. Problematisch sind bei Gleitkommazahlen eigentlich nur Additionen und Subtraktionen.

Zitat von Danny92:

Im Noch extremeren Falle von D=(1/100)^40=1/100^40 wäre die Determinante fälschlicherweise gleich Null!

Nein, immer noch nicht. (1/100)^40 = 1*10^-80

Kannst es gerne testen:

Delphi-Quellcode:

procedure TForm1.Button1Click(Sender: TObject);
var
  c,x: Extended;
  i: integer;
begin
  c := 1/100;
  x := 1;
  for i := 1 to 40 do
    x := c*x;
  ShowMessage(FloatToStr(x));
end;

[Danny92]

Zitat von Namenloser:

Zitat von Danny92:

Sollte also die Determinante D also einmal D=1/40^40 groß sein, wäre diese durch den Datentypen extended nur unzureichend erfasst, denn alle folgenden Nachkommastellen würden durch Null ersetzt (ja das ist so).

Nein, das ist eben nicht so. Es gibt immer 19 signifikante Stellen. 1/40 wird mit der gleichen relativen Präzision repräsentiert wie 1/40^100. Es ist egal in welcher Größenordnung die Zahl ist. Wenn man 1/40^100 naiv berechnet als 1/40*1/40*..., verliert man zwar u.U. durch die Berechnung etwas Genauigkeit, aber das macht kaum etwas aus. Problematisch sind bei Gleitkommazahlen eigentlich nur Additionen und Subtraktionen.

Das habe ich auch nie bestritten. Im Gegenteil. Ich sagte ja, extended bietet 19 signifikante Stellen, aber eben nicht mehr und damit nicht alle! Und da liegt der Hund begraben. Ich möchte lediglich, dass der Bruch 1/3 eineindeutig und exakt dargestellt und vor allem mit ihm gerechnet wird, und nicht mit 1/3=0,333333333333333333300000000000000000000.....
Wenngleich Brüche stets abbrechend oder periodisch sein müssen (per Definition), wäre bei einem großen Bruch die Länge der Periode unter Umständen größer als 19, und genau dann wäre es nicht mehr exakt da extended alle weiteren Nachkommastellen durch Null ersetzt!

Code:

procedure TForm1.Button1Click(Sender: TObject);
var
  zahl: extended;
begin
  zahl:=0.3333333333333333;  //16 Nachkommastellen
  showmessage(floattostr(3*zahl)) //result=1 statt 0.9999999999999999
end;

Und diese minimalen Fehler können sich bei der weiteren Berechnung gravierend potenzieren - was ich zu vermeiden versuchen möchte, deshalb die Sache mit den Strings.
Man bräuchte also einen Gleitkommadatentypen der unendlich viele Nachkommastellen erfasste (Beispiel 1/3), da das aber natürlich nicht geht, kann ich mit Gleitkommazahlen nichts anfangen, also benötige ich einen Ganzzahligen in dieser Größenordnung.

Ich wünschte lediglich, es gäbe einen ganzzahligen Datentypen wie Integer, der aber bisschen weiter hinaus reicht, denn 2^31-1 ist nun wirklich nichts! (Im mathematischen Sinne ist jede Zahl klein, aber z. B. 2^500 wäre ja wirklich schon eher brauchbar als nur 2,1 Milliarden, sogar Banker können sich größere Zahlen vorstellen)

[Namenloser]

Du wirst dich von der Idee der Exaktheit sowieso früher oder später verabschieden müssen. Rationale Zahlen kannst du zwar noch exakt darstellen, aber spätestens bei den reellen Zahlen hört es eh auf. Was machst du, wenn in deinem Gleichungsssystem z.B. die Zahl Pi oder √2 vorkommt? Auch musst du bedenken, dass oft die Eingangsdaten schon nicht exakt sind, weil es sich z.B. um Messwerte handelt.

Zitat:

Und diese minimalen Fehler können sich bei der weiteren Berechnung gravierend potenzieren

Wenn man es ungeschickt anstellt, ja. Wenn man es geschickt anstellt, lassen sich Ungenauigkeiten minimieren. QR-Zerlegung und LUP-Zerlegung (letztere kenne ich nicht) wurden bereits genannt, aber du bist bisher nicht darauf eingegangen.

Das wäre ein besserer Ansatz als Unmengen von Performance zu verschwenden, um sich in scheinbarer Genauigkeit zu wiegen, die einem am Ende sowieso nichts nützt.

Zitat:

Ich wünschte lediglich, es gäbe einen ganzzahligen Datentypen wie Integer, der aber bisschen weiter hinaus reicht, denn 2^31-1 ist nun wirklich nichts

Dann nimm Int64.

[Dejan Vu]

Diese Diskussion à la 'Extended ist genau genug' ist Quark. Entschuldigung. Beim Rechnen mit sehr kleinen UND sehr großen Zahlen kommen falsche Ergebnisse heraus, außer! man rechnet genau. Und das kann man nun einmal nur mit exakten Brüchen machen. Ergo benötigt man als 'Zahl' Datentyp z.B. einen Record (oder ne Klasse) mit Zähler und Nenner. Dann definiert man noch die Operationen auf diesem Zahl-Datentyp und -wupps- werden die Ergebnisse genau.

Ich habe -wie erwähnt- mit Gauß und Extened bei Berechnungen zur Kälteleistung von Kühlkompressoren keine guten Erfahrungen gemacht. Es ist blöd, wenn die Kennlinien statt im Bereich von 600W dann bei -170 liegen, weil klitzekleine Ungenauigkeiten beim Rechnen passiert sind. Gut, ich hatte Glück, das LUP Dekomposition hier geholfen hat. Aber grundsätzlich sollte man das schon anders lösen.

Ich unterstelle dem TE im übrigen, sich in der Zahlentheorie und mit Genauigkeitsrechnungen auszukennen. Ich denke, da muss man nicht mehr den Besserwisser raushängen lassen. Natürlich ist niemand gemeint, ich habe das nur präventiv gesagt

[Danny92]

Zitat von Namenloser:

Du wirst dich von der Idee der Exaktheit sowieso früher oder später verabschieden müssen. Rationale Zahlen kannst du zwar noch exakt darstellen, aber spätestens bei den reellen Zahlen hört es eh auf.

Ich lasse nur rationale Zahlen als Eingabe zu. Sprich ganze Zahl / ganze Zahl. Kein pi oder E oder sonstiges. Und egal welche Rechenoperation (+,-,*,/), das Ergebnis ist immer auch eine rationale Zahl.
Und Int64 kann ich leider nicht nehmen, da dies kein Ordinaldatentyp im Gegensatz zu Integer ist.
Mit dem Methoden QR und LUP bin ich bisher noch nicht vertraut, deswegen hatte ich mir die Berechnung erst einmal ohne weiteren Aufwand erhofft.
Ich werde das mit dem BigInt mal probieren.