Toleranzen von Double und Float

[SMaag]

Ein Auszug aus der PB Hilfe
Der genaue Wertebereich, in dessen Rahmen beim Rechnen mit Floats und Doubles korrekte Ergebnisse erzielt werden, sieht wie folgt aus:
Float: +- 1.175494e-38 bis +- 3.402823e+38
Double: +- 2.2250738585072013e-308 bis +- 1.7976931348623157e+308

D.h. man müsste bei Double auf Toleranzen von < 3e-308 kommen


Versuche ich, dass direkt berechnen zu lassen, d.h. ab wann eine 1 noch als 1 erkannt wird, kommt man nur auf etwa 1e-16.
genau 0.000 000 000 000 000 111 022 3025 ; ca. 1.11e-16

Wo liegt da das Problem? Übrigens die Methode, die Genauigkeit zu berechnen stammt als alten Fortran Codes als die
x87 Coprozessoren aufkamen.

Code: Alles auswählen

Procedure.d CalcDoublePrecision()
    Protected.i I, N
    Protected.d FPr, res
    
    ; IT CARRIES OUT AN APPROXIMATION OF MACHINE PRECISION

    FPr = 1.0  
    N = 2
    While (1.0 + FPR) > 1.0
      FPr = 0.5 * FPr
      
      ; The For Next is only to confuse the copiler optimation
      ; especally for the C-Backend. It's to force the compiler
      ; to copy the FPr variable to the memory. Because the internal 
      ; precision of floating Point registers might be higher.
      ; x32 and x64 (Floating Point calculation precision is 80 Bit)
      ; We need the precsion what can be saved in memory
;       For I = 1 To N
;         res = res + FPR * I
;       Next
    Wend
    res + 1 
    
    ProcedureReturn FPr
  EndProcedure
  
  Procedure.i EQd(V1.d, V2.d, Tol.d) ; As Boolean
         
    If Abs(V1 - V2) <= Tol
      ProcedureReturn #True
    EndIf
    
    ProcedureReturn #False
  EndProcedure

  
  
  Debug CalcDoublePrecision()
  Debug StrD(CalcDoublePrecision(),32)
  
  #Tol = 0.0000000000000002
;        0.0000000000000001110223025 ; 1.11e-16

  Define TOL.d = CalcDoublePrecision()
  
  Define x.d = 1
 ;Define y1.d = 1.0000000000000001
  Define y1.d = 1.0000000000000001
  Define y2.d = 1.0000000000000002
  Define y3.d = 1.0000000000000003
  Define y4.d = 1.0000000000000004
  Define y5.d = 1.0000000000000005
  
  Debug Bool(EQd(x,y1,TOL))
  Debug Bool(EQd(x,y2,TOL))
  Debug Bool(EQd(x,y3,TOL))
  Debug Bool(EQd(x,y4,TOL)) 
  Debug Bool(EQd(x,y5,TOL))

[STARGÅTE]

Die 2.2e-308 bezieht sich nicht auf die Toleranz oder Genauigkeit sonder auf die kleinste darstellbare Zahl.
Die Genauigkeit von Doubles ist etwa 1e-16. Das heißt die Zahlen haben nur 16 darstellbare Dezimalzahlen.

[DarkDragon]

Die Genauigkeiten sind dynamisch.

Wie beschreibt man natürliche Zahlen binär?

x = 2^0 * b_1 + 2^1 * b_2 + 2^2 * b_3 + ... 2^n * b_n

wobei b_k jeweils für ein Bit in der Binärdarstellung steht. Binär also b_1, b_2, ..., b_n.

Wie beschreibt man Fließkommazahlen binär?

x = 2^g * b_1 + 2^(g-1) * b_2 + 2^(g-2) * b_3 + ... 2^(g-n) * b_n

oder

x = 2^g * (2^0 * b_1 + 2^(-1) * b_2 + 2^(-2) * b_3 + ... 2^(-n) * b_n)

wobei die ersten 2^g die Verschiebung des Exponenten (und damit das Komma von Fließkommazahlen im binären) darstellt. Wenn man die 2^k auflöst steht da also

x = 2^g * (1 * b_1 + 0.5 * b_2 + 0.25 * b_3 + ...)

