Toleranzen von Double und Float

[SMaag]

@STARGATE
Das ist worauf ich hinaus wollte! Das Problem scheint aber nicht all zu bekannt zu sein. Mir ist das erst vor kurzem untergekommen.
Hab mittlerweile auch einen Artikel dazu gefunden.
https://www.elektronikpraxis.de/verglei ... -a-657052/

Jetzt kommt aber nochmals die Frage auf, warum man davon üblicher weise nichts weis und warum das dann nicht wenigstens in Programmiersprachen standardmäßig integriert ist! Noch besser wäre eigentlich eine Rundung auf Genauigkeit gleich im Prozessor.

ich hab hier jetzt noch mal eine Code gebastelt, der das Problem mit dem Floatvergleich mit Epsilopn 2e-7 zeigt. Das geht nur für Werte < 4.0

Code: Alles auswählen

; EPSILON und TINY aus Fortran Codes für diverse numerische Methoden
;   EPSLON = 0.0000001
;   Eps = 0.0000001
;   TINY=1.0e-20
;   TINY=1.e-25
;   EPS=1.e-6
;   EPS=3.e-7 MIN=1.e-30
;   EPS=6.e-8
;   EPS=1.e-10
;   TINY=1.5e-38
;   TINY=1.69e-38, SQRTNY=1.3e-19, BIG=3.E37
;   EPS=1.2e-7
;   EPS=3.e-8
;   EPS=1.e-7, TOLMIN=1.e-6
;   EPS=1.e-14
;   EPS=3.e-14
;   EPS=3.e-8

; und in einem Programm war das automatisch geregelt mit MachineFloatPrecsion

EnableExplicit

Procedure.d FloatPrecision()
  Static.f FPr
  
  ; IT CARRIES OUT AN APPROXIMATION OF MACHINE PRECISION
  FPr = 1.0  
  While (1.0 + FPR) > 1.0
    FPr = 0.5 * FPr    
  Wend
  ProcedureReturn FPr  * 2.0
EndProcedure

Define eps.f  ;

eps = FloatPrecision()

eps = 2e-7		; Test mit festem Epsilon: für automatisch berechnet, diese Zeile auskommentieren!!!!

Debug "eps = " + StrF(eps,20)
Define.i I
Define.f n, x 

n = 1
For I = 0 To 16
  
  x = n + eps
  
  If x = n
    Debug Str(I) + " : für n = " + StrF(n,1) + " ist (n + eps) = n"   
  Else
    Debug Str(I) + " : für n = " + StrF(n,1) + " ist (n + eps) <> n"
  EndIf
  n * 2
Next

und es taucht noch ein gravierendes Problem auf. Die automatische Berechung der Floatprecision funktioniert so überhaupt nicht.
Was da rauskommt scheint die Precision von doubles zu sein.
Wenn man Epsilon automatisch berechnen lässt und das dann auf eine Float Zahl addiert ändert sich diese nicht.
Das ermittelete Epsilon ist also viel zu klein. Wie kann das sein?

[SMaag]

Ok! Das erste Problem hab ich gelöst!
Das korrekte Epsilon für 32 Bit Floats. Man muss eine explizite Addition über Variablen einführen, so dass das Ergebnis
im Speicher abgelegt wird. Dann geht's!

Code: Alles auswählen

Procedure.d FloatPrecision()
  Protected.f FPr, s
  
  ; IT CARRIES OUT AN APPROXIMATION OF MACHINE PRECISION
  FPr = 1.0  
  s = FPr + 1.0
  While s > 1.0
    FPr = 0.5 * FPr
    s = FPr + 1.0
  Wend
  ProcedureReturn FPr  * 2.0
EndProcedure

Define eps.f

eps = FloatPrecision()

Debug "eps = " + eps

[STARGÅTE]

@SMaag:

Warum rechnest du das Epsilon aus? Wie ich weiter oben geschrieben hatte sind die Epsilons (im Sinne der Mantisse) konstante Werte:

Code: Alles auswählen

#PB_FloatEpsilon  = 1.1920928955078125e-7
#PB_DoubleEpsilon = 2.2204460492503130808e-16

