Double-double floating-point package for Forth

\ Modified for KForth

[UNDEFINED] DXFORTH [IF]

: SYSTEM ;
: APPLICATION ;

\ : \\  SOURCE-ID IF  BEGIN REFILL 0= UNTIL  THEN ;

: END  POSTPONE EXIT  POSTPONE THEN ; IMMEDIATE

: /CHAR ( a u -- a2 u2 c )  1 /STRING OVER 1- C@ ;

: >DIGIT ( char base -- u true | char 0 )
  OVER  [CHAR] 0  -  DUP 10 17 WITHIN 0= IF  DUP 16 U>
  IF  7 -  THEN  DUP ROT U< IF  NIP TRUE  END  DUP  THEN
  2DROP FALSE ;

: >CHAR ( u -- char )  DUP 9 > IF  7 +  THEN  [CHAR] 0  + ;

: (NUMBER)  0 0 2SWAP >NUMBER ;

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

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

: TOKEN  BL WORD COUNT ;

[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-05 changes:
\
\ - renamed package to DDFLOAT so you know who to blame
\ - fix bugs in DDFS.
\ - replace input routine >DD
\ - 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 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 ----
SYSTEM

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

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

: DDLITERAL
    FSWAP  POSTPONE FLITERAL  POSTPONE FLITERAL
; IMMEDIATE

APPLICATION
\ ---------------------------------------------------

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
134217729 S>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
;

: 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
    /CHAR  10 >DIGIT WHILE
    S>F 0.0E  tmp DD@  dd10. DD*  DD+  tmp DD!
    dpf @  exp +!
  REPEAT  DROP  -1 /STRING  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  END
  tmp DD@  R> IF  DDNEGATE  THEN  TRUE ;

SYSTEM

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

APPLICATION

[DEFINED] DXFORTH [IF] BEHEAD exp getexp [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
    getsign IF  DDNEGATE  THEN   ( |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>
;

CREATE digits  32 CELLS ALLOT

: copydigits ( -- )  \ to HOLD buffer
  1 CELLS digits +  31 CELLS digits +  DO
    I @ 0< IF  10 I +!  -1 I 1 CELLS - +!  THEN  \ fix negative
    I @ >CHAR HOLD
  -1 CELLS +LOOP  [char] .  HOLD  digits @ >CHAR HOLD
;

: peeldigits ( y yy -- )  \ to array
  digits  32 0 DO  >R
    FOVER F>S  R> 2DUP !  1 CELLS + >R
    S>F 0e0  DD-  dd10. DD*  R>
  LOOP  DROP  DDDROP
;

: (DDFS.) ( dd -- c-addr u )  \ string double-double in E-format
  FOVER F0= IF
    0 0 >R
  ELSE
    getsign >R
    getpower        ( x xx n)
    normalize       ( y yy n')
  THEN
  DUP ABS S>D                \ convert exponent
  <#  #S  2DROP  SIGN        \ append  digits of exponent
  [char] E HOLD
  peeldigits  copydigits
  R> IF  [char] - HOLD  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 dd1. dd10.
  BEHEAD getsign peeldigits
[THEN]


1 [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.1112222233333444445555566666" test
\ 0.0000000000000000000000000000000 dd 0  ok
\ <- Must have exponent field  (N/A in this version)

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# 2 DD# 3 DD/ CR DDFS.
DD# 9.9999999999999999999999987654321e CR DDFS.
DD# -9.9999999999999900000000123456789e CR DDFS.
DD# -9.999999999999e CR DDFS.

[THEN]