Ceil Floor Log10 - Rundungsfehler, Überlauf o.ä. (?)

[hansmaad]

Hallo zusammen,
ich steh gerade total auf dem Schlauch und kann mir den folgenden Fehler nicht erklären:

Delphi-Quellcode:

procedure Test(min, max : Double);
var
    x : Integer;
    a, b : Double;
begin
    x := Ceil(Log10(max)) - Floor(Log10(min)) + 1;
    WriteLn(x);

    a := Log10(max);
    b := Log10(min);
    x := Ceil(a) - Floor(b) + 1;
    WriteLn(x);
end;

begin
    Test(1e-2, 1e8);
    Test(1e-2, 1e7);
    Test(1e-2, 1e6);
end.

Code:

11
11
11
10
9
9

Die Ausgabe sollte natürlich sein:

Code:

11
11
10
10
9
9

[shmia]

Was passiert, wenn du den Variablen a und b den Datentyp Extended gibst?

[hansmaad]

Code:

11
11
11
11
9
9

Jetzt bin ich noch verwirrter.

[Jens01]

Ich weiß nicht, ob es hilft, aber guck auch mal hier:
http://home.netsurf.de/wolfgang.ehrh...tml#amath_unit

[hansmaad]

Zitat:

Zitat von Jens01 (Beitrag 1101528)

Ich weiß nicht, ob es hilft, aber guck auch mal hier:
http://home.netsurf.de/wolfgang.ehrh...tml#amath_unit

Ohne nachvollzogen zu haben, ob das jetzt mit dem Fehler hier zu tun hat...:shock::shock::shock:

[Jens01]

Naja, ich hab mir Dein Problem nur kurz angeguckt -ist ja nicht mein Problem:wink:.
Aber ich weiss, dass in amath einige Fehler von math ausgeräumt wurden, auch Rundungsprobleme. Vielleicht fällt das bei Dir ja auch darunter.
Ausprobieren bzw gucken muß Du aber selbst, oder es lassen.

[hansmaad]

Ich glaub du hast meine Antwort falsch verstanden. Ich meinte ich habe zwar nicht nachvollzogen, ob das hier mit dem Thread zusammenhängt, aber ich bin absolut geschockt. Hab mich zwar inzwischen daran gewöhnt, ab und zu mal einen halben Tag wegen irgendwelchen blöden Delphi Bugs zu verschwenden, aber dass solche elementaren Dinge seit 15 Jahren so falsch sind, ist ja echt krass.

[rollstuhlfahrer]

Zitat:

Zitat von hansmaad (Beitrag 1101539)

Hab mich zwar inzwischen daran gewöhnt, ab und zu mal einen halben Tag wegen irgendwelchen blöden Delphi Bugs zu verschwenden, aber dass solche elementaren Dinge seit 15 Jahren so falsch sind, ist ja echt krass.

Eigentlich wollte ich dich gerade bremsen, nicht über die Bugs bei Fließkommazahlen herzuziehen und habe mal einen Test mit C# gemacht. Siehe da: Dort kommt das richtige Ergebnis raus. Da intern eigentlich CPU-Zahlen (oder FPU) verwendet werden sollten, hätte das doch eigentlich stimmen müssen, was das Delphi-Programm ausgibt, aber scheinbar tut es das nicht. Komisch.

Code:

namespace Test_Rundung
{
   class Program
   {
      
      
      static void Test(double min, double max) {
         int x;
         double a;
         double b;
         
         x = Convert.ToInt32(Math.Ceiling(Math.Log10(max))) - Convert.ToInt32(Math.Floor(Math.Log10(min))) + 1;
         Console.WriteLine(x);
         
         a = Math.Log10(max);
          b = Math.Log10(min);
          
          x = Convert.ToInt32(Math.Ceiling(a)) - Convert.ToInt32(Math.Floor(b)) + 1;
         Console.WriteLine(x);
      }

   
   
      public static void Main(string[] args)
      {
         Console.WriteLine("Hello World!");
         
          Test(1e-2, 1e8);
          Test(1e-2, 1e7);
          Test(1e-2, 1e6);

         
         Console.Write("Press any key to continue . . . ");
         Console.ReadKey(true);
      }
   }
}

Versuche mal folgendes: Delphi neu starten, Programm neu kompilieren und wenn das nicht reicht: Rechner neu starten.