Jetzt kommt aber nochmals die Frage auf, warum man davon üblicher weise nichts weis und warum das dann nicht wenigstens in Programmiersprachen standardmäßig integriert ist! Noch besser wäre eigentlich eine Rundung auf Genauigkeit gleich im Prozessor.

Weil man bei Fließkommazahlen üblicher Weise keine "ist gleich"-Vergleiche durchführt. Dieses Runden auf eine bestimmte Genauigkeit ist ja auch wieder nur in der 2er-Basis sinnvoll. In der für uns üblichen 10er-Basis kommt da weiterhin "Murks" raus.

Wenn man Epsilon automatisch berechnen lässt und das dann auf eine Float Zahl addiert ändert sich diese nicht.
Das ermittelete Epsilon ist also viel zu klein. Wie kann das sein?

An der Stelle ist mir immer noch nicht klar, was dein Ziel ist?
Du möchtest ein Epsilon was immer passend für die aktuelle Zahl ist oder was?
Bei 1.0 soll Epsilon 1.19e-7 sein und bei 8.0 dann 9.54e-7?
Aber selbst dafür brauchst du keine Funktion. Das passende Epsilon für eine Zahl die nicht 1 ist, ist einfach die Multiplikation dieser Zahl mit Epsilon, also 8.0 + 8.0*#PB_FloatEpsilon wäre von 8.0 verschieden, usw.

[SMaag]

Aber selbst dafür brauchst du keine Funktion. Das passende Epsilon für eine Zahl die nicht 1 ist, ist einfach die Multiplikation dieser Zahl mit Epsilon, also 8.0 + 8.0*#PB_FloatEpsilon wäre von 8.0 verschieden, usw.

ja genau! das hatte ich nicht gecheckt.


Jetzt muss ich nochmal zurück, zu dem was ich eigentlich wollte. Das ganze kam aus der Konvertierung von FOTRAN Code nach PureBasic und da ist mir 1 mal statt einem fixen Epsion, die automatische Berechnung untergekommen. Und ich wollte wissen ob man das so machen kann und ob die angegebenen fixen Epsilon in anderen Codes so überhaupt funktionieren, wenn man nicht genau weis auf welche Werte man vergleicht.

Aber ich glaub ich hab es jetzt kapiert! Ich müsste im Fortran Code nachsehen, ob eine Multiplikation von Epsilon durchgeführt wird um das korrekte Epsilon für die verglichenen Werte zu erhalten.

[SMaag]

hier nun mal meine Interpretation von einem IsEqual für Float und Double, wenn man das Espsilon automatisch bestimmt, so
dass man eine minimale Toleranz nach oben und unten hat, was Rundungsfehler ausgleicht.
#Eps1F und #Eps1D sind jeweils die minimalen Werte, die man zu einer 1 addieren muss, damit sie nicht mehr als 1 erkannt wird.

Zur Diskussion!

Code: Alles auswählen

  Procedure.i IsEqualF(A.f, B.f)
    Protected.f eps
    #Eps1F = 1.192e-7 ; 0.0000001192
    
    If A > B
      eps = A * #Eps1F  
    Else
      eps = B * #Eps1F
    EndIf
    
    If Abs(A-B) < eps
      ProcedureReturn #True
    EndIf
    
    ProcedureReturn #False
   
  EndProcedure
  
  Procedure.i IsEqualD(A.d, B.d)
    Protected.d eps
    #Eps1D =2.220446e-16  ; 0.0000000000000002220446
    
    If A > B
      eps = A * #Eps1D  
    Else
      eps = B * #Eps1D
    EndIf
    
    If Abs(A-B) < eps
      ProcedureReturn #True
    EndIf
    
    ProcedureReturn #False
   
  EndProcedure

[DarkDragon]

Jetzt muss ich nochmal zurück, zu dem was ich eigentlich wollte. Das ganze kam aus der Konvertierung von FOTRAN Code nach PureBasic und da ist mir 1 mal statt einem fixen Epsion, die automatische Berechnung untergekommen. Und ich wollte wissen ob man das so machen kann und ob die angegebenen fixen Epsilon in anderen Codes so überhaupt funktionieren, wenn man nicht genau weis auf welche Werte man vergleicht.

