128-bit Forth floating point arithmetic package 2026-09-09

\ Tested on Gforth, Win32Forth (8-byte floats).

\ Misc utilities

: DIGIT? ( char -- u true | 0 )
  [CHAR] 0 -  DUP 10 17 WITHIN 0= IF
    DUP 16 U> IF  7 -  THEN
    DUP BASE @ U< IF  TRUE  EXIT  THEN
  THEN  DROP  FALSE ;

[DEFINED] DXFORTH [IF]  AKA SHOLD HOLDS  [ELSE]

: /SIGN ( a u -- a' u' f )
  DUP IF  OVER C@  DUP [CHAR] + =  SWAP [CHAR] - =
  DUP >R  OR  NEGATE /STRING  R>  EXIT  THEN  0 ;

: /NUMBER ( a u -- a' u' d|ud )
  /SIGN >R 0 0 2SWAP >NUMBER 2SWAP R> IF DNEGATE THEN ;

[THEN]


\ DDFLOAT.F

\ Doubled-precision floating point arithmetic
\
\ Reference
\
\  JVN's 128-bit Forth floating point arithmetic package:
\  http://galileo.phys.virginia.edu/classes/551.jvn.fall01/ddarith.htm
\  http://galileo.phys.virginia.edu/classes/551.jvn.fall01/dd_io.htm
\
\ 2026-09-09 changes:
\
\ - renamed package to DDFLOAT so you know who to blame
\ - fix bugs in DDFS.
\ - replace input routine >DD  Note: now returns a flag
\ - rename r128! r128@ real*16 to DD! DD@ DDVARIABLE
\ - add DD< DD> DD= DDLITERAL (DDFS.) DD#
\ - remove separate float stack requirement
\
\ Material not (c) Julian V. Noble is hereby placed
\ in the PUBLIC DOMAIN. Use at your own risk.