Bernhard

PS: Die Einrückung der Forensoftware ist wieder brutal toll.

[hansmaad]

Ich hatte den gleichen Test auch vorher mit vc und gcc gemacht. Aber da gibt ceil und floor ja (wie c#) eine Fließkommazahl zurück, daher dachte ich, ich würde bei Delphi irgendwas übersehen.

[gammatester]

Deine Idee ist zwar ein Disaster-Konzept, weil Differenzen von zwei Sprungfunktionen geradezu danach schreien, solche Effekte hervorzurufen; in Deinem Fall für Werte, die keine Zehnerpotenzen sind.

Interessant ist jedoch an Deinem Problem folgendes:

Es sieht so aus als wenn Du auf eine Inkonsistenz der FPU gestoßen bist. log10(x) wird im Prinzip als log10(2)*log2(x) berechnet. Wenn man das so programmiert, verschwindet auch Dein Problem. Etwas mehr Hintergrund-Info nach einiger Knobelei:

Delphi-Quellcode:

{---------------------------------------------------------------------------}
function log10(x: extended): extended; assembler;
  {-Return base 10 logarithm of x}
asm
  fldlg2
  fld     [x]
  fyl2x
  fwait
end;

{---------------------------------------------------------------------------}
function log10a(x: extended): extended; assembler;
  {-Return base 10 logarithm of x}
asm
  fld1
  fld     [x]
  fyl2x
  fldlg2
  fmul
end;

Die beiden Funktionen sollten eigentlich dieselben Ergebisse bringen, log10(x) rechnet log10(2) * log2(x) via FPU-Befehl fyl2x und log10a(x) rechnet (1*log2(x)) * log10(2). Für x=1e7 rechnet log10a richtig und log10 (die Delphi-Implementation) um 1 ulp zu hoch (die $-Zahlen sind die internen extended Darstellungen via AMath und MPArith)

Code:

log10a(1e7) = $4001E000000000000000 = 7.0
log10(1e7)  = $4001E000000000000001 = 7.00000000000000000043368086899420177360298112034797668457031

Gruß Gammatester

[idefix2]

Zitat:

Deine Idee ist zwar ein Disaster-Konzept, weil Differenzen von zwei Sprungfunktionen geradezu danach schreien, solche Effekte hervorzurufen

Das ist der springende Punkt: Gleitkommaarithmetik kann von der Systematik her keine ganz exakten Werte liefern, und schon gar nicht, wenn wie hier eine transzendentale Funktion wie der Dezimallogarithmus aufgerufen wird. MINIMALE Abweichungen vom mathematisch exakten Ergebnis sind da nicht überraschend, sonden völlig normal und zu erwarten. Diese Abweichungen sind in einer Grössenordnung, die normalerweise keine Rolle spielt.

Wenn Du jetzt aber als Ergebnis einer solchen Funktion eine ganze Zahl erwartest, und es kommt etwas heraus, das im Rahmen der erlaubten Abweichung kleiner ist als die erwartete Zahl, dann wirkt sich die danach angewendete Funktion Floor natürlich katastrophal aus, weil der Dezimalteil, in dem Fall vielleicht 0.99999999999999999999999999, abgeschnitten wird. Daß das aber kein Fehler von Delphi, sondern ein Fehler im Programmierkonzept ist, muss jedem klar sein, der sich auch nur im Entfertnesten mit numerischer Mathematik und der Rundungsproblematik beschäftigt hat.

Zitat:

Siehe da: Dort kommt das richtige Ergebnis raus.

Ob das Ergebnis der log10 Funktion 7.000000000000000001 oder 7.000000000000000000 oder 6.999999999999999999999 liefert, ist ein reines Lotteriespiel - in den beiden ersten Fällen liefert floor(log10(x)) dann die erwartete 7, im anderen Fall eben 6. Das Ergebnis ist deshalb nicht in einem Fall richtig und im anderen Fall falsch (in Hinblick auf die Umsetzung durch den Compiler), sondern in allen Fällen im Rahmen der von einer Gleitkommaarithmetik erwartbaren Genauigkeit.

Das Problem wirst Du auch nicht los, indem Du die Rechengenauigkeit erhöhst, weil, egal auf wieviele Stellen Du rechnest, in der allerletzten Binärstelle eine Abweichung sein wird, und ganz egal wie klein die Abweichung ist, wenn sie zufällig nach unten geht, schneidet floor den Dezimalteil weg.

[hansmaad]

Also erstmal Danke für die hilfreichen Hinweise.
Vielleicht können wir dann ja noch klären, was die richtige Implementierung für das Problem ist. Gesucht ist eine Funktion, die die nächst niedrige und eine Funktion, die die nächst höhere Zehnerpotenz zu einer Gleitkommazahl zurückgibt. Und zwar als Ganzzahl, die den Exponenten angibt.

[Jens01]

Delphi-Quellcode:

Floor(Round(Log10(min)))

[idefix2]

Es gibt zu Deiner Frage zwei Probleme, die man berücksichtigen muss:
Das erste hast Du im Ursprungsposting angesprochen, nämlich die Genauigkeit der Funktion selbst. Die Funktion ist nichtlinear, sie hat bei jeder Zehnerpotenz eine Sprungstelle. Nachdem auch Gleitkommaarithmetik exakt rechnet, solange nur ganze Zahlen und die Grundrechenarten involviert sind, kannst Du in einer Schleife abfragen, ob Dein x grösser als 1,10,100 etc. ist, eine einfache Multiplikation einer ganzen Realzahl mit 10 wird sicher keine Rundungsfehler produzieren. Bei den Dezimalzahlen 0.1, 0.01 ist die Angabe des Sprungpunkts aber nicht exakt möglich, weil diese Zahlen in Binärdarstellung eine periodische Mantisse haben.

Damit hast Du aber das zweite Problem noch nicht gelöst, nämlich die Genauigkeit des Funktionsparameters x selbst: Wo kommt der her? Wenn der das Ergebnis einer Geitkommaberechnung ist, dann kann sich auch dort der kleinste Rundungsfeher dahingehend auswirken, dass bei einer Zahl, die z.B. ganz knapp kleiner als 100 sein sollte, etwas herauskommt, was ganz knapp grösser als 100 ist und deshalb das Ergebnis der Gesamtberechnung einen unerwarteten Wert liefert. Oder anders herum, der Parameter, den Du übergibst, hätte genau 10 sein sollen, ist aber als Ergebnis einer Gleitkommarechnung 9.999999999999999999999999999999994, und das Ergebnis Deiner Funktion stimmt deshalb nicht mit dem, was Du erwartest, überein, obwohl die Funktion selbst ganz genau und ohne Rundungsfehler rechnet.

Zitat:

Floor(Round(Log10(min)))

Das wird ganz falsche Ergebnisse bringen, floor(round(0.7))=1 statt der erwarteten 0. Floor auf einen gerundeten Wert anzuwenden ist übrigens sinnlos, weils der gerundete Wert ja schon eine ganze Zahl ist und der nachfolgende floor-Aufruf nichts mehr verändert.

[gammatester]

Hier ein Vorschlag, der mit einer Toleranz für den log10 arbeitet. Die Funktionen schauen nach, ob x nahe bei einer Zehnerpotenz liegt. Wenn nein, werden die nächst niedrigeren bzw. höheren Potenzen genommen. Die Nähe-Toleranz tol sollte sinnvollerweise nicht kleiner als ca 1e-15 gewählt werden.

Delphi-Quellcode:

{---------------------------------------------------------------------------}
function minpower10(x: extended): integer;
  {-Berechnet für x>0 den maximalen Integerwert mit 10^result <= x}
const
  tol = 1e-9;
var
  lx: extended;
  li: integer;
begin
  lx := log10(x);
  li := round(lx);
  if abs(li - lx) > tol then li := floor(lx);
  minpower10 := li;
end;


{---------------------------------------------------------------------------}
function maxpower10(x: extended): integer;
  {-Berechnet für x>0 den minimalen Integerwert mit 10^result >= x}
const
  tol = 1e-9;
var
  lx: extended;
  li: integer;
begin
  lx := log10(x);
  li := round(lx);
  if abs(li - lx) > tol then li := ceil(lx);
  maxpower10 := li;
end;

[idefix2]

Das geht einfacher, indem Du für die Funktion floor einfach 1E-9 zum Parameter addierest und für die Funktion Ceil einfach 1E-9 vom Parameter abziehst.