Double Extended-precision floating point arithmetic - 2026-09-27

\ Last revised: 2026-09-27

DECIMAL

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

[UNDEFINED] DXFORTH [IF]

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

[ELSE]  AKA SHOLD HOLDS  [THEN]


\ XDFLOAT.F

\ Double Extended-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
\
\ Changes:
\
\ - updated to use Intel 80-bit extended precision
\ - renamed package to DDFLOAT so you know who to blame
\ - fix bugs in DDFS.
\ - replace input routine >DD  Note: now returns a flag!!!
\ - Add DD# for ease of number input
\ - rename r128! r128@ real*16 to DD! DD@ DDVARIABLE
\ - add DD< DD> DD= DDLITERAL (DDFS.)
\ - 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 80-bit internal operations
\               ^^^^

CR .( Double Extended-precision FP ) CR

1 FLOATS 10 - [IF]
  CR .( ERROR: Requires Extended-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]

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

\ 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

4294967297.0000000000E 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< ( F: x xx y yy -- ) ( -- f )  ddcomp -1 = ;
: DD> ( F: x xx y yy -- ) ( -- f )  ddcomp 1 = ;
: DD= ( F: x xx y yy -- ) ( -- f )  ddcomp 0= ;

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

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

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

: DD^N ( F: 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 var

: /digs ( a u -- a' u' )
  BEGIN  DUP  WHILE
    OVER C@  DIGIT? WHILE
    S>F 0.0E  var DD@  dd10 DD*  DD+  var 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 @ +
  >R  var DD@  dd10 R> DD^N  DD*  var DD!
;

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

[DEFINED] DXFORTH [IF] SYSTEM [THEN]

\ Interpret/compile dd number
: DD# ( "number" )
  BL WORD COUNT  >DD 0= ABORT" DD-BAD"
  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.


38 CONSTANT maxdig

CREATE peeled  maxdig 1+ CELLS ALLOT
CREATE dig$    maxdig 1+ ALLOT

VARIABLE sgn  \ sign

: getsign ( F: x xx -- x xx ) ( -- )
    FOVER  F0<  sgn ! ;

: getpower ( F: x xx -- x xx ) ( -- )
    FOVER   FABS  FLOG
    ( FDUP  F0< IF  1e0 F-  THEN )
    F>S  exp ! ;

: normalize ( F: x xx -- |y yy| ) ( -- f )
    DDABS

    \ an error here results in NaN NaN
    dd10 exp @ DD^N  DD/    ( [|x|+|xx|]/10^n )

    FOVER FDUP F= IF

    \ make sure it's between 1 and 10
      BEGIN   FOVER   10e0  F>
      WHILE   dd10  dd/  1 exp +!  REPEAT

      BEGIN   FOVER   1e0   F<
      WHILE   dd10  dd*  -1 exp +!  REPEAT

      TRUE  EXIT

    THEN  0
;

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

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

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

: (DDFS.) ( F: x xx -- ) ( -- a u )  \ string double-double in E-format
  <#  FOVER F0= IF  DDDROP  S" 0.E0" HOLDS  ELSE
    getsign
    getpower
    normalize IF
      peeldigits  copydigits
      dig$ maxdig ( a u)  0digit-fix ( a' u')
      exp @  DUP ABS 0 #S 2DROP  SIGN  [char] E HOLD  \ exp
      1-  OVER 1+ SWAP  HOLDS  [char] . HOLD  C@ HOLD  \ mant
      sgn @ IF  [char] - HOLD  THEN
    ELSE  DDDROP  S" **DD**" HOLDS  THEN
  THEN  0 0 #>
;

: DDFS. ( F: x xx -- )  \ display double-double in E-format
  (DDFS.) TYPE SPACE ;


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

\ pi = 3.1415926535897932384626433832795028841971693993751...

DD# 3.141592653589793238462643383279502884197e DDCONSTANT DDPI

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


[DEFINED] DXFORTH [IF]
  BEHEAD split exactmul
  BEHEAD dtemp ddtemp1
  BEHEAD a3 ddcomp
  BEHEAD exp /exp
  BEHEAD dd1 dd10
  BEHEAD maxdig 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

s" +11.1112222233333444445555566666d" test

s" +11.1112222233333444445555566666D" test

s" 1111122.222333D" test

CR

DD# 0 CR DDFS.
DD# 1.2345678901234567890123456789012345678e CR DDFS.
DD# 2 DD# 3 DD/ CR DDFS.
DD# 9.9999999999999999999999999999987654321e CR DDFS.
DD# -9.9999999999999900000000000000123456789e CR DDFS.
DD# -9.999999999999e CR DDFS.
DD# 99999999999999999999999999999999999999.e CR DDFS.

[THEN]