Code-Beispiel
quadruple precisions module
| Lizenz: | Erster Autor: | Letzte Bearbeitung: |
| k. A. |  Volta |  21.10.2010 |
Dies sind einige Module die es erlauben in FreeBASIC mit einer Genauigkeit von 30 Stellen zu rechnen. Das ist eine Verdopplung gegenüber den Double-Variablen.
Ob schon alles fehlerfrei funktioniert kann ich nicht 100% sagen. Einige Testprogramme laufen aber einwandfrei.
Die Idee dazu ist nicht von mir, 
hier hat srvaldez diese Module aus Fortran übersetzt und für FB angepasst.
Eine echte Fleißarbeit. Die Umsetzung in FPU-Befehlen ist Ihm dann aber nicht so ganz gelungen, FBC setzt aus den Basic-Befehlen bessere Assembleranweisungen um.
Deshalb habe ich alles wieder in Basic-Befehlen zurück gewandelt. Leider produzierte das aber eine deutlich schlechtere "precision". :(
Diese Probleme hatte srvaldez wohl auch schon und hatte sie in der Funktion 'exactmul2' vermutet.
Nach allen Rechenregeln ist a1 = (a - t) + t auch a1 = a - t + t also a1 = a !!!
Nicht so für die FPU. Ich kann es immer noch nicht einsehen aber einige Tests zeigen es, 'a1 = (a - t) + t' ergibt deutlich bessere Ergebnisse. Wichtig war auch das Controlword der FPU auf &h27f zu setzen. Ok, also diese Funktion doch in Assembleranweisungen.
Auszug aus 'quad_1x.bi'
....
Function exactmul2(ByRef a As Double, ByRef c As Double) As quad 'Result(ac)
' Procedure exactmul2, translated from Pascal, from:
' Linnainmaa, Seppo (1981). Software for doubled-precision floating-point
' computations. ACM Trans. on Math. Software (TOMS), 7, 272-283.
Dim As quad ac
Dim As Double a1, a2, c1, c2, t
Dim As UShort ocw, ncw = &h27f
Asm
fstcw Word [ocw]
fldcw Word [ncw]
't = constant * a
mov eax, Dword Ptr [a]
fld qword Ptr [eax]
fmul qword Ptr [constant]
fstp qword Ptr [t]
'a1 = (a - t) + t ' Lahey's optimization removes the brackets
' ' and sets a1 = a which defeats the whole point.
'mov eax, dword ptr [a]
fld qword Ptr [eax]
fsub qword Ptr [t]
fadd qword Ptr [t]
fstp qword Ptr [a1]
'a2 = a - a1
'mov eax, dword ptr [a]
fld qword Ptr [eax]
fsub qword Ptr [a1]
fstp qword Ptr [a2]
't = constant * c
mov eax, Dword Ptr [c]
fld qword Ptr [eax]
fmul qword Ptr [constant]
fstp qword Ptr [t]
'c1 = (c - t) + t
'mov eax, dword ptr [c]
fld qword Ptr [eax]
fsub qword Ptr [t]
fadd qword Ptr [t]
fstp qword Ptr [c1]
'c2 = c - c1
'mov eax, dword ptr [c]
fld qword Ptr [eax]
fsub qword Ptr [c1]
fstp qword Ptr [c2]
'ac.hi = a * c
'mov eax, dword ptr [c]
fld qword Ptr [eax]
mov eax, Dword Ptr [a]
fmul qword Ptr [eax]
fstp qword Ptr [ac]
'ac.lo = (((a1 * c1 - ac.hi) + a1 * c2) + c1 * a2) + c2 * a2
fld qword Ptr [c1]
fmul qword Ptr [a1]
fsub qword Ptr [ac]
fld qword Ptr [c2]
fmul qword Ptr [a1]
fxch st(1)
faddp
fld qword Ptr [a2]
fmul qword Ptr [c1]
fxch st(1)
faddp
fld qword Ptr [a2]
fmul qword Ptr [c2]
fxch st(1)
faddp
fstp qword Ptr [ac+8]
fldcw Word [ocw]
End Asm
Return ac
End Function 'exactmul2
....Hier die Dateien:
t_cubert.bas
t_logexp.bas
t_cst.bas
test_quad.bas
quad_1x.bi
quad_2.bi
quad_3.bi
quad_4.bi
| Zusätzliche Informationen und Funktionen |
|---|
- Das Code-Beispiel wurde am 30.09.2010 von Volta angelegt.
- Die aktuellste Version wurde am 21.10.2010 von Volta gespeichert.
|
|
Wer ist online?
Buchempfehlung
FreeBASIC-Chat
FreeBASIC-Nachrichten jetzt auch über Twitter erhalten. Follow us!
Bearbeiten
Versionen