Double Ungenauigkeit

High all,

Ich habe folgendes Problem: Nachdem ich eine Klasse zum Halten und Berechnen von komplexen Zahlen erstellt habe, gibt es (zumindest) eine merkwürdige Ausgabe. exp(new Complex(0,Math.PI)) sollte eigentlich (Double) -1 bzw. (Complex) (-1,0) ergeben, und das 'E-16' verrät ja auch die Affinität des Imaginärteils zu Null, aber zufrieden bin ich damit nicht.

Relevant dafür sind folgende Code-Abschnitte:

-------------- the testing code --------------

Complex z1 = new Complex();
Console.WriteLine(z1.ToString());
z1.Imag = Math.PI;
Console.WriteLine("{0}\t{1}",z1.ToString(), z1.Exp().ToString());
z1.DoExp();
Console.WriteLine(z1.ToString());
Console.WriteLine("Cos(PI)\t\t{0}", Math.Cos(Math.PI));
Console.WriteLine("Sin(0)\t\t{0}", Math.Sin(0));
Console.WriteLine("Exp(0)\t\t{0}", Math.Exp(0));
Console.WriteLine("Sin(0)*Exp(0)\t{0}", Math.Sin(0)*Math.Exp(0));

----------------- the output -----------------

0
3,14159265358979i       -1+1,22460635382238E-16i
-1+1,22460635382238E-16i
Cos(PI)         -1
Sin(0)          0
Exp(0)          1
Sin(0)*Exp(0)   0

Die Ausgabe '0' statt '0+0i' ist so gewollt in der ToString-Methode implementiert.

---------------- constructors ----------------

public Complex() : this(0.0d, 0.0d) { }
public Complex(Double aReal, Double aImag) { fReal = aReal; fImag = aImag; }

------------------ property ------------------

public double Imag
{
  get { return fImag; }
  set { fImag = value; }
}

---------- Exp() and DoExp() methods ---------

public Complex DoExp()
{
  Complex tmpZ = new Complex(exp(this));
  fReal = tmpZ.Real;
  fImag = tmpZ.Imag;
  return this;
}
public Complex Exp()
{
  return exp(this);
}

------------ exp() static function -----------

public static Complex exp(Complex aZ1)
{
  return new Complex(Math.Cos(aZ1.Imag),Math.Sin(aZ1.Imag))
               * Math.Exp(aZ1.Real);
}

---------- Complex * Double operator ---------

public static Complex operator *(Complex aZ1, Double aR1)
{
  return new Complex(aR1 * aZ1.Real, aR1 * aZ1.Imag);
}

Soweit ich weiß, ist die exp() Funktion korrekt implementiert...

Auch die Ergebnisse von Math.Sin(), Math.Cos(), Math.Exp() sollten korrekt sein, also wieso in Teufels Namen gibt mein Testcode nicht einfach '-1' aus?

Oh, nebenbei: ich benutze MS VC# 2k5 EE. MSVS 2k5 reports:

------------- version information ------------
Microsoft Visual Studio 2005
Version 8.0.50727.42 (RTM.050727-4200)
Microsoft .NET Framework
Version 2.0.50727

Installierte Edition: C# Express

Microsoft Visual C# 2005 76544-000-0000011-00126 Microsoft Visual C# 2005

Der Thread trägt genau den richtigen Titel: Double Ungenauigkeit.
Du hast im Double-Format nun mal nur 53Bit für die Mantisse zur Verfügung. Und eine Ungenauigkeit im 53. Bit, wenn das erste den Wert 1 hat (wegen der -1), ergibt nun mal 2^(-52) = 2,22e-16. Dies ist die gleiche Größenordnung, wie deine Abweichung und zeigt uns damit, dass es sich dabei einfach um die zwangsläufig auftretenden Rundungsfehler handelt.
Das ist aber ein prinzipielles Problem von Fließkommazahlen. Die einzige Möglichkeit, so etwas bei der Ausgabe zu unterbinden, ist das Wegrunden der Rundungsfehler. So könnte man davon ausgehen, das nur 10 Dezimalstellen genau wären. Mit einer Stelle vor dem Komma muss man also auf die neunte Nachkommastelle runden.
Danach hätte man genau -1.

