\ Doubled-precision floating point arithmetic \ \ --------------------------------------------------- \ (c) Copyright 2006 Julian V. Noble. \ \ Permission is granted by the author to \ \ use this software for any application pro- \ \ vided this copyright notice is preserved. \ \ --------------------------------------------------- \ This is an ANS Forth program requiring the \ FLOAT, FLOAT EXT, FILE and TOOLS EXT wordsets. \ \ Environmental dependences: \ Assumes independent floating point stack \ FPU must be set for IEEE 64-bit internal arithmetic \ ^^^^ \ based on "Software for Doubled-Precision Floating-Point Computations" \ by Seppo Linnainmaa \ ACM Transactions on Mathematical Software, \ Vol 7, No 3, September 1981, Pages 272-283 marker -double FALSE [IF] \ use these to compile the versions using include ftran202.f \ the FORmula TRANslator and Baden's flocals include flocals.f [THEN] : fvariables: 0 DO FVARIABLE LOOP ; [undefined] fvalue [IF] : fvalue CREATE HERE F! 1 DFLOATS allot DOES> F@ ; [THEN] [undefined] fto [IF] : fto ' >BODY STATE @ IF POSTPONE FLITERAL POSTPONE F! ELSE F! THEN ; IMMEDIATE [THEN] [undefined] f-rot [IF] : f-rot FROT FROT ; [THEN] [undefined] ftuck [IF] : ftuck FSWAP FOVER ; [THEN] [undefined] f>s [IF] : f>s f>d d>s ; [THEN] [undefined] s>f [IF] : s>f s>d d>f ; [THEN] \ ---------------------------------------- LOAD, STORE : r128@ DUP F@ FLOAT+ F@ ; : r128! DUP FLOAT+ F! F! ; \ ------------------------------------ END LOAD, STORE \ ----------------------------------- data types ---- : real*16 \ create a double-double variable CREATE 2 FLOATS ALLOT ; : ddconstant CREATE HERE r128! 2 FLOATS ALLOT DOES> r128@ ; \ --------------------------------------------------- 3.1415926535897931e0 1.2246467991473532e-16 ddconstant ddpi \ based on "Software for Doubled-Precision Floating-Point Computations" \ by Seppo Linnainmaa \ ACM Transactions on Mathematical Software, \ Vol 7, No 3, September 1981, Pages 272-283 : fvariables: 0 DO FVARIABLE LOOP ; \ determine base and precision of fpu \ (( 4e0 3e0 F/ 1e0 F- 3e0 F* 1e0 F- FVALUE u u 2e0 F/ 1e0 F+ 1e0 F- FVALUE r r F0= 0= [IF] r FTO u [THEN] 2e0 3e0 F/ 0.5e0 F- 3e0 F* 0.5e0 F- FVALUE uu uu 2e0 F/ 0.5e0 F+ 0.5e0 F- FVALUE rr rr F0= 0= [IF] rr FTO uu [THEN] u uu F/ ( f: -- beta) uu FLN FOVER FLN F/ FNEGATE 0.5e0 F+ F>S F>S CR CR .( base = ) . .( precision = ) . \ FORGET u \ )) \ ---------------- Exact multiplication ------------- 134217729 S>F FCONSTANT split 9 fvariables: t a1 a2 b1 b2 b21 b22 ta tb 8 fvariables: q qq x xx y yy z1 z2 2 fvariables: aaa bbb FALSE [IF] \ these versions use ftran202.f and flocals.f : exactmul ( f: a b -- x xx) \ multiply 2 fp#'s to get ddfp# FRAME| bb aa | f" t = bb*[split]" f" a1 = (bb-t) + t" f" a2 = bb - a1" f" t = aa*[split]" f" b1 = (aa-t) + t" f" b2 = aa - b1" f" t = b2*[split]" f" b21 = (b2-t) + t" f" b22 = b2 - b21" f" bb*aa" FDUP f" t=" ( f: x) f" ((((a1*b1 - t)+a1*b2)+b1*a2)+b21*a2)+b22*a2" |FRAME ; ( f: x xx) : dd/ ( f: x xx y yy -- [x+xx]/[y+yy] ) f" yy=" f" y=" f" xx=" f" x=" f" abs(y)" F0= ABORT" Can't divide by 0!" f" z1=x/y" f" exactmul(y,z1)" f" qq=" f" q=" f" z2=((((x-q)-qq)+xx)-z1*yy)/(y+yy)" f" q=z1+z2" f" q" f" (z1-q)+z2" ; : dd* ( f: x xx y yy -- [x+xx]*[y+yy] ) f" yy=" f" y=" f" xx=" f" x=" f" exactmul(x,y)" f" qq=" f" z1=" f" z2=((x+xx)*yy+xx*y)+qq" f" q=z1+z2" f" q" f" (z1-q)+z2" ; : dd+ ( f: x xx y yy -- [x+xx] + [y+yy] ) f" yy=" f" y=" f" xx=" f" x=" f" z1=x+y" f" q=x-z1" f" z2=(((q+y)+(x-(q+z1)))+xx)+yy" f" q=z1+z2" f" q" f" (z1-q)+z2" ; \ Square root based on T.J. Dekker, \ "A Floating-Point Technique for Extending the Available Precision" \ Numerische Mathematik 18 (1971) 224-242. : ddsqrt ( f: x xx -- ddsqrt[x+xx]) f" xx=" f" x=" f" x" F0< ABORT" Can't take sqrt of negative number!" f" q = sqrt(x)" f" exactmul(q,q)" f" yy=" f" y=" f" qq = (((x-y)-yy)+xx)*0.5/q" f" y=q + qq" f" y" f" (q-y)+qq" ; [THEN] : exactmul ( f: a b -- x xx) \ multiply 2 fp#'s to get ddfp# aaa F! bbb F! bbb F@ split F* t F! t F@ bbb F@ FOVER F- F+ FDUP a1 F! ( f: a1) FNEGATE bbb F@ F+ a2 F! aaa F@ FDUP split F* t F! ( f: aaa) t F@ ftuck F- F+ b1 F! aaa F@ b1 F@ F- FDUP FDUP b2 F! ( f: b2 b2) split F* t F! t F@ ftuck F- F+ FDUP b21 F! ( f: b21) FNEGATE b2 F@ F+ b22 F! bbb F@ aaa F@ F* FDUP t F! a1 F@ b1 F@ F* t F@ F- a1 F@ b2 F@ F* F+ b1 F@ a2 F@ F* F+ b21 F@ a2 F@ F* F+ b22 F@ a2 F@ F* F+ ; : dd/ ( f: x xx y yy -- [x+xx]/[y+yy] ) yy F! y F! xx F! x F! y F@ FABS F0= ABORT" Can't divide by 0!" x F@ y F@ F/ FDUP z1 F! y F@ exactmul qq F! FDUP q F! ( f: q) FNEGATE x F@ F+ qq F@ F- xx F@ F+ z1 F@ yy F@ F* F- y F@ yy F@ F+ F/ FDUP z2 F! z1 F@ F+ FDUP FNEGATE z1 F@ F+ z2 F@ F+ ; : dd* ( f: x xx y yy -- [x+xx]*[y+yy] ) yy F! y F! xx F! x F! x F@ y F@ exactmul qq F! z1 F! x F@ xx F@ F+ yy F@ F* xx F@ y F@ F* F+ qq F@ F+ FDUP z2 F! z1 F@ F+ FDUP FNEGATE z1 F@ F+ z2 F@ F+ ; : dd+ ( f: x xx y yy -- [x+xx] + [y+yy] ) yy F! y F! xx F! x F! x F@ y F@ F+ z1 F! x F@ z1 F@ F- FDUP q F! y F@ F+ ( f: q+y) x F@ q F@ z1 F@ F+ F- ( f: q+y x-[q+z1]) F+ xx F@ F+ yy F@ F+ FDUP z2 F! z1 F@ F+ FDUP ( f: z1+z2 z1+z2) FNEGATE z1 F@ F+ z2 F@ F+ ; : ddnegate ( f: x xx -- -x -xx) FNEGATE FSWAP FNEGATE FSWAP ; : dd- ( f: x xx y yy -- [x+xx] - [y+yy] ) ddnegate dd+ ; \ ddsqrt based on T.J. Dekker, \ "A Floating-Point Technique for Extending the Available Precision" \ Numerische Mathematik 18 (1971) 224-242. : ddsqrt ( f: x xx -- ddsqrt[x+xx]) xx F! FDUP x F! F0< ABORT" Can't take sqrt of negative number!" x F@ FSQRT FDUP q F! FDUP exactmul yy F! FDUP y F! FNEGATE x F@ F+ yy F@ F- xx F@ F+ 0.5e0 F* q F@ F/ FDUP qq F! q F@ F+ FDUP FNEGATE q F@ F+ qq F@ F+ ; 1e0 0e0 ddconstant dd=1 \ ----------------------------------- stack ops ----- fvariable dtemp real*16 ddtemp real*16 ddtemp1 : ddswap ( f: x xx y yy -- y yy x xx ) dtemp F! F-ROT dtemp F@ F-ROT ; : dddup FOVER FOVER ; : dddrop FDROP FDROP ; : ddover ddtemp r128! dddup ddtemp1 r128! ddtemp r128@ ddtemp1 r128@ ; : ddtuck ddtemp r128! ddtemp1 r128! ddtemp r128@ ddtemp1 r128@ ddtemp r128@ ; \ --------------------------------------------------- : dd^2 dddup dd* ; \ : dd^n ( n -- ) ( f: x xx -- [x+xx]^n ) \ raise dd to positive or negative integer power \ return 1 if n=0, dd^{-|n|} if n<0 \ dd=1 ddswap ( f: 1e0 0e0 x xx ) \ DUP 0= IF dddrop dd=1 drop EXIT THEN \ DUP 0< SWAP ABS ( -- sign |n| ) \ BEGIN DUP 0> WHILE \ DUP 1 AND IF ddtuck dd* ddswap THEN z^2 \ 2/ \ REPEAT dddrop DROP \ IF dd=1 ddswap dd/ THEN \ ; \ Tests for DoubleDouble \ February 11th, 2006 MARKER -ddtest 15 set-precision 1e0 0e0 3e0 0e0 dd/ ddconstant dd1/3 1e0 0e0 6e0 0e0 dd/ ddconstant dd1/6 dd1/3 dd1/6 dd+ FSWAP CR CR .( Addition: 1/3 + 1/6) CR 5 spaces .( 5.00000000000000E-1 -1.54074395550979E-33 <- should be this) CR 5 spaces FS. FS. CR dd1/3 dd1/6 dd- FSWAP CR .( Subtraction: 1/3 - 1/6) CR 5 spaces .( 1.66666666666667E-1 9.25185853854297E-18 <- should be this) CR 5 spaces FS. FS. CR 6e0 0e0 dd1/6 dd* FSWAP CR .( Multiplication: 6 * 1/6) CR 5 spaces .( 1.00000000000000E0 0.00000000000000E-1 <- should be this) CR 5 spaces FS. FS. CR .( 6 * 1/3) 6e0 0e0 dd1/3 dd* FSWAP CR 5 spaces .( 2.00000000000000E0 0.00000000000000E-1 <- should be this) CR 5 spaces fs. fs. CR dd1/3 dd1/6 dd/ FSWAP CR .( Division: [1/3] / [1/6] ) CR 5 spaces .( 2.00000000000000E0 0.00000000000000E-1 <- should be this) CR 5 spaces FS. FS. CR CR .( [1/6] / [1/3]) dd1/6 dd1/3 dd/ FSWAP CR 5 spaces .( 5.00000000000000E-1 0.00000000000000E-1 <- should be this) CR 5 spaces FS. FS. CR .( Square root: sqrt[1/324] ) dd1/6 dd1/3 dd* FOVER FOVER dd* ddsqrt FOVER FOVER FSWAP CR 5 spaces FS. FS. .( <- sqrt[1/324] = 1/18) CR 5 spaces .( 5.55555555555556E-2 3.08395284618099E-18 <- should be this) CR 5 spaces .( multiply by 18) 18e0 0e0 dd* FSWAP CR 5 spaces .( 1.00000000000000E0 0.00000000000000E-1 <- should be this) CR 5 spaces FS. FS.