Du hattest irgendwo was von einem anderen Prozessor geschrieben. Das blöde mit Floats ist, sie haben einen Standard, aber es gibt durchaus Unterschiede zwischen den Prozessoren. Mal gibt es einen Bug, mal ist der Rundungsmodus anders. Ich hatte mal etwas für x86 und armv6 geschrieben und beide drifteten voneinander weg in den Berechnungen.

Hier ein bekanntes Problem, wo genau das auftritt

Das wäre vielleicht ein valider Grund es dynamisch zu berechnen.

[SMaag]

Du hattest irgendwo was von einem anderen Prozessor geschrieben.

Ja, das mit den anderen Prozessor war aus dem Fortran Code. Der war wahrscheinlich für PC und Großrechner.

Für PureBasic kommen aber nur Prozessoren in Frage, die von PB auch unterstützt werden! Da wäre es noch interessant, wie
sich das mit den Flaots auf ARM verhält.
Ist dort das Epsilon identisch und werden bei gleichen Berechnungen Bit identische Werte erzeugt?
Folgendes hab ich in einer ARM Dokumentation gefunden

Note that ARM64 normally treats C long doubles as 128-bit, but Apple doesn’t, and retains them as 64-bit. For the sake of simplicity here, wherever appropriate I’ll use D versions to work in 64-bit Doubles.

Ich gehe mal davon aus, dass das nicht ganz kompatibel mit den x87 FPUs ist!

Es kommt noch ein weiteres Problem hinzu: Abs()
"Gibt den absoluten Wert (ohne Vorzeichen) der angegebenen Fließkomma-Zahl zurück."
Das hab ich bisher ohne viel nachdenken immer auch für Integer verwendet.
Auf x87 FPUs funktioniert das ohne Probleme auch für sehr große Zahlen, da die FPU Register intern eine 64Bit Mantisse haben und
64Bit Integer somit bis ins letzte Bit darstellen können.
Wie das auf ARM aussieht weis ich leider nicht!
Meist ist das aber eh kein Problem, da man so große Zahlen nur seltenst hat! Nur wissen sollte man es!

[STARGÅTE]

Es kommt noch ein weiteres Problem hinzu: Abs()
"Gibt den absoluten Wert (ohne Vorzeichen) der angegebenen Fließkomma-Zahl zurück."
Das hab ich bisher ohne viel nachdenken immer auch für Integer verwendet.
Auf x87 FPUs funktioniert das ohne Probleme auch für sehr große Zahlen, da die FPU Register intern eine 64Bit Mantisse haben und
64Bit Integer somit bis ins letzte Bit darstellen können.

Das mit der internen 64 Bit Mantisse stimmt schon, aber im Moment der Übergabe eben nicht.
Somit kann Abs() nicht für den gesamten Wertebereich von Quads genutzt werden! Genauso wie Mod() oder Pow(). Bei all diesen Fließkomma-Funktionen sollte man immer auf ein Quad-Pendant zurückgreifen.

[SMaag]

Bei x64 funktioniert ABS ausser bei der größten negativen Zahl -9223372036854775808

Code: Alles auswählen

Define x.q = -9223372036854775808

Define.q I,K

For I = x To x + 100
  K = Abs(I)
  Debug K
Next

[STARGÅTE]

Bei x64 funktioniert ABS ausser bei der größten negativen Zahl -9223372036854775808

Code: Alles auswählen

Define x.q = -9223372036854775808

Define.q I,K

For I = x To x + 100
  K = Abs(I)
  Debug K
Next

Da würde ich mich aber nicht drauf verlassen. Die Hilfe schreibt eindeutig: "Diese Funktion verarbeitet und gibt Float oder Double Werte zurück.".
Hier zwei Gegenbeispiele:

Code: Alles auswählen

Define quad.q = Abs(-100000000000000005)
Debug quad

Define quad.q = -100000000000000005
Debug Abs(quad)

Scheinbar wurde für Quad = Abs(Quad) eine alternative ABS Funktion überladen, das ist aber nicht dokumentiert.