(* * LANGUAGE : ANS Forth with extensions * PROJECT : Forth Environments * DESCRIPTION : Doubled-precision floating point arithmetic * CATEGORY : Tools * AUTHOR : Marcel Hendrix * AUTHOR : (c) Copyright 2006 Julian V. Noble. * : Permission is granted by the author to use this software for * : any application provided this copyright notice is preserved. * LAST CHANGE : Monday, February 20, 2006, 23:03 PM, Marcel Hendrix, added D. Bailey's functions * LAST CHANGE : February 12, 2006, Marcel Hendrix *) NEEDS -cplx_fsl REVISION -dd-fp "ÄÄÄ 128-bit FP words Version 1.05 ÄÄÄ" PRIVATES DOC (* Environmental dependencies: FPU *must* be set for IEEE 64-bit internal arithmetic Based on "Software for Doubled-Precision Floating-Point Computations", Seppo Linnainmaa, ACM Transactions on Mathematical Software, Vol 7, No 3, September 1981, Pages 272-283 Note: When a number is >= 1e23 (more than 53 bits), it must be split in high and low parts. *) ENDDOC 53-bits! NESTING @ 1 = [IF] CR .( Warning: 53-bit FP mode set.) CR .( Note: this mode is reset after an abort!) [THEN] 4e 3e F/ F1- 3e F* F1- FVALUE w w F2/ F1+ F1- FVALUE r r F0<> [IF] r TO w [THEN] 2e 3e F/ 0.5e F- 3e F* 0.5e F- FVALUE uu uu F2/ 0.5e F+ 0.5e F- FVALUE rr rr F0<> [IF] rr TO uu [THEN] w uu F/ ( F: -- beta) uu FABS FLN FOVER FABS FLN F/ FNEGATE 0.5e F+ NESTING @ 1 = [IF] CR .( base = ) FSWAP F>S . CR .( precision = ) F>S . [ELSE] F2DROP [THEN] FORGET w -- Basic fetching and storing : DD@ ( addr -- ) ( F: -- r1 r2 ) DF2@ ; : DD! ( addr -- ) ( F: r1 r2 -- ) DF2! ; : DD@+ ( addr1 -- addr2 ) ( F: -- r1 r2 ) DF2@+ ; : DD!+ ( addr1 -- addr2 ) ( F: r1 r2 -- ) DF2!+ ; -- We need these for the methods array in fsl_util.frt : _DD! ( addr -- ) ( F: r1 r2 -- ) DD! ; : _DD@ ( addr -- ) ( F: -- r1 r2 ) DD@ ; 2 DFLOATS =: =DDFLOAT : DDFLOATS ( n1 -- n2 ) =DDFLOAT * ; -- Create a double-double variable : DDVARIABLE CREATE ( -- ) =DDFLOAT ALLOT DOES> ( -- addr ) ; : DDF, ( F: ddouble -- ) DF2, ; -- Create a double-double constant : DDCONSTANT CREATE ( F: r1 r2 -- ) DDF, DOES> ( F: -- r1 r2 ) DD@ ; 1 [IF] : DDCOPY ( a1 incx a2 incy n -- ) >S 2>R 2>R S> 2R> 2R> 5 BLAS-LIB: izcopy FOREIGN DROP ; [ELSE] : DDCOPY ( a1 incx a2 incy n -- ) 0 0 =DDFLOAT LOCALS| sz incx*sz incy*sz n incy a2 incx a1 | incx 1 = incy 1 = AND IF a1 a2 n sz * MOVE EXIT ENDIF incx sz * TO incx*sz incy sz * TO incy*sz incx 1 = IF a1 a2 n incy*sz * BOUNDS ?DO DD@+ I DD! incy*sz +LOOP DROP EXIT ENDIF incy 1 = IF a2 a1 n incx*sz * BOUNDS ?DO I DD@ DD!+ incx*sz +LOOP DROP EXIT ENDIF a1 a2 n incy*sz * BOUNDS ?DO DUP DD@ I DD! incx*sz + incy*sz +LOOP DROP ; [THEN] 80-bits! ( turn it off so we can compile DDCONSTANT accurately) 0e 0e DDCONSTANT dd=0 1e 0e DDCONSTANT dd=1 10e 0e DDCONSTANT dd=10 6.283185307179586232e+00 2.449293598294706414e-16 DDCONSTANT dd=PI*2 3.141592653589793116e+00 1.224646799147353207e-16 DDCONSTANT dd=PI 1.570796326794896558e+00 6.123233995736766036e-17 DDCONSTANT dd=PI/2 7.853981633974482790e-01 3.061616997868383018e-17 DDCONSTANT dd=PI/4 1.963495408493620697e-01 7.654042494670957545e-18 DDCONSTANT dd=PI/16 2.356194490192344837e+00 9.1848509936051484375e-17 DDCONSTANT dd=3*PI/4 2.718281828459045091e+00 1.445646891729250158e-16 DDCONSTANT dd=E 6.931471805599452862e-01 2.319046813846299558e-17 DDCONSTANT dd=LN2 2.302585092994045901e+00 -2.170756223382249351e-16 DDCONSTANT dd=LOG10 1.23259516440783e-32 FCONSTANT ddeps -- 2^-106 -- Table of sin(k * pi/16) and cos(k * pi/16) -- CREATE sin_table PRIVATE 1.950903220161282758e-01 DF, -7.991079068461731263e-18 DF, 3.826834323650897818e-01 DF, -1.005077269646158761e-17 DF, 5.555702330196021776e-01 DF, 4.709410940561676821e-17 DF, $667F3BCD,3FE6A09E 2, $13B26456,BC8BDD34 2, CREATE cos_table PRIVATE 9.807852804032304306e-01 DF, 1.854693999782500573e-17 DF, 9.238795325112867385e-01 DF, 1.764504708433667706e-17 DF, 8.314696123025452357e-01 DF, 1.407385698472802389e-18 DF, $667F3BCD,3FE6A09E 2, $13B26456,BC8BDD34 2, 134217729e FCONSTANT split PRIVATE -- For exactmul 0e FVALUE a1 PRIVATE 0e FVALUE a2 PRIVATE 0e FVALUE a3 0e FVALUE b1 PRIVATE 0e FVALUE b2 PRIVATE 0e FVALUE aaa PRIVATE 0e FVALUE bbb PRIVATE 0e FVALUE q PRIVATE 0e FVALUE qq PRIVATE 0e FVALUE x PRIVATE 0e FVALUE xx PRIVATE 0e FVALUE y PRIVATE 0e FVALUE yy PRIVATE 0e FVALUE z1 PRIVATE 0e FVALUE z2 PRIVATE DDVARIABLE 'pp PRIVATE DDVARIABLE 'rr PRIVATE 53-bits! -- Designed to round through memory. : quick_two_sum ( 'err -- ) ( F: a b -- s ) F2DUP F+ FLOCAL s ( F: -- a b ) s FROT F- F- DF! s ; PRIVATE -- Multiply 2 fp#'s to get ddfp# : exactmul ( F: a b -- x xx ) TO bbb TO aaa bbb split F* ( t) bbb FOVER F- F+ TO a1 bbb a1 F- TO a2 aaa split F* ( t) aaa FOVER F- F+ TO b1 aaa b1 F- TO b2 b2 split F* ( t) b2 FOVER F- F+ b2 FOVER F- F+ a2 F* TO a3 bbb aaa F* FDUP ( F: -- x x ) a1 b1 F* FSWAP F- a1 b2 F* F+ b1 a2 F* F+ a3 F+ ; -- Multiply 2 ddfp#'s to get ddfp# : DD* ( F: x xx y yy -- [x+xx] * [y+yy] ) TO yy TO y TO xx TO x x y exactmul TO qq TO z1 x xx F+ yy F* xx y F* F+ qq F+ TO z2 z1 z2 F+ FDUP FNEGATE z1 F+ z2 F+ ; : DDF* ( F: x xx y -- [x+xx] * y ) TO y TO xx ( x) y exactmul TO qq TO z1 xx y F* qq F+ TO z2 z1 z2 F+ FDUP FNEGATE z1 F+ z2 F+ ; -- Add 2 ddfp#'s to get ddfp# : DD+ ( F: x xx y yy -- [x+xx] + [y+yy] ) TO yy TO y TO xx TO x x y F+ FDUP TO z1 ( F: -- z1 ) x FOVER F- TO q ( F: -- z1 ) y q F+ ( F: -- z1 q+y ) FSWAP q F+ x FSWAP F- F+ ( F: -- q+y+x-[q+z1] ) xx F+ yy F+ ( F: -- z2 ) z1 FOVER F+ FDUP FNEGATE ( F: -- z2 z1+z2 -[z1+z2] ) z1 F+ FROT F+ ; : DD+!+ ( addr1 -- addr2 ) ( F: rr r -- ) DUP DD@ DD+ DD!+ ; -- Add ddfp# to fp# to get ddfp# : DDF+ ( F: x xx y -- [x+xx] + y ) TO y TO xx TO x x y F+ TO z1 x z1 F- TO q y q F+ ( F: -- q+y ) x q z1 F+ F- F+ ( F: -- q+y+x-[q+z1] ) xx F+ ( F: -- z2 ) z1 FOVER F+ FDUP FNEGATE z1 F+ FROT F+ ; ( F: -- z2 z1+z2 -[z1+z2] ) -- Negate ddfp# : DDNEGATE ( F: x xx -- -x -xx ) FNEGATE FSWAP FNEGATE FSWAP ; : DDABS ( F: x xx -- |x+xx| ) FOVER F0< IF ddnegate ENDIF ; -- Subtract 2 ddfp#'s to get ddfp# : DD- ( F: x xx y yy -- [x+xx] - [y+yy] ) ddnegate dd+ ; : DD/ ( F: x xx y yy -- c cc ) TO yy TO y TO xx TO x x y F/ TO z1 z1 y exactmul TO qq TO q x q F- qq F- xx F+ z1 yy F* F- y yy F+ F/ TO z2 z1 z2 F+ FDUP FNEGATE z1 F+ z2 F+ ; -- Divide ddfp# by fp# to get ddfp# : DDF/ ( F: x xx y -- [x+xx]/y ) TO y TO xx TO x x y F/ TO z1 z1 y exactmul TO qq TO q x q F- qq F- xx F+ y F/ TO z2 z1 z2 F+ FDUP FNEGATE z1 F+ z2 F+ ; -- 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] ) TO xx TO x x FSQRT TO q q q exactmul TO yy TO y x y F- yy F- xx F+ F2/ q F/ TO qq qq q F+ FDUP FNEGATE q F+ qq F+ ; -- ----------------------------------------------------- : DD+! ( addr -- ) ( F: xx x -- ) DUP DD@ DD+ DD! ; : DD-! ( addr -- ) ( F: xx x -- ) DDNEGATE DD+! ; : DD0! ( addr -- ) 0e 0e DD! ; : DDVALUE ( F: xx x "str" -- ) CREATE IMMEDIATE FALIGN HERE 2 DFLOATS ALLOT DD! DOES> FALIGNED %VAR @ 0 %VAR ! CASE ( -to) -2 OF ALITERAL EVAL" DD-! " ENDOF ( +to) -1 OF ALITERAL EVAL" DD+! " ENDOF ( from) 0 OF ALITERAL EVAL" DD@ " ENDOF ( to) 1 OF ALITERAL EVAL" DD! " ENDOF ( 0to) 2 OF ALITERAL EVAL" DD0! " ENDOF ( 'of) 3 OF ALITERAL ENDOF ( sizeof) 4 OF DROP 2 DFLOATS ILITERAL ENDOF ( /of) 5 OF DROP 1 ILITERAL ENDOF DUP #-21 ?ERROR" DDVALUE: undefined message" ENDCASE ; -- ----------------------------------- stack ops ------- : DDSWAP ( F: x xx y yy -- y yy x xx ) F2SWAP ; : DDDUP ( F: x xx -- x xx x xx ) F2DUP ; : DDDROP ( F: x xx -- ) F2DROP ; : DDNIP ( F: dd1 dd2 -- dd2 ) F2SWAP F2DROP ; : DDOVER ( F: x xx y yy -- x xx y yy x xx ) F2OVER ; : DD2DUP ( F: x xx y yy -- x xx y yy x xx y yy ) F2OVER F2OVER ; : DDTUCK ( F: x xx y yy -- y yy x xx y yy ) F2SWAP F2OVER ; -- ----------------------------------------------------- : DD^2 ( F: x xx -- y yy ) dddup dd* ; : 1/DD ( F: x xx -- y yy ) dd=1 ddswap dd/ ; : DD1+ ( F: x xx -- y yy ) 1e ddf+ ; : DD1- ( F: x xx -- y yy ) -1e ddf+ ; : DD2/ ( F: x xx -- y yy ) 2e ddf/ ; : DD2* ( F: x xx -- y yy ) 2e ddf* ; -- Raise dd to positive or negative integer power -- Return 1 if n=0, dd^{-|n|} if n<0 : DD^n ( n -- ) ( F: x xx -- [x+xx]^n ) DUP #64 = IF DROP dd^2 dd^2 dd^2 dd^2 dd^2 dd^2 EXIT ENDIF DUP 0= IF DROP dddrop dd=1 EXIT ENDIF DUP 1 = IF DROP EXIT ENDIF DUP 2 = IF DROP dd^2 EXIT ENDIF DUP 3 = IF DROP dddup dd^2 dd* EXIT ENDIF DUP 4 = IF DROP dd^2 dd^2 EXIT ENDIF dd=1 ddswap DUP 0< SWAP ABS ( F: -- 1e 0e ) ( -- sign |n| ) BEGIN DUP 0> WHILE DUP 1 AND IF ddtuck dd* ddswap ENDIF dd^2 2/ REPEAT dddrop DROP IF 1/dd ENDIF ; -- Nearest integer : DDNINT ( F: x xx -- rr r ) FOVER FROUND 0e FLOCALS| lo hi a.lo a.hi | a.hi hi F= IF a.lo FROUND TO lo ( renormalize ) hi lo 'OF lo quick_two_sum ( F: -- hi ) lo EXIT ENDIF hi a.hi F- FABS 0.5e F= a.lo F0< \ Had to tweak it back for unknown reason ??? AND IF hi F1- FDUP F0< IF F1+ ENDIF lo EXIT ENDIF hi lo ; PRIVATE : DDREM ( F: aa a bb b -- cc c ) ddover ddover dd/ ddnint dd* dd- ; PRIVATE : DIVREM ( F: aa a bb b -- cc c ) ( -- div ) ddover ddover dd/ ddnint FOVER F>S dd* dd- ; PRIVATE DOC (* Exponential. Computes exp(x) in double-double precision. We first reduce the size of x by noting that exp(kr + m) = exp(m) * exp(r)^k Thus by choosing m to be a multiple of log(2) closest to x, we can make |kr| <= log(2) / 2 = 0.3466. Now we can set k = 64, so that |r| <= 0.000542. Then exp(x) = exp(kr + s log 2) = (2^s) * [exp(r)]^64 Then exp(r) is evaluated using the familiar Taylor series. Reducing the argument substantially speeds up the convergence. *) ENDDOC 0e FCONSTANT F>DD ( F: xx -- xx x ) : DD>F ( F: xx x -- xx ) FDROP ; : DD>S ( F: xx x -- ) ( -- n ) FDROP F>S ; : pp*rr->pp ( -- ) 'pp DF@ 'rr DF@ exactmul TO qq TO z1 'pp DD@ F+ 'rr DFLOAT+ DF@ F* 'pp DFLOAT+ DF@ 'rr DF@ F* F+ qq F+ TO z2 z1 z2 F+ FDUP FNEGATE z1 F+ z2 F+ 'pp DD! ; PRIVATE : DDEXP ( F: a.hi a.lo -- r.hi r.lo ) FOVER -709e F<= IF F2DROP dd=0 EXIT ENDIF FOVER 709e F>= IF F2DROP +INF +INF EXIT ENDIF F2DUP F0= F0= AND IF F2DROP dd=1 EXIT ENDIF 2e 2e FLOCALS| m f | 0 LOCAL z dddup dd=LN2 dd/ ddnint DD>S TO z ( a) dd=LN2 z S>F ddf* dd- 64e ddf/ 'rr DD! 'rr DD@ dd^2 'pp DD! 'rr DD@ dd1+ 'pp DD@ f ddf/ ( F: -- accu tt ) #16 >R BEGIN dd+ ( F: -- accu+tt ) pp*rr->pp 1e +TO m f m f* TO f 'pp DD@ f ddf/ dddup ddabs DD>F 1e-32 F<= -1 R+! R@ 0< OR UNTIL -R dd+ ( F: -- accu+tt ) #64 dd^n z S>F F2^x ddf* ; DOC (* Logarithm. Computes ln(x) in double-double precision. The Taylor series for ln converges much more slowly than that of exp, due to the lack of the factorial term in the denominator. Hence this routine instead tries to determine the root of the function f(x) = exp(x) - a using Newton iteration. The iteration is given by x' = x - f(x)/f'(x) = x - (1 - a * exp(-x)) = x + a * exp(-x) - 1. Only one iteration is needed, since Newton's iteration approximately doubles the number of digits per iteration. *) ENDDOC : DDLN ( F: a.hi a.lo -- rhi rlo ) DDDUP F0= 1e F= AND IF FNIP 0e EXIT ENDIF dd=0 FLOCALS| x.lo x.hi a.lo a.hi | a.hi FLN TO x.hi \ Initial approximation x.hi x.lo ddnegate ddexp a.hi a.lo dd* x.hi x.lo dd+ dd1- ; : DDLOG ( F: a.hi a.lo -- rhi rlo ) DDLN dd=log10 dd/ ; DOC (* Computes sin(a) and cos(a) using Taylor series. Assumes |a| <= pi/32. *) ENDDOC : sincos_taylor ( 'sin 'cos -- ) ( F: a aa -- ) dddup F0= F0= AND IF dd=1 DD! DD! EXIT ENDIF FOVER FABS 1e-35 F* FLOCAL thresh ZLOCAL a dd=0 ZLOCAL t \ Term being added. a ZLOCAL p \ Current power of a. 1e FLOCAL f \ Denominator. a dd^2 ddnegate ZLOCAL x \ -sqr(a) 1e FLOCAL m a \ Current partial sum. #32 >R BEGIN p x dd* TO p 2e +TO m m m F1- F* f F* TO f p f ddf/ TO t t dd+ t DD>F FABS thresh F< -1 R+! R@ 0< OR UNTIL -R dddup SWAP DD! dd=1 ddswap dd^2 dd- ddsqrt DD! ; PRIVATE DOC (* To compute sin(x), we choose integers a, b so that x = s + a * (pi/2) + b * (pi/16) and |s| <= pi/32. Using the fact that sin(pi/16) = 0.5 * sqrt(2 - sqrt(2 + sqrt(2))) we can compute sin(x) from sin(s), cos(s). This greatly increases the convergence of the sine Taylor series. *) ENDDOC : DDSINCOS ( F: aa a -- ss s cc c ) dddup F0= F0= AND IF dd=1 EXIT ENDIF dd=pi*2 ddrem ZLOCAL r ( First reduce modulo 2*pi so that |r| <= pi ) dd=0 ZLOCAL t dd=0 ZLOCAL sin_t dd=0 ZLOCAL cos_t dd=0 ZLOCAL s dd=0 ZLOCAL c dd=0 ZLOCAL u dd=0 ZLOCAL v 0 0 0 0 LOCALS| abs_j j abs_k k | \ Now reduce by modulo pi/2 and then by pi/16 so that we obtain numbers a, b, and t r dd=pi/2 divrem TO t DUP ABS TO abs_j TO j t dd=pi/16 divrem TO t DUP ABS TO abs_k TO k abs_j 2 > IF +INF +INF +INF +INF EXIT ENDIF abs_k 4 > IF +INF +INF +INF +INF EXIT ENDIF t 'OF sin_t 'OF cos_t sincos_taylor abs_k 0= IF sin_t TO s cos_t TO c ELSE cos_table abs_k 1- COMPLEX[] DD@ TO u sin_table abs_k 1- COMPLEX[] DD@ TO v k 0> IF u sin_t dd* v cos_t dd* dd+ TO s u cos_t dd* v sin_t dd* dd- TO c ELSE u sin_t dd* v cos_t dd* dd- TO s u cos_t dd* v sin_t dd* dd+ TO c ENDIF ENDIF abs_j 0= IF s c ELSE j 1 = IF c s ddnegate ELSE j -1 = IF c ddnegate s ELSE s ddnegate c ddnegate ENDIF ENDIF ENDIF ; : DDSIN ( F: aa a -- ss s ) DDSINCOS dddrop ; : DDCOS ( F: aa a -- cc c ) DDSINCOS ddswap dddrop ; : DDTAN ( F: aa a -- tt t ) DDSINCOS dd/ ; DOC (* Instead of using Taylor series to compute arctan, we instead use Newton's iteration to solve the equation sin(z) = y/r or cos(z) = x/r where r = sqrt(x^2 + y^2). The iteration is given by z' = z + (y - sin(z)) / cos(z) (for equation 1) z' = z - (x - cos(z)) / sin(z) (for equation 2) Here, x and y are normalized so that x^2 + y^2 = 1. If |x| > |y|, then first iteration is used since the denominator is larger. Otherwise, the second is used. *) ENDDOC : DDATAN2 ( F: yy y xx x -- tt t ) ZLOCAL x ZLOCAL y dd=0 ZLOCAL r dd=0 ZLOCAL z dd=0 ZLOCAL sin_z dd=0 ZLOCAL cos_z x F0= F0= AND IF y F0= F0= AND IF +INF +INF ( both x,y = 0 ) EXIT ELSE y DD>F F0> IF dd=pi*2 ELSE dd=pi*2 ddnegate ENDIF EXIT ENDIF ELSE y F0= F0= AND IF x DD>F F0> IF dd=0 ELSE dd=pi ENDIF EXIT ENDIF ENDIF x y Z= IF y DD>F F0> IF dd=pi/4 EXIT ELSE dd=3*pi/4 EXIT ENDIF ENDIF x y ddnegate Z= IF y DD>F F0> IF dd=3*pi/4 EXIT ELSE dd=pi/4 ddnegate EXIT ENDIF ENDIF x dd^2 y dd^2 dd+ ddsqrt TO r x r dd/ ( xx ) y r dd/ ( yy ) \ Compute double precision approximation to atan. y DD>F x DD>F FATAN2 F>DD TO z ddover ddover ( F: xx yy xx yy -- ) dd- DD>F F0> IF \ Use Newton iteration 1. z' = z + (y - sin(z)) / cos(z) z ddsincos TO cos_z dd- cos_z dd/ ddswap dddrop z dd+ ELSE \ Use Newton iteration 2. z' = z - (x - cos(z)) / sin(z) dddrop z ddsincos ddswap TO sin_z dd- sin_z dd/ z ddswap dd- ENDIF ; : DDATAN ( F: aa a -- tt t ) dd=1 DDATAN2 ; : DDASIN ( aa a -- rr r ) dddup ddabs dddup dd1- DD>F F0> IF dddrop dddrop +INF +INF EXIT ENDIF dd1- F0= F0= AND IF DD>F F0> IF dd=pi/2 EXIT ELSE dd=pi/2 ddnegate EXIT ENDIF ENDIF dd=1 ddover dd^2 dd- ddsqrt ddatan2 ; : DDACOS ( aa a -- rr r ) dddup ddabs dddup dd1- DD>F F0> IF dddrop dddrop +INF +INF EXIT ENDIF dd1- F0= F0= AND IF DD>F F0> IF dd=0 EXIT ELSE dd=pi EXIT ENDIF ENDIF dd=1 ddover dd^2 dd- ddsqrt ddswap ddatan2 ; : 0.5*mul_pwr2 ( F: aa a -- cc c ) F2/ FSWAP F2/ FSWAP ; PRIVATE : DDSINH ( F: aa a -- rr r ) dddup F0= F0= AND ?EXIT dddup ddabs 0.05e F>DD dd- DD>F F0> IF ddexp dddup 1/dd dd- 0.5*mul_pwr2 EXIT ENDIF \ since a is small, using the above formula gives \ a lot of cancellation. So use Taylor series. ZLOCAL a a ZLOCAL t a dd^2 ZLOCAL r 1e FLOCAL m a ddeps ddf* DD>F FLOCAL thresh #32 >R a ( accu) BEGIN 2e +TO m t r dd* ( t) m F1- m F* ddf/ dddup TO t dd+ ( accu) t ddabs DD>F thresh F< -1 R+! R@ 0< OR UNTIL -R ; : DDCOSH ( F: aa a -- rr r ) dddup F0= F0= AND IF dddrop dd=1 EXIT ENDIF ddexp dddup 1/dd dd+ 0.5*mul_pwr2 ; : DDTANH ( F: aa a -- rr r ) dddup F0= F0= AND ?EXIT ddexp ZLOCAL ea ea 1/dd ZLOCAL inv_ea ea inv_ea dd- ea inv_ea dd+ dd/ ; : DDASINH ( F: aa a -- rr r ) dddup dd^2 dd1+ ddsqrt dd+ ddln ; : DDACOSH ( F: aa a -- rr r ) dddup dd1- DD>F F0< IF dddrop +INF +INF EXIT ENDIF dddup dd^2 dd1- ddsqrt dd+ ddln ; : DDATANH ( F: aa a -- rr r ) dddup ddabs dd1- DD>F F0> IF dddrop +INF +INF EXIT ENDIF dddup dd1+ ddswap dd=1 ddswap dd- dd/ ddln 0.5*mul_pwr2 ; -- ------------------------------------------------------------------------------------------------------------------------ : DD= ( F: dd1 dd2 -- ) ( -- bool ) dd- F0= F0= AND ; : DD<> ( F: dd1 dd2 -- ) ( -- bool ) dd- F0<> F0<> OR ; -- Untested. : DD< ( F: dd1 dd2 -- ) ( -- bool ) dd- FDROP F0< ; : DD> ( F: dd1 dd2 -- ) ( -- bool ) ddswap dd- FDROP F0< ; : DD>= ( F: dd1 dd2 -- ) ( -- bool ) dd- FDROP F0>= ; : DD<= ( F: dd1 dd2 -- ) ( -- bool ) ddswap dd- FDROP F0>= ; : DD0= ( F: dd -- ) ( -- bool ) dd=0 DD= ; : DD0< ( F: dd -- ) ( -- bool ) dd=0 DD< ; : DD0> ( F: dd -- ) ( -- bool ) dd=0 DD> ; : DD0>= ( F: dd -- ) ( -- bool ) dd=0 DD>= ; : DD0<= ( F: dd -- ) ( -- bool ) dd=0 DD<= ; : DDMIN ( F: dd1 dd2 -- ) ( -- bool ) DD2DUP DD< IF dddrop ELSE ddnip ENDIF ; : DDMAX ( F: dd1 dd2 -- ) ( -- bool ) DD2DUP DD< IF ddnip ELSE dddrop ENDIF ; -- ------------------------------------------------------------------------------------------------------------------------ #31 VALUE DDPRECISION : SET-DDPRECISION ( u -- ) #31 MIN 1 MAX TO ddprecision ; : getSign ( F: x xx -- x xx ) ( -- flag ) FOVER F0< ; PRIVATE : getPower ( F: x xx -- x xx ) ( -- n ) FOVER FABS FLOG F>S ; PRIVATE : normalize ( F: x xx -- |y + yy| ) ( n -- n' ) ddabs ( F: -- |x + xx| ) DUP dd=10 dd^n dd/ ( F: [|x+xx|]/10^n -- ) ( -- n ) BEGIN FOVER 10e F>= \ make sure it's between 1 and 10 WHILE dd=10 dd/ 1+ REPEAT BEGIN FOVER 1e F< \ make sure it's between 1 and 10 WHILE dd=10 dd* 1- REPEAT ; PRIVATE : peelDigit ( F: y yy -- y' yy' ) ( -- digit ) FOVER 1e-10 F+ \ fix inadvertently rounding down problem F>D 2DUP D>F F>DD dd- D>S ; PRIVATE : shiftBy10 10e ddf* ; PRIVATE \ shift dd# left one place -- Handle carry : cnvrt ( ... dn-1 dn -- dn-1' ) S>D # 2DROP ; PRIVATE : (ddout) ( F: |x+xx| -- y yy ) ( sign power -- caddr u ) S>D SWAP OVER DABS \ convert exponent <# FELEN @ 1- 0 MAX 0 ?DO # LOOP #S ROT IF '-' HOLD ELSE FESIGN @ IF '+' HOLD ENDIF ENDIF \ append digits of exponent FECHAR @ HOLD \ append 'e' ddprecision 0 ?DO peelDigit shiftBy10 LOOP peelDigit ddprecision 0 ?DO cnvrt LOOP FDEC @ HOLD cnvrt ROT IF '-' HOLD ELSE CASE FMSIGN @ 1 OF '+' HOLD ENDOF 2 OF BL HOLD ENDOF ENDCASE ENDIF #> ; PRIVATE : (DDFE.) ( F: x xx -- ) ( -- c-addr u ) \ create memory string for double-double in E-format FOVER F0= IF F2DROP S" 0e" EXIT ENDIF getSign ( -- s ) getPower ( F: x xx -- ) ( -- s n ) normalize ( F: y yy -- ) ( -- s n' ) (ddout) F2DROP ; -- Like FE. +E. : DDFE. ( F: x xx -- ) (DDFE.) TYPE SPACE ; \ display double-double in E-format : DDFE.R ( fieldwidth -- ) ( F: x xx -- ) >S (DDFE.) S> OVER - 0 MAX SPACES TYPE SPACE ; : DD+E.R ( mantissawidth -- ) ( F: rr r -- ) DDPRECISION >S >S SAVE-FFORMAT 2 FMSIGN ! 1 FESIGN ! '.' FDEC ! 'e' FECHAR ! 3 FELEN ! \ _/-x.xxxxxxxxxxxxxxE+/-yyy S SET-DDPRECISION S> 8 + DDFE.R RESTORE-FFORMAT S> SET-DDPRECISION ; : DD+E. ( F: rr r -- ) #29 DD+E.R ; : +DDCE.R ( n -- ) ( F: rr r -- ) DD>F E.R ; -- ------------------------------------------------------------------------------------------------------------------------ DOC (* Algorithm: Convert a string representing a dd# to 2 IEEE 64-bit floats internally 0. initialize buffers hi$ and lo$ 1. skip leading + sign, leave - flag on stack. 2. copy digits to hi$ buffer until dp found; let L be # digits to left of dp. 3. continue copying until #digits = MaxPlaces/2 or string is exhausted. 4. IF (string is exhausted) set lo$ to 0e0 ELSE append remaining digits to lo$ until #digits = MaxPlaces THEN 5. hi$ >FLOAT 1e0 exactmul lo$ >FLOAT 1e0 exactmul dd+ 6. get p = power of 10, let p' = p + L - MaxPlaces/2 7. p' dd=10 dd^n dd* *) ENDDOC #30 CONSTANT MaxPlaces CREATE lo# PRIVATE MaxPlaces CHARS ALLOT -- numerical string buffers CREATE hi# PRIVATE MaxPlaces CHARS ALLOT : isdigit ( c -- f ) '0' '9' 1+ WITHIN ; PRIVATE : dp? ( c -- f ) '.' = ; PRIVATE : dDeE? ( c -- f ) >UPC DUP 'D' = SWAP 'E' = OR ; PRIVATE : |power| ( c-addr u -- |p| ) \ abs value of pwr of 10 0. 2SWAP >NUMBER ( -- ud2 c-addr2 u2 ) 0= IF DROP D>S ELSE -2 ABORT" Bad exponent" ENDIF ; PRIVATE : do_exponent ( end beg -- p end' ) \ peel exponent from right LOCALS| beg end | 0 \ count on stack BEGIN end beg > ( -- n f1) end C@ isdigit ( -- n f1 f2) AND WHILE 1 -TO end \ dec adr 1+ \ inc count REPEAT end 1+ SWAP |power| ( |p|) end C@ '-' = IF NEGATE 1 -TO end ENDIF \ get sign of p end C@ '+' = IF 1 -TO end ENDIF \ bypass + sign end C@ dDeE? -1 AND +TO end \ advance past dDeE end C@ dDeE? -1 AND end + \ ( -- p end') DUP C@ DUP isdigit SWAP dp? OR \ digit|dp are legal 0= ABORT" Bad exponent" ; PRIVATE \ anything else Not Good : [dig|dp] ( char -- col# ) \ 0 = "other", 1 = digit, 2 = dp DUP isdigit 1 AND SWAP dp? 2 AND + ; PRIVATE : init_buffer ( c-addr -- ) \ initialize a string buffer DUP C0! 1+ MaxPlaces 1- '0' FILL ; PRIVATE : +CHAR ( c-addr c -- ) \ append char to counted $ OVER COUNT + C! 1 SWAP C+! ; PRIVATE 0 VALUE pre_dp PRIVATE \ # of digits to left of dp 0 VALUE post_dp PRIVATE \ # of digits to right of dp 0 VALUE (state) PRIVATE : +hi ( c -- ) hi# SWAP +CHAR \ append character 1 +TO pre_dp ; PRIVATE \ increment count : +hi|lo ( c -- ) pre_dp post_dp + MaxPlaces 2/ < IF hi# ELSE lo# ENDIF SWAP +CHAR \ append character 1 +TO post_dp ; PRIVATE \ increment post_dp : <> ( c col# -- ) (state) 3 * + ( -- c cell# ) CASE 0 OF ABORT" Non-digit in significand" ENDOF 1 OF +hi 0 TO (state) ENDOF 2 OF drop 1 TO (state) ENDOF 3 OF ABORT" Non-digit in significand" ENDOF 4 OF +hi|lo 1 TO (state) ENDOF 5 OF ABORT" One dp to a customer" ENDOF ENDCASE ; PRIVATE : leadingSign ( end beg -- sgn end beg' ) DUP C@ '-' OVER = ( -- end beg c f1 ) >R '+' = R@ OR ( -- end beg f2+f1 ) 1 AND + R> -ROT ; PRIVATE : initialize ( -- ) hi# init_buffer \ buffers lo# init_buffer lo# '.' +CHAR CLEAR pre_dp CLEAR post_dp \ counts CLEAR (state) ; PRIVATE -- digits to hi# and lo# : >(hi/lo) ( end beg -- ) BEGIN 2DUP >= ( -- end beg f ) pre_dp post_dp + MaxPlaces <= AND WHILE DUP C@ DUP [dig|dp] <> 1+ ( -- end beg+1 ) REPEAT 2DROP ; PRIVATE -- convert buffer to dd : buf->dd ( c_addr -- ) ( F: -- x xx ) COUNT >FLOAT 0= ABORT" Float conversion failed" 1e exactmul ; PRIVATE : adjustExponent ( p -- ) ( F: |y+yy| -- |x+xx| ) pre_dp + ( p+n1 ) pre_dp post_dp + ( p+n1 n1+n2) MaxPlaces 2/ MIN - ( p'=p+n1-MIN[n1+n2,MaxP/2] ) dd=10 dd^n dd* ; PRIVATE : >DD ( c-addr u -- ) ( F: -- x xx ) -TRAILING -LEADING OVER + 1- SWAP ( -- end beg ) initialize ( -- end beg ) leadingSign ( -- sgn end beg' ) TUCK do_exponent ROT ( -- sgn pwr end' beg' ) >(hi/lo) hi# buf->dd lo# buf->dd dd+ ( F: -- |y+yy| ) ( s p -- ) adjustExponent IF ddnegate ENDIF ; -- ------------------------------------------------------------------------------------------------------------------------ NESTING @ 1 = [IF] dd=1 3e ddf/ DDCONSTANT dd=1/3 dd=1 6e ddf/ DDCONSTANT dd=1/6 2e F>DD DDCONSTANT dd=2 : TEST ( -- ) PRECISION LOCAL old-prec #15 SET-PRECISION dd=1/3 dd=1/6 dd+ FSWAP CR CR ." Addition: 1/3 + 1/6" CR 6 spaces ." 5.00000000000000e-001 -1.54074395550979e-033 <- should be this" CR 5 spaces +E. SPACE +E. CR dd=1/3 dd=1/6 dd- FSWAP CR ." Subtraction: 1/3 - 1/6" CR 6 spaces ." 1.66666666666667e-001 9.25185853854297e-018 <- should be this" CR 5 spaces +E. SPACE +E. CR 6e F>DD dd=1/6 dd* FSWAP CR ." Multiplication: 6 * 1/6" CR 6 spaces ." 1.00000000000000e+000 0.00000000000000e+000 <- should be this" CR 5 spaces +E. SPACE +E. CR ." 6 * 1/3" 6e F>DD dd=1/3 dd* FSWAP CR 6 spaces ." 2.00000000000000e+000 0.00000000000000e+000 <- should be this" CR 5 spaces +E. SPACE +E. CR dd=1/3 dd=1/6 dd/ FSWAP CR ." Division: [1/3] / [1/6]" CR 6 spaces ." 2.00000000000000e+000 0.00000000000000e+000 <- should be this" CR 5 spaces +E. SPACE +E. CR CR ." [1/6] / [1/3]" dd=1/6 dd=1/3 dd/ FSWAP CR 6 spaces ." 5.00000000000000e-001 0.00000000000000e+000 <- should be this" CR 5 spaces +E. SPACE +E. CR ." Square root: sqrt[1/324]" dd=1/6 dd=1/3 dd* F2DUP dd* ddsqrt F2DUP FSWAP CR 6 spaces ." 5.55555555555556e-002 3.08395284618099e-018 <- should be this" CR 5 spaces +E. SPACE +E. ." <- sqrt[1/324] = 1/18" CR 6 spaces ." ... multiplied by 18" 18e F>DD dd* FSWAP CR 6 spaces ." 1.00000000000000e+000 0.00000000000000e+000 <- should be this" CR 5 spaces +E. SPACE +E. -9 10 SWAP ?DO dd=1 6e ddf/ I dd=10 dd^n dd* CR I . Tab EMIT Tab EMIT ddfe. LOOP CR dd=PI ddfe. old-prec SET-PRECISION ; : .TEST+ ( -- ) S" -11.1112222233333444445555566666e-45" >DD CR DDFE. ." ( -1.11112222233333444445555566666e-044 )" S" +11.1112222233333444445555566666e+45" >DD CR DDFE. ." ( 1.11112222233333444445555566666e+046 )" S" +11.1112222233333444445555566666dd" >DD CR DDFE. ." ( 1.11112222233333444445555566666e+001 )" S" +11.1112222233333444445555566666D" >DD CR DDFE. ." ( 1.11112222233333444445555566666e+001 )" S" 1111122.222333D" >DD CR DDFE. ." ( 1.11112222233300000000000000000e+006 )" ; : .TEST- ( -- ) S" +11.1112222233333444445555566666" >DD CR DDFE. ." ( 0.00000000000000000000000000000e-648 )" S" +11.1112222.233333444445555566666dd-45" >DD CR DDFE. ." ( Error: >dd One dp to a customer! )" S" +-11.1112222233333444445555566666dd-45" >DD CR DDFE. ." ( Error: >dd Non-digit in significand! )" ; : PRINT-TEST ( -- ) s" 1e-23" >DD ZLOCAL small s" 1e+23" >DD ZLOCAL big #34 #-34 DO dd=10 I dd^n I 0> IF big ELSE small ENDIF dd+ CR ddfe. LOOP CR CR #34 #-34 DO dd=10 I dd^n CR ddfe. LOOP ; : TEST-DD+ ( -- ) CR ." 10,000,000 executions of dd+ : " TIMER-RESET #10000000 0 DO dd=1 dd=1 dd+ DDDROP LOOP .ELAPSED ; : BENCH ( -- ) CR ." -- double-doubles (53) --" 53-bits! CR ." 10,000,000 executions of dd* : " TIMER-RESET #10000000 0 DO dd=2 dd=2 dd* DDDROP LOOP .ELAPSED CR ." 10,000,000 executions of dd/ : " TIMER-RESET #10000000 0 DO dd=10 dd=2 dd/ DDDROP LOOP .ELAPSED TEST-DD+ CR ." 10,000,000 executions of dd- : " TIMER-RESET #10000000 0 DO dd=1 dd=1 dd- DDDROP LOOP .ELAPSED CR ." 10,000,000 executions of ddsqrt : " TIMER-RESET #10000000 0 DO dd=1 ddsqrt DDDROP LOOP .ELAPSED CR ." 10,000,000 executions of dd^n : " TIMER-RESET #10000000 0 DO dd=1 3 dd^n DDDROP LOOP .ELAPSED CR ." 1,000,000 executions of ddexp : " TIMER-RESET #1000000 0 DO dd=2 ddexp DDDROP LOOP .ELAPSED CR ." 1,000,000 executions of ddln : " TIMER-RESET #1000000 0 DO dd=10 ddln DDDROP LOOP .ELAPSED CR ." 1,000,000 executions of ddsin : " TIMER-RESET #1000000 0 DO dd=2 ddsin DDDROP LOOP .ELAPSED CR ." 1,000,000 executions of ddcos : " TIMER-RESET #1000000 0 DO dd=10 ddcos DDDROP LOOP .ELAPSED CR ." 1,000,000 executions of ddtan : " TIMER-RESET #1000000 0 DO dd=2 ddtan DDDROP LOOP .ELAPSED CR ." 1,000,000 executions of ddasin : " TIMER-RESET #1000000 0 DO dd=1/3 ddasin DDDROP LOOP .ELAPSED CR ." 1,000,000 executions of ddacos : " TIMER-RESET #1000000 0 DO dd=1/3 ddacos DDDROP LOOP .ELAPSED CR ." 1,000,000 executions of ddatan : " TIMER-RESET #1000000 0 DO dd=1/3 ddatan DDDROP LOOP .ELAPSED CR ." 1,000,000 executions of ddasinh : " TIMER-RESET #1000000 0 DO dd=1/3 ddasinh DDDROP LOOP .ELAPSED CR ." 1,000,000 executions of ddacosh : " TIMER-RESET #1000000 0 DO dd=2 ddacosh DDDROP LOOP .ELAPSED CR ." 1,000,000 executions of ddatanh : " TIMER-RESET #1000000 0 DO dd=1/3 ddatanh DDDROP LOOP .ELAPSED CR ." -- doubles (53) --" 53-bits! CR ." 10,000,000 executions of f* : " TIMER-RESET #10000000 0 DO 1e 1e f* FDROP LOOP .ELAPSED CR ." 10,000,000 executions of f/ : " TIMER-RESET #10000000 0 DO 1e 1e f/ FDROP LOOP .ELAPSED CR ." 10,000,000 executions of f+ : " TIMER-RESET #10000000 0 DO 1e 1e f+ FDROP LOOP .ELAPSED CR ." 10,000,000 executions of f- : " TIMER-RESET #10000000 0 DO 1e 1e f- FDROP LOOP .ELAPSED CR ." 10,000,000 executions of fsqrt : " TIMER-RESET #10000000 0 DO 1e fsqrt FDROP LOOP .ELAPSED CR ." 10,000,000 executions of f^n : " TIMER-RESET #10000000 0 DO 1e 3 fpow FDROP LOOP .ELAPSED CR ." 1,000,000 executions of fexp : " TIMER-RESET #1000000 0 DO 2e fexp FDROP LOOP .ELAPSED CR ." 1,000,000 executions of fln : " TIMER-RESET #1000000 0 DO 10e fln FDROP LOOP .ELAPSED CR ." 1,000,000 executions of fsin : " TIMER-RESET #1000000 0 DO 2e fsin FDROP LOOP .ELAPSED CR ." 1,000,000 executions of fcos : " TIMER-RESET #1000000 0 DO 10e fcos FDROP LOOP .ELAPSED CR ." 1,000,000 executions of ftan : " TIMER-RESET #1000000 0 DO 2e ftan FDROP LOOP .ELAPSED CR ." 1,000,000 executions of fasin : " TIMER-RESET #1000000 0 DO 0.3e fasin FDROP LOOP .ELAPSED CR ." 1,000,000 executions of facos : " TIMER-RESET #1000000 0 DO 0.3e facos FDROP LOOP .ELAPSED CR ." 1,000,000 executions of fatan : " TIMER-RESET #1000000 0 DO 0.3e fatan FDROP LOOP .ELAPSED CR ." 1,000,000 executions of fasinh : " TIMER-RESET #1000000 0 DO 0.3e fasinh FDROP LOOP .ELAPSED CR ." 1,000,000 executions of facosh : " TIMER-RESET #1000000 0 DO 2e facosh FDROP LOOP .ELAPSED CR ." 1,000,000 executions of fatanh : " TIMER-RESET #1000000 0 DO 0.3e fatanh FDROP LOOP .ELAPSED CR ." -- doubles (80) --" 80-bits! CR ." 10,000,000 executions of f* : " TIMER-RESET #10000000 0 DO 1e 1e f* FDROP LOOP .ELAPSED CR ." 10,000,000 executions of f/ : " TIMER-RESET #10000000 0 DO 1e 1e f/ FDROP LOOP .ELAPSED CR ." 10,000,000 executions of f+ : " TIMER-RESET #10000000 0 DO 1e 1e f+ FDROP LOOP .ELAPSED CR ." 10,000,000 executions of f- : " TIMER-RESET #10000000 0 DO 1e 1e f- FDROP LOOP .ELAPSED CR ." 10,000,000 executions of fsqrt : " TIMER-RESET #10000000 0 DO 1e fsqrt FDROP LOOP .ELAPSED CR ." 10,000,000 executions of f^n : " TIMER-RESET #10000000 0 DO 1e 3 fpow FDROP LOOP .ELAPSED CR ." 1,000,000 executions of fexp : " TIMER-RESET #1000000 0 DO 2e fexp FDROP LOOP .ELAPSED CR ." 1,000,000 executions of fln : " TIMER-RESET #1000000 0 DO 10e fln FDROP LOOP .ELAPSED CR ." 1,000,000 executions of fsin : " TIMER-RESET #1000000 0 DO 2e fsin FDROP LOOP .ELAPSED CR ." 1,000,000 executions of fcos : " TIMER-RESET #1000000 0 DO 10e fcos FDROP LOOP .ELAPSED CR ." 1,000,000 executions of ftan : " TIMER-RESET #1000000 0 DO 2e ftan FDROP LOOP .ELAPSED CR ." 1,000,000 executions of fasin : " TIMER-RESET #1000000 0 DO 0.3e fasin FDROP LOOP .ELAPSED CR ." 1,000,000 executions of facos : " TIMER-RESET #1000000 0 DO 0.3e facos FDROP LOOP .ELAPSED CR ." 1,000,000 executions of fatan : " TIMER-RESET #1000000 0 DO 0.3e fatan FDROP LOOP .ELAPSED CR ." 1,000,000 executions of fasinh : " TIMER-RESET #1000000 0 DO 0.3e fasinh FDROP LOOP .ELAPSED CR ." 1,000,000 executions of facosh : " TIMER-RESET #1000000 0 DO 2e facosh FDROP LOOP .ELAPSED CR ." 1,000,000 executions of fatanh : " TIMER-RESET #1000000 0 DO 0.3e fatanh FDROP LOOP .ELAPSED 53-bits! ; DOC (* Win32Forth 4.2.0671 on 3 GHz PIV 10,000,000 executions of dd* : Elapsed time: 00:00:20.454 10,000,000 executions of dd/ : Elapsed time: 00:00:20.453 10,000,000 executions of dd+ : Elapsed time: 00:00:09.562 10,000,000 executions of dd- : Elapsed time: 00:00:11.453 10,000,000 executions of ddsqrt : Elapsed time: 00:00:21.094 ok MPE VFX Forth for Windows IA32 ( same machine) © MicroProcessor Engineering Ltd, 1998-2005 Version: 3.80 [build 1937] Build date: 29 September 2005 bench ( rebuild for 64bit float, cw@ $300 NOT AND $200 OR cw! ( 53-bit mode)) 10,000,000 executions of dd* : 1047 ms elapsed. 10,000,000 executions of dd/ : 1062 ms elapsed. 10,000,000 executions of dd+ : 391 ms elapsed. 10,000,000 executions of dd- : 391 ms elapsed. 10,000,000 executions of ddsqrt : 1625 ms elapsed. ok FORTH> ( iForth 2.0, same machine ) bench 10,000,000 executions of dd* : 0.512 seconds elapsed. 10,000,000 executions of dd/ : 0.512 seconds elapsed. 10,000,000 executions of dd+ : 0.244 seconds elapsed. 10,000,000 executions of dd- : 0.259 seconds elapsed. 10,000,000 executions of ddsqrt : 0.674 seconds elapsed. ok *) ENDDOC : TEST-ddexp ( -- ) -1e38 FLOCAL fmaxi 0e FLOCAL arg #1000000 0 DO #100 FCHOOSE FDUP F0= IF 1e-6 F+ ENDIF TO arg arg FEXP F>DD dddup arg F>DD ddexp dd- ddswap dd/ ddabs DD>F fmaxi FMAX TO fmaxi LOOP CR ." ddexp max delta = " fmaxi +E. ; : TEST-ddln ( -- ) -1e38 FLOCAL fmaxi 0e FLOCAL arg #1000000 0 DO #100 FCHOOSE FDUP F0= IF 1e-6 F+ ENDIF TO arg arg FLN F>DD dddup arg F>DD ddln dd- ddswap dd/ ddabs DD>F fmaxi FMAX TO fmaxi LOOP CR ." ddln max delta = " fmaxi +E. ; : TEST-ddlog ( -- ) -1e38 FLOCAL fmaxi 0e FLOCAL arg #1000000 0 DO #100 FCHOOSE FDUP F0= IF 1e-6 F+ ENDIF TO arg arg FLOG F>DD dddup arg F>DD ddlog dd- ddswap dd/ ddabs DD>F fmaxi FMAX TO fmaxi LOOP CR ." ddlog max delta = " fmaxi +E. ; : TEST-ddsin ( -- ) -1e38 FLOCAL fmaxi 0e FLOCAL arg #1000000 0 DO #100 FCHOOSE 50e F- FDUP F0= IF 1e-6 F+ ENDIF TO arg arg FSIN F>DD dddup arg F>DD ddsin dd- ddswap dd/ ddabs DD>F fmaxi FMAX TO fmaxi LOOP CR ." ddsin max delta = " fmaxi +E. ; : TEST-ddcos ( -- ) -1e38 FLOCAL fmaxi 0e FLOCAL arg #1000000 0 DO #100 FCHOOSE 50e F- FDUP F0= IF 1e-6 F+ ENDIF TO arg arg FCOS F>DD dddup arg F>DD ddcos dd- ddswap dd/ ddabs DD>F fmaxi FMAX TO fmaxi LOOP CR ." ddcos max delta = " fmaxi +E. ; : TEST-ddtan ( -- ) -1e38 FLOCAL fmaxi 0e FLOCAL arg #1000000 0 DO #100 FCHOOSE 50e F- FDUP F0= IF 1e-6 F+ ENDIF TO arg arg FTAN F>DD dddup arg F>DD ddtan dd- ddswap dd/ ddabs DD>F fmaxi FMAX TO fmaxi LOOP CR ." ddtan max delta = " fmaxi +E. ; : TEST-ddasin ( -- ) -1e38 FLOCAL fmaxi 0e FLOCAL arg #1000000 0 DO 2 FCHOOSE F1- FDUP F0= IF 1e-6 F+ ENDIF TO arg arg FASIN F>DD dddup arg F>DD ddasin dd- ddswap dd/ ddabs DD>F fmaxi FMAX TO fmaxi LOOP CR ." ddasin max delta = " fmaxi +E. ; : TEST-ddacos ( -- ) -1e38 FLOCAL fmaxi 0e FLOCAL arg #1000000 0 DO 2 FCHOOSE F1- FDUP F0= IF 1e-6 F+ ENDIF TO arg arg FACOS F>DD dddup arg F>DD ddacos dd- ddswap dd/ ddabs DD>F fmaxi FMAX TO fmaxi LOOP CR ." ddacos max delta = " fmaxi +E. ; : TEST-ddatan ( -- ) -1e38 FLOCAL fmaxi 0e FLOCAL arg #1000000 0 DO #100 FCHOOSE 50e F- FDUP F0= IF 1e-6 F+ ENDIF TO arg arg FATAN F>DD dddup arg F>DD ddatan dd- ddswap dd/ ddabs DD>F fmaxi FMAX TO fmaxi LOOP CR ." ddatan max delta = " fmaxi +E. ; : TEST-ddsinh ( -- ) -1e38 FLOCAL fmaxi 0e FLOCAL arg #1000000 0 DO #10 FCHOOSE 5e F- FDUP F0= IF 1e-6 F+ ENDIF TO arg arg FSINH F>DD dddup arg F>DD ddsinh dd- ddswap dd/ ddabs DD>F fmaxi FMAX TO fmaxi LOOP CR ." ddsinh max delta = " fmaxi +E. ; : TEST-ddcosh ( -- ) -1e38 FLOCAL fmaxi 0e FLOCAL arg #1000000 0 DO #10 FCHOOSE 5e F- FDUP F0= IF 1e-6 F+ ENDIF TO arg arg FCOSH F>DD dddup arg F>DD ddcosh dd- ddswap dd/ ddabs DD>F fmaxi FMAX TO fmaxi LOOP CR ." ddcosh max delta = " fmaxi +E. ; : TEST-ddtanh ( -- ) -1e38 FLOCAL fmaxi 0e FLOCAL arg #1000000 0 DO #100 FCHOOSE 50e F- FDUP F0= IF 1e-6 F+ ENDIF TO arg arg FTANH F>DD dddup arg F>DD ddtanh dd- ddswap dd/ ddabs DD>F fmaxi FMAX TO fmaxi LOOP CR ." ddtanh max delta = " fmaxi +E. ; : TEST-ddasinh ( -- ) -1e38 FLOCAL fmaxi 0e FLOCAL arg #1000000 0 DO #100 FCHOOSE 50e F- FDUP F0= IF 1e-6 F+ ENDIF TO arg arg FASINH F>DD dddup arg F>DD ddasinh dd- ddswap dd/ ddabs DD>F fmaxi FMAX TO fmaxi LOOP CR ." ddasinh max delta = " fmaxi +E. ; : TEST-ddacosh ( -- ) -1e38 FLOCAL fmaxi 0e FLOCAL arg #1000000 0 DO 100 FCHOOSE F1+ TO arg arg FACOSH F>DD dddup arg F>DD ddacosh dd- ddswap dd/ ddabs DD>F fmaxi FMAX TO fmaxi LOOP CR ." ddacosh max delta = " fmaxi +E. ; : TEST-ddatanh ( -- ) -1e38 FLOCAL fmaxi 0e FLOCAL arg #1000000 0 DO 100e 101 FCHOOSE F- FDUP F0= IF 1e-6 F+ ENDIF TO arg arg FATANH F>DD dddup arg F>DD ddatanh dd- ddswap dd/ ddabs DD>F fmaxi FMAX TO fmaxi LOOP CR ." ddatanh max delta = " fmaxi +E. ; : TESTS ( -- ) TEST-ddexp TEST-ddln TEST-ddlog TEST-ddsin TEST-ddcos TEST-ddtan TEST-ddasin TEST-ddacos TEST-ddatan TEST-ddsinh TEST-ddcosh TEST-ddtanh TEST-ddasinh TEST-ddacosh TEST-ddatanh ; : exactmulbench ( -- ) TIMER-RESET #10000000 0 ?DO PI PI exactmul dddrop LOOP .ELAPSED ; : PTEST ( F: eps -- ) dd=10 -2 dd^n ZLOCAL tiny 0 LOCAL errors FLOCAL eps dd=0 #1000 0 DO tiny dd+ dddup dddup 1/dd (DDFE.) >dd dd* dd1- DD>F FABS eps F> 1 AND +TO errors LOOP dddrop CR ." Found " errors DEC. ." errors." ; :ABOUT CR ." Try: TEST .TEST+ .TEST- PRINT-TEST TESTS BENCH exactmulbench or ( F: eps -- ) PTEST" ; [ELSE] :ABOUT CR ." *** DD-FP words ***" CR ." dd@ dd! dd0! dd=1 dd=10 ddconstant ddvalue ddvariable" CR ." ddswap ddover dddup dddrop " CR ." dd+ dd- dd+! dd-! ddf+ " CR." dd/ dd* ddf* ddf/ " CR ." ddabs dd^2 ddsqrt dd^n ddnegate " CR ." dd+e. dd+e.r (ddfe.) ddfe. ddfe.r >dd " CR ." dd= dd< dd<= dd> dd>= " CR ." dd0= dd0< dd0>= dd0> dd0>= " CR ." ddexp ddln ddlog ddsincos ddatan2 " CR ." ddsin ddcos ddtan ddasin ddacos ddatan " CR ." ddsinh ddcosh ddtanh ddasinh ddacosh ddatanh " ; [THEN] NESTING @ 1 = [IF] .ABOUT -dd-fp CR [THEN] DEPRIVE (* End of Source *)