High nochmal,

sorry, wenn ich immer etwas später reagiere, aber ich hab nicht immer Gelegenheit, das Internet zu nutzen.

Der Thread trägt genau den richtigen Titel: Double Ungenauigkeit.

Ja, du hast ja recht... ich versuche stets, mich möglichst präzise auszudrücken 🙂 und danke, dass du so ausführlich geantwortet hast, aber das war eigentlich unnötig, da ich die Interna von IEEE 754 so im Groben schon kenne.

Punkt is: ich drücke mich doch nicht stets präzise aus 😕 Die Kernfrage war eher:

Wieso erhalte ich für die einzelnen Komponenten (sin, cos, exp von double) exakte Werte 1 und 0, nicht aber im Kompositum (exp von complex), das sich ja aus exakten Werten 0 und 1 zusammensetzt?

Eine Frage, die ich mir schon gestellt habe, ist, ob double.ToString() die Werte nicht "schönrechnet", sprich: rundet. Da allerdings complex.ToString() für Real- und Imaginärteil auf double.ToString() zurückgreift, dürfte sich dieses "Phänomen" ebenso nicht zeigen.

Das Problem ist ein mathematisches. Du berechnest Werte der Sinus- und Kosinusfunktion nahe der Stelle Pi.
Jetzt ist die Steigung von cos an dieser Stelle null, so dass geringe Abweichungen der Stelle nur geringe Abweichungen des Funktionswertes verursachen.
Im Gegensatz zum Sinus. Dieser hat an der Stelle Pi seine größte Steigung, so dass leicht abweichende Stellen eine große Abweichung des Funktionswertes verursachen.

Da du ja nicht genau mit Pi rechnest, sondern mit dem Näherungswert 3,14159265358979 rechnest, erhältst du (FPU-intern als Extended, 64 Bit Mantisse):

cos 3,14159265358979 = -0,99999999999999999999999999999476
binär: Sign: -1, Exponent: -1, Mantisse: 1x Hidden Bit + 64 gesetzte Bit
beim Berechnen in der FPU fallen zusätzliche Rundungsbits an, diese geben in diesem Fall ein Aufrunden an, daher wird auf [Sign: -1, Exponent: 0, Mantisse: 1x Hidden Bit + 64 gelöschte Bits] gerundet. Dies ergibt genau -1.

sin 3,14159265358979 = 0,0000000000000032384626433832795028841970499567
binär: Sign: +1, Exponent: -49, Mantisse: 1x Hidden Bit + [1101001010110110001...]
Hier bringen aber die internen Rundungsbits keine Veränderung in den MSB der Mantisse. Die Berechnung der ersten Bits der Mantisse von 12^(-49) + 12^(-50) + 1*2^(-51) liefert rund 3,1e-15, was ja in der Größenordnung deines Fehlers liegt.

Solcherart Probleme kann man nur dadurch beheben, dass man von einem berechnenden System zu einem auswertenden System (also eine Algebrasystem) übergeht, welches intern nicht berechnet, sondern als algrbraischen Ausdruck zusammenfasst.
Ansonsten kann man die Probleme nur umgehen, in dem man die double-Mantisse auf beispielsweise 50 Bit rundet. Leider kann man nicht direkt angeben, auf wieviel Bits gerundet werden soll, daher muss man zuerst die Dezimalstelle des 50. Bits ausrechnen:

// Dezimalstelle ausrechnen
int ValidDigits = (int)Math.Round(Math.Log10(unrounded) - 50 * Math.Log10(2) );
// leider wirft Math.Round eine Exception, wenn Stelle nicht im Intervall [0;15] => anpassen
ValidDigits = (ValidDigits < 0) ? 0 : ((ValidDigits > 15) ? 15 : ValidDigits);
double rounded = Math.Round( unrounded, ValidDigits );