\ 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

: 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


\ 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)

: 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  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 integer power
    \ return 1 if n=0, dd^{-|n|} if n<0
       dd=1   ddswap        ( f: 1e0 0e0 x xx )
       DUP  0=   IF  dddrop  drop  EXIT  THEN
       DUP  0<   SWAP  ABS   ( -- sign |n| )
         BEGIN   DUP  0>  WHILE
                 DUP  1 AND   IF ddtuck  dd*  ddswap THEN dd^2
         REPEAT  dddrop  DROP
       IF  dd=1  ddswap  dd/  THEN