# 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

: 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  THEN )
F>S
;

: 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: y yy y)
1e-10 F+                        ( f: y yy y+epsilon)
\ fixes rounding problem
F>D   2DUP  D>F  0e0  dd-
D>S
;

\ ******** bug in peelDigit fixed Wednesday, March 08 2006 by JVN
\          thanks to M. Hendrix for finding it

: 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      ( include)
\ needs fsm2.f         ( include)

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.

\        end C@  dDeE?  -1  AND  ( DUP  >R)
\                end +  TO  end  \ advance past dDeE
\        end C@  dDeE?  -1  AND  ( DUP)  end +           ( -- p end')
\            ( SWAP  R>  OR  0=  ABORT" Bad exponent" )  \ no dDeE !

;

: [dig|dp]  ( char -- col#)     \ 0 = "other", 1 = digit, 2 = dp
DUP  digit?  1 AND  SWAP  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 \$
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!"  ;

FALSE [IF]
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
[THEN]

\ FALSE [IF]
0 VALUE  (state)        : >state  TO  (state)  ;

: <<hi/lo>>     ( c col# --)
(state) 3 *  +          ( c cell#)
case
0   OF  err1   2 >state     ENDOF
1   OF  +hi    0 >state     ENDOF
2   OF  drop   1 >state     ENDOF
3   OF  err1   2 >state     ENDOF
4   OF  +hi|lo 1 >state     ENDOF
5   OF  err2   3 >state     ENDOF
endcase
;
\ [THEN]

: 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)
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]