```dd_io.f \ Doubled-precision arithmetic  I/O words \ \ --------------------------------------------------- \     (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 \ 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 \      Email: dhbailey@lbl.gov MARKER  -ddio \ 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. 10e0 0e0  ddconstant dd10 : ddabs    ( f: x xx -- |x+xx|)     FOVER  F0<   IF  ddnegate  THEN  ; : getSign   ( f: x xx -- x xx)  ( -- f)     FOVER  F0<  ; : getPower   ( f: x xx -- x xx) ( --  n)     FOVER   FABS  FLOG     FDUP  F0<     IF  1e0 F- F>S  ELSE  F>S  THEN ; : normalize  ( f: x xx -- y yy) ( n -- n')     ddabs           ( f: |x| |xx| )     DUP  dd10  dd^n  dd/        ( f: [|x|+|xx|]/10^n ) ( n)   \ make sure it's between 1 and 10     BEGIN   FOVER   10e0  F>     WHILE   dd10  dd/   1 +  REPEAT     BEGIN   FOVER   1e0   F<     WHILE   dd10  dd*   1 -  REPEAT ;   ( f: [y+yy] = [x+xx]/10^n' )    ( n') : peelDigit ( f: y yy -- y' yy')    ( digit)     FOVER  F>D   2DUP  D>F  0e0  dd-     D>S ; : shiftBy10       dd10  dd*  ; : <ddout>   ( f: |x+xx| -- )   ( sgn power)     S>D  SWAP  OVER  DABS       \ convert exponent     <#  #S  ROT SIGN            \ append  digits of exponent        BL       HOLD        [char] d HOLD        [char] d HOLD        BL       HOLD            \ append dd        31 0 DO   peelDigit  shiftBy10  LOOP  peelDigit        31 0 DO   S>D   #  2DROP   LOOP        [char] .  HOLD  S>D  #  2DROP        ROT  IF  [char] -  ELSE  [char] +  THEN   HOLD     #>  TYPE SPACE ; : ddfs.     ( f: x xx -- )  \ display double-double in E-format     FOVER   F0=  IF  ." 0.0 DD 0"  EXIT  THEN   \ handle 0     getSign                     ( s)     getPower        ( f: x xx ) ( s n)     normalize       ( f: y yy)  ( s n')     <ddout>                     \ display     FDROP  FDROP                \ clean up fstack ; 15 SET-PRECISION : test     SWAP  DO  1e0 0e0 3e0 0e0  dd/  I dd10 dd^n  dd*               CR I .  TAB TAB  ddfs. LOOP ; : test1     SWAP  DO  1e0 0e0 6e0 0e0  dd/  I dd10 dd^n  dd*               CR I .  TAB TAB  ddfs. LOOP ; \ Convert a string representing a dd# to internal form \   as 2 IEEE 64-bit floats \ Algorithm: \ \   0. initialize buffers hi\$ and lo\$ \   1. skip leading + sign, leave - flag on stack. \   2. copy digits to hi\$ buffer until dp found; \        let L be # digits to left of dp. \   3. continue copying until #digits = MaxPlaces/2 \        or string is exhausted. \   4. IF   (string is exhausted) set lo\$ to 0e0 \      ELSE append remaining digits to lo\$ \           until #digits = MaxPlaces  THEN \   5. hi\$ >FLOAT   1e0  exactmul \      lo\$ >FLOAT   1e0  exactmul \      dd+ \   6. get p = power of 10, let p' = p + L - MaxPlaces/2 \   7. p' dd10 dd^n  dd* \ : \$ends   ( c-adr u -- end beg)   \ convert c-adr u to ends \    OVER  +   1-  SWAP ;          ( end beg) \ : ends->count   ( end beg -- c-adr u)  TUCK  -  1+  ; needs vector1.f needs fsm2.f 30 CONSTANT  MaxPlaces CREATE lo#  MaxPlaces  CHARS ALLOT  \ numerical string buffers CREATE hi#  MaxPlaces  CHARS ALLOT : digit?    ( c -- f)   \ test for digit     [ CHAR 0 ] LITERAL  [ CHAR 9  1+ ]  LITERAL  WITHIN ; : dp?   [CHAR] .  =    ; : dDeE?     ( c -- f)     [CHAR] d  OVER  =               OVER  [CHAR] D  =  OR               OVER  [CHAR] e  =  OR               SWAP  [CHAR] E  =  OR ; : |power|     ( c-addr count -- |p|) \ abs value of pwr of 10     0 0  2SWAP  >NUMBER   ( -- ud2 c-addr2 u2)     0=  IF  DROP  D>S  ELSE  ." Bad exponent"  ABORT   THEN ; : do_exponent  ( end beg -- p end')    \ peel exponent from right     LOCALS| beg end |         0                               \ count on stack             BEGIN   end beg >           ( n f1)                     end C@  digit?      ( n f1 f2)                     AND             WHILE   end 1-  TO end      \ dec adr                     1+                  \ inc count             REPEAT         end 1+  SWAP   |power|          ( |p|)         end C@  [char] -  =                 \ get sign of p         IF  NEGATE  end  1-  TO end   THEN  ( -- p)         end C@  [char] +  =             \ bypass + sign         IF  end  1-  TO  end  THEN         end C@  dDeE?  -1  AND  end +  TO  end  \ advance past dDeE         end C@  dDeE?  -1  AND  end +           ( -- p end')         DUP  C@  DUP  digit?  SWAP  dp?  OR  \ digit|dp are legal             0=  ABORT" Bad exponent"         \ anything else n.g. ; : [dig|dp]  ( char -- col#)     \ 0 = "other", 1 = digit, 2 = dp     digit?  1 AND  OVER  dp?  2 AND  +  ; : init_buffer   ( c-addr -- )   \ initialize a string buffer     0 OVER  C!  1+  MaxPlaces  1-  [CHAR] 0  FILL ; : +char     ( c-addr c -- )     \ append char to counted \$     OVER  COUNT  +              ( c-addr c c-addr+u)     C!   DUP C@   1+  SWAP  C!  \ store c, inc \$len ; 0 VALUE  pre_dp         \ # of digits to left of dp 0 VALUE  post_dp        \ # of digits to right of dp : +hi   ( c -- )     hi#  SWAP  +char        \ append character     pre_dp  1+  TO  pre_dp  \ increment count ; : +hi|lo   ( c -- )     pre_dp  post_dp +  MaxPlaces 2/  <     IF    hi#   ELSE   lo#   THEN     SWAP  +char                 \ append character     post_dp  1+  TO  post_dp    \ increment post_dp ; : err1   TRUE  ABORT" Non-digit in significand!"  ; : err2   TRUE  ABORT" One dp to a customer!"  ; 3 wide  fsm:  <<hi/lo>>     ( c col# --) \ input       other   |    digit      |      dp    | \ state  -------------------------------------------   ( 0)   ||  err1 >2  || +hi     >0   ||  DROP >1   ( 1)   ||  err1 >2  || +hi|lo  >1   ||  err2 >3   ( 2)   ( abnormal termination w/ error1 )   ( 3)   ( abnormal termination w/ error2 ) ;fsm : leadingSign   ( end beg -- sgn end beg')     DUP C@  [CHAR] -  OVER =    ( end beg c f1)     >R  [CHAR] +  =  R@  OR     ( end beg f2+f1)     1 AND +  R>  -ROT ; : initialize     hi#  init_buffer            \ buffers     lo#  init_buffer     lo#  [CHAR] .  +char     0 TO pre_dp   0 TO post_dp  \ counts     0 >state  <<hi/lo>>         \ fsm ; : >(hi/lo)  ( end beg -- )  \ digits to hi# and lo#     BEGIN   2DUP >=         ( end beg f)             pre_dp  post_dp  +  MaxPlaces  <=  AND     WHILE   DUP  C@  DUP  [dig|dp]  <<hi/lo>>             1+              ( end beg+1)     REPEAT  2DROP ; : buf->dd   ( c_addr )  ( f: x xx)  \ convert buffer to dd     COUNT  >FLOAT  0=  ABORT" Float conversion failed"     1e0  exactmul ; : adjustExponent    ( pwr)  ( f: |y+yy| -- |x+xx|)     pre_dp  +               ( p+n1)     pre_dp post_dp +        ( p+n1 n1+n2)     MaxPlaces 2/  MIN  -    ( p'=p+n1-MIN[n1+n2,MaxP/2] )     dd10  dd^n  dd* ; : >dd   ( c-addr u -- )  ( f: -- x xx)     OVER  +   1-  SWAP          ( end beg)  \ \$ends     initialize                  ( end beg)     leadingSign                 ( sgn end beg')     TUCK  do_exponent   ROT     ( sgn pwr end' beg')     >(hi/lo)                    \ digits -> hi/lo buffers     hi#  buf->dd    lo#  buf->dd    dd+     ( f: |y+yy|)     ( s p)  adjustExponent      ( s)     IF  ddnegate  THEN          \ adjust sign ; FALSE [IF] Examples: s" -11.1112222233333444445555566666dd-45" >dd cr ddfs. -1.1111222223333344444555556666603 dd -44  ok s" +-11.1112222233333444445555566666dd-45" >dd cr ddfs. Error: >dd Non-digit in significand! s" +11.1112222.233333444445555566666dd-45" >dd cr ddfs. Error: >dd One dp to a customer! s" +11.1112222233333444445555566666dd+45" >dd cr ddfs. +1.1111222223333344444555556666604 dd 46  ok s" +11.1112222233333444445555566666" >dd cr ddfs. +0.0000000000000000000000000000000 dd 0  ok  <- Must have exponent field s" +11.1112222233333444445555566666dd" >dd cr ddfs. +1.1111222223333344444555556666603 dd 1  ok s" +11.1112222233333444445555566666D" >dd cr ddfs. +1.1111222223333344444555556666603 dd 1  ok s" 1111122.222333D" >dd cr ddfs. +1.1111222223330000000000000000000 dd 6  ok [THEN] ```