Bei 32bit floats sind 1 Bit für das Vorzeichen, 8 Bit für den Exponenten (g) und 23 Bits für den eigentlichen Wert (b_1, ... b_n) vorgesehen nach IEEE754.

Numerisch gesehen sind bei einigen Operationen die Fehler größer als bei anderen. Der relative Fehler (nicht der absolute) bei der Subtraktion zweier ungefähr gleich großer Zahlen ist z.B. unlimitiert. Kannst ja mal schauen wie du die miteinander verrechnen würdest in 32bit o.ä. und was die Schritte sind. Exponenten angleichen z.B. und du hast nicht unbegrenzt bits dafür. Aber es gibt keine absolute Genauigkeit bei Fließkommazahlen. Dazu müsste man Festkommazahlen verwenden.

[SMaag]

Die 2.2e-308 bezieht sich nicht auf die Toleranz oder Genauigkeit sonder auf die kleinste darstellbare Zahl.
Die Genauigkeit von Doubles ist etwa 1e-16. Das heißt die Zahlen haben nur 16 darstellbare Dezimalzahlen.

Ooohhh! Klar! Die Anzahl der Dezimalstellen haben ja nichts mit dem Exponenten zu tun, sondern hängen nur von der Mantisse ab!

Dann macht aber die Berechnung der Genauigkeit wenig Sinn. Die ändert sich ja je nach Expontenten. D.h. in etwa in der 16 Stelle liegt die Genauigkeit.
Diese hat je nach Exponenten einen anderen Wert.

Jetzt hab ich das mit der Genauigkeitsberechung aus professionellem FORTRAN CODE. Prinzipiell macht das aber keine Sinn, so eine Genauigkeit für Vergleich von Floats zu bestimmen, ohne den Exponenten zu berücksichtigen.
Ich hab aber auch schon Code gesehen, der einfach ein #Epsilon angibt mit 1e-99. Das macht dann aber auch wenig Sinn!

[STARGÅTE]

Ooohhh! Klar! Die Anzahl der Dezimalstellen haben ja nichts mit dem Exponenten zu tun, sondern hängen nur von der Mantisse ab!
Dann macht aber die Berechnung der Genauigkeit wenig Sinn. Die ändert sich ja je nach Expontenten. D.h. in etwa in der 16 Stelle liegt die Genauigkeit.
Diese hat je nach Exponenten einen anderen Wert.

Richtig. Genauer wäre es also von der Genauigkeit der Mantisse zu reden, die bei Floats und Doubles konstant ist. Hier kann man sich in PureBasic folgendes definieren:

Code: Alles auswählen

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

Das Epsilon ist hierbei die kleinstmögliche Zahl damit 1 + ε ≠ 1 ist.

Code: Alles auswählen

