ddarith.f
\ 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
\ the fpu must be set to 64-bit internal operations
\ ^^^^
MARKER -double
\ ---------------------------------------- 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
\ = 3.141592653589793238462643383279
\ pi = 3.1415926535897932384626433832795028841971693993751...
\ 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 ;
FALSE [IF] \ 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= NOT [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= NOT [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
[THEN]
\ Exact multiplication
134217729 S>F FCONSTANT split
9 fvariables: t a1 a2 b1 b2 b21 b22 ta tb
: ftuck FSWAP FOVER ;
: f-rot FROT FROT ;
8 fvariables: q qq x xx y yy z1 z2
2 fvariables: aaa bbb
: 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+
;
\ 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])
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