\ ---------------------------------------------------
\     (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
\               ^^^^

CR .( Double-Double Precision FP) CR

1 FLOATS 8 - [IF]
  CR .( ERROR: Requires Double-Precision FP system) ABORT
[THEN]

MARKER -double

[UNDEFINED] F>S   [IF] : F>S  F>D D>S ; [THEN]
[UNDEFINED] S>F   [IF] : S>F  S>D D>F ; [THEN]
[UNDEFINED] FTUCK [IF] : FTUCK  FSWAP FOVER ; [THEN]
[UNDEFINED] F-ROT [IF] : F-ROT  FROT FROT ; [THEN]


\ ---------------------------------------- LOAD, STORE
: DD@ ( adr -- dd )  DUP >R  F@  R> FLOAT+  F@ ;
: DD! ( dd adr -- )  DUP >R  FLOAT+ F!  R> F! ;
\ ------------------------------------ END LOAD, STORE

\ ----------------------------------- data types ----

[DEFINED] DXFORTH [IF] SYSTEM [THEN]

: DDVARIABLE   \ create a double-double variable
    FALIGN HERE  2 FLOATS ALLOT  CONSTANT ;

: DDCONSTANT   \ create a double-double constant
    CREATE  FALIGN HERE  2 FLOATS ALLOT  DD!
    DOES>  FALIGNED DD@ ;

: DDLITERAL
    FSWAP  POSTPONE FLITERAL  POSTPONE FLITERAL
; IMMEDIATE

[DEFINED] DXFORTH [IF] APPLICATION [THEN]

\ ---------------------------------------------------

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

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. D>F FCONSTANT split

fvariable t     fvariable a1    fvariable a2    fvariable b1
fvariable b2    fvariable b21   fvariable b22   fvariable ta
fvariable tb

fvariable q     fvariable qq    fvariable x     fvariable xx
fvariable y     fvariable yy    fvariable z1    fvariable z2

fvariable aaa   fvariable 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
;

: DDABS    ( f: x xx -- |x+xx|)
    FOVER  F0<   IF  ddnegate  THEN
;

: 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+
;

\ ----------------------------------- stack ops -----
FVARIABLE dtemp

DDVARIABLE  ddtemp   DDVARIABLE  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  DD!  dddup  ddtemp1 DD!
           ddtemp  DD@         ddtemp1 DD@  ;

: DDTUCK   ddtemp  DD!  ddtemp1 DD!
           ddtemp  DD@  ddtemp1 DD@
           ddtemp  DD@  ;

fvariable a3    fvariable a4    fvariable b3    fvariable b4

: ddcomp ( a3 a4 b3 b4 -- -1|0|1 )
  b4 F!  b3 F!  a4 F!  a3 F!
  a3 F@ b3 F@ F< IF  -1 EXIT  THEN
  a3 F@ b3 F@ F= IF
    a4 F@ b4 F@ F< IF  -1 EXIT  THEN
    a4 F@ b4 F@ F= IF  0 EXIT  THEN
  THEN  1 ;

: DD< ( dd1 dd2 -- f )  ddcomp -1 = ;
: DD> ( dd1 dd2 -- f )  ddcomp 1 = ;
: DD= ( dd1 dd2 -- f )  ddcomp 0= ;

\ ---------------------------------------------------

1e0  0e0 DDCONSTANT dd1.
10e0 0e0 DDCONSTANT dd10.

: DD^2      ( x xx -- [x+xx]^2 )
     dddup  dd*  ;

: DD^N      ( x xx n -- [x+xx]^n )
    \ raise dd to integer power
    \ return 1 if n=0, dd^{-|n|} if n<0
       >R
       dd1.   DDSWAP        ( f1e0 0e0 x xx )
       R> DUP 0= IF  DROP  DDDROP  EXIT  THEN
       DUP 0< >R ( sign)  ABS ( |n| )
         BEGIN   DUP  0>  WHILE
                 DUP >R  1 AND IF DDTUCK  DD*  DDSWAP THEN dd^2
                 R> 2/
         REPEAT  DROP  DDDROP
       R> IF  dd1.  DDSWAP  DD/  THEN
;

\ Input routine adapted from:
\ http://dxforth.webhop.org/finput.html

VARIABLE exp  \ exponent
VARIABLE dpf  \ decimal point

DDVARIABLE tmp

: /digs ( a u -- a' u' )
  BEGIN  DUP  WHILE
    OVER C@  DIGIT? WHILE
    S>F 0.0E  tmp DD@  dd10. DD*  DD+  tmp DD!
    dpf @  exp +!  1 /STRING
  REPEAT  THEN ;

: /mant ( a u -- a' u' flag )
  TUCK  /digs  DUP IF
    OVER C@ [CHAR] . = IF
      -1 dpf !  1 /STRING  /digs
    THEN
  THEN  ROT OVER - dpf @ + ;

: /exp ( a u -- a' u' )
  DUP IF
    OVER C@  33 OR  [CHAR] e = ( 'D' 'E' 'd' 'e')
    1 AND /STRING
  THEN
  /NUMBER DROP exp @ +
  BEGIN  ?DUP WHILE  DUP 0<
    IF    1+  tmp DD@  dd10.  DD/
    ELSE  1-  tmp DD@  dd10.  DD*  THEN  tmp DD!
  REPEAT ;

\ Like >FLOAT but returns a dd number
: >DD ( c-addr u -- dd true | false )
  0.0E 0.0E tmp DD!  0 exp !  0 dpf !
  -TRAILING  /SIGN >R  DUP IF
    /mant IF  /exp  THEN
  ELSE  0=  THEN
  NIP IF  R> DROP  0  EXIT  THEN
  tmp DD@  R> IF  DDNEGATE  THEN  TRUE ;

[DEFINED] DXFORTH [IF] SYSTEM [THEN]

\ Interpret/compile dd number
: DD# ( "number" )
  BL WORD COUNT  >DD 0= ABORT" bad float"
  STATE @ IF  POSTPONE DDLITERAL  THEN ; IMMEDIATE

[DEFINED] DXFORTH [IF] APPLICATION [THEN]

\ This is an ANS Forth program requiring the
\   FLOAT, FLOAT EXT, FILE and TOOLS EXT wordsets.
\
\ Environmental dependences:
\       Assumes independent floating point stack

\ Based on I/O subroutines from
\
\   DDFUN: A Double-Double Floating Point Computation Package
\   IEEE Fortran-90 version
\   Version Date:  2005-01-26
\
\   Author:
\
\      David H. Bailey
\      NERSC, Lawrence Berkeley Lab
\      Mail Stop 50B-2239
\      Berkeley, CA 94720

\ Output a dd# .

\ Algorithm:
\
\  1. Determine sign and power of 10
\  2. Normalize so that significand lies between 1 and <10
\  3. Peel digits off from left by converting to integer
\       multiplying by dd10 and repeating
\  4. Construct output string digit by digit, append
\       sign and exponent, and display.

: getsign   ( x xx -- x xx f )
    FOVER  F0<  ;

: getpower   ( x xx -- x xx n)
    FOVER  FABS  FLOG
    FLOOR  F>S  ;

: normalize  ( x xx n -- y yy n')
    >R
    DDABS                       ( |x| |xx| )
    dd10. R@ DD^N  DD/          ( [|x|+|xx|]/10^n )

  \ make sure it's between 1 and <10
    BEGIN   DDDUP dd10. DD< 0=
    WHILE   dd10. DD/  R> 1+ >R  REPEAT

    BEGIN   DDDUP dd1. DD<
    WHILE   dd10. DD*  R> 1- >R  REPEAT
    R>
;

32 CONSTANT maxdig

CREATE peeled  maxdig CELLS ALLOT

: peeldigits ( y yy -- )  \ to integer array
  maxdig 0 DO
    FOVER  F>S  DUP  peeled  maxdig 1- I - CELLS + !
    S>F 0e0  DD-  dd10. DD*
  LOOP  DDDROP
;

CREATE dig$  maxdig ALLOT

: copydigits ( -- )  \ to character array
  peeled  maxdig 0 DO
    DUP CELL+  SWAP @  DUP 0< IF ( -ve fix-up)
      10 +  OVER  -1 SWAP +!
    THEN  [CHAR] 0 +  dig$ maxdig 1- I - + C!
  LOOP DROP
;

: 0digit-fix ( a u -- a' u' n )  \ leading digit fix-up
  OVER C@  CASE
    [CHAR] 0    OF  1 /STRING  -1  ENDOF
    [CHAR] 9 1+ OF  OVER  [CHAR] 1  SWAP C!  1  ENDOF
    0 SWAP
  ENDCASE
;

: (DDFS.) ( dd -- c-addr u )  \ string double-double in E-format
  <#  FOVER F0= IF  DDDROP  S" 0.E0"  HOLDS  ELSE
    getsign >R
    getpower        ( x xx n)
    normalize >R    ( y yy)
    peeldigits  copydigits
    dig$ maxdig ( a u)  0digit-fix  R> +  \ adjust exp
    ( p) DUP ABS 0 #S 2DROP  SIGN  \ exponent
    [char] E HOLD  ( a u)
    1-  OVER 1+ SWAP  HOLDS  [char] .  HOLD  C@ HOLD
    R> ( s) IF  [char] - HOLD  THEN
  THEN  0 0 #>
;

: DDFS. ( dd -- )  \ display double-double in E-format
  (DDFS.) TYPE SPACE ;

[DEFINED] DXFORTH [IF]
  BEHEAD split exactmul
  BEHEAD dtemp ddtemp1
  BEHEAD a3 ddcomp
  BEHEAD exp /exp
  BEHEAD dd1. dd10.
  BEHEAD getsign 0digit-fix
[THEN]


0 [IF]  CR .( Testing ) CR

: test ( a u -- )  2DUP >DD
  CR IF DDFS. 2DROP ELSE TYPE ."  failed" THEN ;

s" -11.1112222233333444445555566666d-45" test

s" +-11.1112222233333444445555566666d-45" test
\ Error: >dd Non-digit in significand!

s" +11.1112222.233333444445555566666d-45" test
\ Error: >dd One dp to a customer!

s" +11.1112222233333444445555566666d+45" test
\ +1.1111222223333344444555556666604 dd 46  ok

s" +11.1112222233333444445555566666d" test
\ +1.1111222223333344444555556666603 dd 1  ok

s" +11.1112222233333444445555566666D" test
\ +1.1111222223333344444555556666603 dd 1  ok

s" 1111122.222333D" test
\ +1.1111222223330000000000000000000 dd 6  ok

CR

DD# 0 CR DDFS.
DD# 1.2345678901234567890123456789012e CR DDFS.
DD# 2 DD# 3 DD/ CR DDFS.
DD# 9.9999999999999999999999987654321e CR DDFS.
DD# -9.9999999999999900000000123456789e CR DDFS.
DD# -9.999999999999e CR DDFS.
DD# 99999999999999999999999999999999.e CR DDFS.

[THEN]