Debug Bool(1+#PB_DoubleEpsilon=1)
Debug Bool(1+#PB_DoubleEpsilon/2=1)

Ist die Zahl selbst großer, ist auch das entsprechende Epsilon größer damit sich die Zahl unterscheidet.
Wird z.B. eine Double Zahl großer als ~4.5e15, wird die Genauigkeit für Operationen schlechter als 1.0, sodass nicht alle ganzen Zahlen mehr darstellbar sind.
Auf der anderen Seite hat die Zahl 1.0e-200 natürlich ein viel kleineres "persönliches" Epsilon von etwa 2.0e-216.

[SMaag]

ich hab jetzt mal die Epsilons für verschiedene Werte durchprobiert. Bis 2^52.

Ich würde sagen, für den Vergleich ob 2 Floats identisch sind müsste ich mein Epsilon je nach Größe der Floats anapssen. Da wird der Vergleich dann aber übel! Dass man das anpassen müsste ist mir aber bisher noch nicht untergekommen. Aber eigentlich ist das klar. Ich frage mich aber, warum es dann keine allgemeine Lösung für das Problem gibt. Sowas wie ein IsEqualF(f1, f2)

Hier mal mein Code für die Liste der Epsilon-Werte 2^0 bis 2^52. Also wann wird eine 1 noch als 1 erkannt, wann eine 2, wann eine 4 als 4 usw.
Als Vergleichsschranke müsste man dann das doppelte Epsilon verwenden.

Code: Alles auswählen

 Procedure.d CalcSinglePrecision(x.f)
     
  Protected.f FPr
    
  ; IT CARRIES OUT AN APPROXIMATION OF MACHINE PRECISION
     
  If x = 0
    FPr = 1.0
    Debug "IsNull"
  Else   
    FPr = x  
  EndIf
  
  While (x + FPR) > x
    FPr = 0.5 * FPr
  Wend
    
  ProcedureReturn FPr ; * 2.0
EndProcedure
  
  Define I, N
  Define Epsilon.f
  
  Debug StrF(CalcSinglePrecision(N),25) + " : N = " + Str(0)
  ; Debug CalcSinglePrecision(1) ; N = 1 mit max Nachkommastellen ergibt 25 Nachkommastellen
  Debug " "
  
  Debug "Liste Epsilon(N)"
  Debug " "
  N = 1
  For I = 0 To 52
    Epsilon = CalcSinglePrecision(N)
    ; bei StrF muss man die Nachkommastellen angeben sonst wird mit 10 gearbeitet
    Debug Str(I) + " : " + StrF(Epsilon,25) + " : N = " + Str(N)
    N * 2
  Next

das kommt dabei raus

Liste Epsilon(N)

0 : 0.0000000000000001110223025 : N = 1
1 : 0.0000000000000002220446049 : N = 2
2 : 0.0000000000000004440892099 : N = 4
3 : 0.0000000000000008881784197 : N = 8
4 : 0.0000000000000017763568394 : N = 16
5 : 0.0000000000000035527136788 : N = 32
6 : 0.0000000000000071054273576 : N = 64
7 : 0.0000000000000142108547152 : N = 128
8 : 0.0000000000000284217094304 : N = 256
9 : 0.0000000000000568434188608 : N = 512
10 : 0.0000000000001136868377216 : N = 1024
11 : 0.0000000000002273736754432 : N = 2048
12 : 0.0000000000004547473508865 : N = 4096
13 : 0.0000000000009094947017729 : N = 8192
14 : 0.0000000000018189894035459 : N = 16384
15 : 0.0000000000036379788070917 : N = 32768
16 : 0.0000000000072759576141834 : N = 65536
17 : 0.0000000000145519152283669 : N = 131072
18 : 0.0000000000291038304567337 : N = 262144
19 : 0.0000000000582076609134674 : N = 524288
20 : 0.0000000001164153218269348 : N = 1048576
21 : 0.0000000002328306436538696 : N = 2097152
22 : 0.0000000004656612873077393 : N = 4194304
23 : 0.0000000009313225746154785 : N = 8388608
24 : 0.0000000018626451492309570 : N = 16777216
25 : 0.0000000037252902984619141 : N = 33554432
26 : 0.0000000074505805969238281 : N = 67108864
27 : 0.0000000149011611938476560 : N = 134217728
28 : 0.0000000298023223876953130 : N = 268435456
29 : 0.0000000596046447753906250 : N = 536870912
30 : 0.0000001192092895507812500 : N = 1073741824
31 : 0.0000002384185791015625000 : N = 2147483648
32 : 0.0000004768371582031250000 : N = 4294967296
33 : 0.0000009536743164062500000 : N = 8589934592
34 : 0.0000019073486328125000000 : N = 17179869184
35 : 0.0000038146972656250000000 : N = 34359738368
36 : 0.0000076293945312500000000 : N = 68719476736
37 : 0.0000152587890625000000000 : N = 137438953472
38 : 0.0000305175781250000000000 : N = 274877906944
39 : 0.0000610351562500000000000 : N = 549755813888
40 : 0.0001220703125000000000000 : N = 1099511627776
41 : 0.0002441406250000000000000 : N = 2199023255552
42 : 0.0004882812500000000000000 : N = 4398046511104
43 : 0.0009765625000000000000000 : N = 8796093022208
44 : 0.0019531250000000000000000 : N = 17592186044416
45 : 0.0039062500000000000000000 : N = 35184372088832
46 : 0.0078125000000000000000000 : N = 70368744177664
47 : 0.0156250000000000000000000 : N = 140737488355328
48 : 0.0312500000000000000000000 : N = 281474976710656
49 : 0.0625000000000000000000000 : N = 562949953421312
50 : 0.1250000000000000000000000 : N = 1125899906842624
51 : 0.2500000000000000000000000 : N = 2251799813685248
52 : 0.5000000000000000000000000 : N = 4503599627370496

[DarkDragon]

Ich würde sagen, für den Vergleich ob 2 Floats identisch sind müsste ich mein Epsilon je nach Größe der Floats anapssen. Da wird der Vergleich dann aber übel! Dass man das anpassen müsste ist mir aber bisher noch nicht untergekommen. Aber eigentlich ist das klar. Ich frage mich aber, warum es dann keine allgemeine Lösung für das Problem gibt. Sowas wie ein IsEqualF(f1, f2)

Die Erklärung dafür habe ich bereits geliefert. Warum sollte es bei einer dynamischen Präzision eine allgemeine Lösung geben? Das ginge nur wenn du NIEMALS anfängst zu rechnen, sondern Werte nur vergleichst, aber wann macht man das schon? Dann kannst du auch den ist gleich Operator verwenden. Sobald du mit den Zahlen irgendwas anderes machst als Vergleichen pflanzt sich ein Fehler fort, der ggf. alles eliminiert.

Hier mal ein kleines Beispiel, wo f2 (= 1.0) gar keinen Einfluss auf die Rechnung hat. Es danach mit dem Zielwert zu vergleichen, den man ggf. erwarten würde, bräuchte ein Epsilon, das entsprechend größer als 1.0 ist. Du kannst damit rumspielen und die Bits der Fließkommazahlen ändern, wie du Lust hast:

Code: Alles auswählen

EnableExplicit

Structure VarType
  StructureUnion
    f.f
    l.l
  EndStructureUnion
EndStructure

Procedure.f SetSign(f.f, bit.i)
  Protected *var.VarType = @f
  
  If Bool(bit)
    *var\l = *var\l | $80000000
  Else
    *var\l = *var\l & $80000000
  EndIf
  
  ProcedureReturn *var\f
EndProcedure

Procedure.f SetExponent(f.f, exponent.a)
  Protected *var.VarType = @f
  
  *var\l = (*var\l & (~$7F800000)) | (exponent << 23)
  
  ProcedureReturn *var\f
EndProcedure

Procedure.f SetMantissa(f.f, mantissa.i)
  Protected *var.VarType = @f
  
  *var\l = (*var\l & (~$007FFFFF)) | (mantissa & $007FFFFF)
  
  ProcedureReturn *var\f
EndProcedure

Procedure.s BinF(f.f)
  Protected *var.VarType = @f
  
  ProcedureReturn RSet(Bin(*var\l), 32, "0")
EndProcedure

Define f1.f = SetMantissa(SetExponent(SetSign(0.0, #False), 254), 0)
Define f2.f = SetMantissa(SetExponent(SetSign(0.0, #False), 127), 1)
Define f3.f = f1 - f2

Debug "f1: " + StrF(f1, 5)
Debug "f2: " + StrF(f2, 5)
Debug "f3: " + StrF(f3, 5)
Debug "f1 (binary): " + BinF(f1)
Debug "f2 (binary): " + BinF(f2)
Debug "f3 (binary): " + BinF(f3)

[STARGÅTE]

Ich würde sagen, für den Vergleich ob 2 Floats identisch sind müsste ich mein Epsilon je nach Größe der Floats anapssen. Da wird der Vergleich dann aber übel! Dass man das anpassen müsste ist mir aber bisher noch nicht untergekommen. Aber eigentlich ist das klar. Ich frage mich aber, warum es dann keine allgemeine Lösung für das Problem gibt. Sowas wie ein IsEqualF(f1, f2)

Bei verschiedenen Größen der Floats vergleicht man dann nicht absolut sondern relativ.
Eine Möglichkeit wäre z.B. Abs(x1-x2)/Abs(x1+x2) < epsilon zu vergleichen, dann reicht "ein" Epsilon.

Wenn du dir das Dividieren sparen willst (weil es langsam ist) kannst du die Exponenten und die Mantissen der beiden Zahlen auch direkt vergleichen und abschätzen, wann Zahlen "fast gleich" sind.
Ich habe mir dafür schon mal n Funktion geschrieben, müsste sich aber erst wieder suchen.

[DarkDragon]

Ich würde sagen, für den Vergleich ob 2 Floats identisch sind müsste ich mein Epsilon je nach Größe der Floats anapssen. Da wird der Vergleich dann aber übel! Dass man das anpassen müsste ist mir aber bisher noch nicht untergekommen. Aber eigentlich ist das klar. Ich frage mich aber, warum es dann keine allgemeine Lösung für das Problem gibt. Sowas wie ein IsEqualF(f1, f2)

Bei verschiedenen Größen der Floats vergleicht man dann nicht absolut sondern relativ.
Eine Möglichkeit wäre z.B. Abs(x1-x2)/Abs(x1+x2) < epsilon zu vergleichen, dann reicht "ein" Epsilon.

Ja, wobei das eigentlich ohne Fließkommazahlen berechnet werden müsste, oder nicht? Sonst haben wir bei der Berechnung des relativen Fehlers ggf. einen hohen relativen Fehler
Gerade x1-x2 macht mir Sorgen (Auslöschung). Dein zweiter Punkt würde mich interessieren, wie du das gelöst hast.

Gute Nacht erstmal.

[STARGÅTE]

Hier die Funktion aus meinem größeren Projekt:

Code: Alles auswählen

Procedure.i CompareDouble(Double1.d, Double2.d, Accuracy.i=46)
	
	Protected *Double1.Quad = @Double1
	Protected *Double2.Quad = @Double2
	Protected Sign.i, Exponent.q, Fraction.q
	Protected Uncertainty = 1<<(52-Accuracy) - 1
	
	Sign        = *Double2\q>>63 & 1    - *Double1\q>>63 & 1
	Exponent    = *Double2\q>>52 & $7FF - *Double1\q>>52 & $7FF
	
	If Sign <> 0 Or Exponent > 1 Or Exponent < -1
		; Die zweite Zahl ist deutlich kleiner oder deutlich größer
		If Double1 > Double2
			ProcedureReturn 1
		ElseIf Double1 < Double2
			ProcedureReturn -1
		EndIf
	ElseIf Exponent = -1
		; Die zweite Zahl ist etwas kleiner
		Sign     = 1 - 2 * (*Double2\q>>63 & 1)
		Fraction = ( (*Double2\q&$FFFFFFFFFFFFF) - (*Double1\q&$FFFFFFFFFFFFF)|1<<52 ) * Sign
		If Fraction < -Uncertainty
			ProcedureReturn 1
		ElseIf Fraction > Uncertainty
			ProcedureReturn -1
		EndIf
	ElseIf Exponent = 1
		; Die zweite Zahl ist etwas größer
		Sign     = 1 - 2 * (*Double2\q>>63 & 1)
		Fraction = ( (*Double2\q&$FFFFFFFFFFFFF)|1<<52 - (*Double1\q&$FFFFFFFFFFFFF) ) * Sign
		If Fraction < -Uncertainty
			ProcedureReturn 1
		ElseIf Fraction > Uncertainty
			ProcedureReturn -1
		EndIf
	Else
		; Die Zahlen sind fast gleich
		Sign     = 1 - 2 * (*Double2\q>>63 & 1)
		Fraction = ( *Double2\q&$FFFFFFFFFFFFF - *Double1\q&$FFFFFFFFFFFFF ) * Sign
		If Fraction < -Uncertainty
			ProcedureReturn 1
		ElseIf Fraction > Uncertainty
			ProcedureReturn -1
		EndIf
	EndIf
	
	ProcedureReturn 0
	
EndProcedure


Debug CompareDouble(1.0, 2.0)
Debug CompareDouble(3.0, 2.0)

Debug "---"

Debug CompareDouble(1.000000000000001,  0.999999999999999)  ; Ist in den höchsten 46 Bit gleich.
Debug          Bool(1.000000000000001 > 0.999999999999999)

Debug "---"

Debug CompareDouble(1.000000000001, 0.999999999999)  ; Ist in den höchsten 46 Bit größer.
Debug CompareDouble(1.000000000001, 0.999999999999, 35) ; Ist in den höchsten 35 Bit gleich.

Debug "---"

Debug CompareDouble(1.001, 1.002, 8)
Debug CompareDouble(1001, 1002, 8)
Debug CompareDouble(1001, 1003, 8)