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