\ FLOATING and FLOATING-EXT words not implemented here are:
\ FLITERAL REPRESENT F** FACOS FACOSH FALOG FASIN FASINH FATAN FATAN2
\ FATANH FCOS FCOSH FEXP FEXPM1 FLN FLNP1 FLOG FSIN FSINCOS FSINH FTAN
\ and FTANH
\ 0: Initial release 11/07/2002, Brad Eckert
\ 1: Removed locals from F*, PAD from F/. Added FSQRT.
\ 2: Fixed FROUND.
\ 3: Fixed FPICK, converted standard words to upper case, sped up stack.
\ 4: Changed F2* and F2/ to use fshift. You can use fshift for quick *2^n.
\ Replaced DO..LOOP with 0 ?DO LOOP structure.
\ 5: Changed 0 ?DO LOOP to use EVALUATE (more optimizer-friendly).
\ 6: Fixed a bug in F/.
\ 7: Removed extra FDUPs from F<. Fixed a bug in F+.
\ 8: Changed -ROT to ROT ROT for strict ANS compliance.
\ 9: Made F# state smart. Fixed bug in FSQRT.
\ 10: Made 4tH version, Hans Bezemer
\ 11: Fixed FLOOR, added F=.
\ 12: Added IEEE 754 exceptions.
\ 13: Fixed FSQRT (didn't clear stack on error).
\ 14: Fixed F+ (did FABS when 2OS is zero).
\ 15: Added F>S.
\ 16: Fixed F~ (didn't compare zero and negative zero properly).
\ 17: Fixed F0< (didn't compare zero and negative zero properly).
\ 18: Added FSIGN? and mmaxdigits.
\ 19: Made &esign more portable.
\ 20: Replaced 3 PICK with >R >R OVER R> R> ROT.
\ 21: Fixed FLOOR again.
\ 22: Fixed t2/, cleared sign bit. Optimized frshift.
\ 23: Made constants of &esign and maxdigits. Optimized FSQRT.
\ 24: Increased FLOATING-STACK
\ 25: Split up library. Made fsp public.
[UNDEFINED] floats [IF]
[UNDEFINED] 2@ [IF] include lib/anscore.4th [THEN]
[UNDEFINED] um* [IF] include lib/mixed.4th [THEN]
[UNDEFINED] FLOATING-STACK [IF]
16 CONSTANT FLOATING-STACK \ size of floating point stack
[THEN]
\ Floating point representation uses three cells:
\ 0 upper mantissa
\ 1 lower mantissa
\ 2 exponent: MSB = sign of mantissa, others = signed exponent
3 CELLS CONSTANT c/float
c/float [NEGATE] CONSTANT -c/float
c/float CONSTANT FLOAT
VARIABLE sigdigits \ # of digits to display (see F., R.)
VARIABLE fsp \ stack pointer
FLOAT FLOATING-STACK [*] ARRAY fstak \ floating point stack
FLOAT 2 [*] ARRAY ftemp \ temporary float variables
/cell 8 [*] CONSTANT bits/cell \ 16-bit or 32-bit Forth
bits/cell -1 [+] 602 [*] 1000 [/] CONSTANT maxdigits
\ max. digits
max-n CONSTANT &unsign
(error) CONSTANT &sign
&sign 2 [/] &sign [+] CONSTANT &esign \ equals &sign 1 rshift
&esign [SIGN] 1 [=] [NOT] [IF] [ABORT] [THEN]
\ abort when not positive
\ A variable called ferror stores error codes in the event of an error.
\ Execution is not interrupted, but a bad result could be generated. You
\ should check for a non-zero error code at appropriate times.
VARIABLE ferror \ last seen error
\ Standard IEEE 754 exceptions
0 enum FE.NOERRORS \ No FP errors
enum FE.DIVBYZERO \ Division by zero
enum FE.OVERFLOW \ Overflow
enum FE.UNDERFLOW \ Underflow
enum FE.INVALID \ Invalid operation
constant FE.INEXACT \ Inexact result
: ud2/ ( d -- d' )
D2/ &unsign AND ;
: frshift ( count 'man -- ) \ right shift mantissa
SWAP 0 MAX bits/cell 2 [*] MIN
>R DUP 2@ R> 0 ?DO ud2/ LOOP ROT 2! ;
: exp>sign ( exp -- exp' sign )
DUP &unsign AND \ mask off exponent
DUP &esign AND 2* OR \ sign extend exponent
SWAP &sign AND ; \ get sign of mantissa
: sign>exp ( exp sign -- exp' )
SWAP &unsign AND OR ;
: t2* ( triple -- triple' )
D2* ROT DUP 0< 1 AND >R 2* ROT ROT R> 0 D+ ;
: t2/ ( triple -- triple' )
OVER 1 AND 0 SWAP IF INVERT THEN &sign AND >R D2/ ROT
1 RSHIFT &unsign AND R> OR ROT ROT ;
\ MHL|C
: t+ ( t1 t2 -- t3 )
2>R >R ROT 0 R> 0 D+ 0 2R> D+
ROT >R D+ R> ROT ROT ;
: tneg ( d flag -- t ) \ conditionally negate d to 3-cell t
0<> >R 2DUP OR 0<> R> AND IF DNEGATE -1 ELSE 0 THEN ;
: lstemp ( -- ) \ 6-cell left shift FTEMP
ftemp 6 CELLS + 0 ( *LSB carry . . )
6 0 ?DO >R -1 CELLS + \ LSB in high memory
DUP @ 0 D2* SWAP R> + ( a co n ) \ ROL
ROT TUCK ! SWAP \ a carry
LOOP 2DROP ;
: *norm ( triple exp -- double exp' ) \ normalize triple
>R BEGIN DUP 0< 0=
WHILE t2* R> 1- >R
REPEAT &sign 0 0 t+ \ round
ROT DROP R> ;
: longdiv ( DA DB -- DA/DB ) \ fractional divide
0 0 ftemp 2!
bits/cell 2* 1+ \ long division
0 ?DO 2OVER 2OVER DU< \ double/double
IF 0
ELSE 2DUP 2>R D- 2R> 1 \ a b --> a-b b
THEN 0 ftemp 2@ D2* D+ ftemp 2!
2SWAP D2* 2SWAP
LOOP DU< 0= 1 AND 0 \ round
ftemp 2@ D+ ;
\ The rest can be left in high level Forth...
: 'm1 ( -- addr ) fsp @ 3 cells - ; \ -> 1st stack mantissa
: 'm2 ( -- addr ) fsp @ 6 cells - ; \ -> 2nd stack mantissa
: 'm3 ( -- addr ) fsp @ 9 cells - ; \ -> 3rd stack mantissa
: 'e1 ( -- addr ) fsp @ 1 cells - ; \ -> 1st stack exponent
: 'e2 ( -- addr ) fsp @ 4 cells - ; \ -> 2nd stack exponent
: fmove ( a1 a2 -- ) c/float 0 DO OVER I + @ OVER I + ! LOOP 2DROP ;
: m1get ( addr -- ) 'm1 fmove ; \ move to M1
: m1sto ( addr -- ) 'm1 SWAP fmove ; \ move from M1
: expx1 ( -- exp sign ) 'e1 @ exp>sign ;
: expx2 ( -- exp sign ) 'e2 @ exp>sign ;
: ftemp' ( -- addr ) ftemp 2 CELLS + ;
: >exp1 ( exp sign -- ) sign>exp 'e1 ! ;
: fshift ( n -- ) expx1 >R + R> >exp1 ; \ F = F * 2^n
\ A normalized float has an unsigned mantissa with its MSB = 1
: normalize ( -- )
'm1 2@ 2DUP OR 0=
IF 2DROP \ Zero is a special case.
ELSE 0 ROT ROT expx1 >R *norm
R> >exp1 'm1 2!
THEN ;
\ *S Floating Point Words
\ *+
: F2* ( fs: r1 -- r3 ) 1 fshift ;
: F2/ ( fs: r1 -- r3 ) -1 fshift ;
: PRECISION ( -- n ) sigdigits @ ; \ max decimal digits/double
: SET-PRECISION ( n -- ) maxdigits MIN sigdigits ! ;
: FCLEAR ( -- ) fstak fsp ! FE.NOERRORS ferror ! ;
: FDEPTH ( -- ) fsp @ fstak - c/float / ;
: FDUP ( fs: r -- r r ) c/float fsp +! 'm2 m1get ;
: FDROP ( fs: r -- ) -c/float fsp +! ;
: FNEGATE ( fs: r1 -- r3 ) 'e1 @ &sign XOR 'e1 ! ;
: D>F ( d -- | -- r ) FDUP DUP 0< IF DNEGATE &sign ELSE 0 THEN
'e1 ! 'm1 2! normalize ;
: F>D ( -- d | r -- ) expx1 >R DUP 1- 0<
IF NEGATE &unsign AND 'm1 frshift 'm1 2@ R> IF DNEGATE THEN
ELSE R> 2DROP -1 &unsign FE.OVERFLOW ferror !
THEN FDROP ; \ overflow: return maxdint
: S>F ( n -- | -- r ) DUP ABS U>D ROT 0< IF DNEGATE THEN D>F ;
: F>S ( r -- | -- n ) F>D 2DUP D0< >R DABS D>U R> IF NEGATE THEN ;
: FLOAT+ ( n1 -- n2 ) c/float + ;
: FLOATS ( n1 -- n2 ) c/float * ;
: F@ ( a -- | -- r ) FDUP m1get ;
: F! ( a -- | r -- ) m1sto FDROP ;
: FSIGN? ( -- t | r -- r) 'e1 @ 0< ;
: F0= ( -- t | r -- ) 'm1 2@ OR 0= FDROP ;
: F0< ( -- t | r -- ) 'm1 2@ OR 0<> 'e1 @ 0< AND FDROP ;
: FABS ( fs: r1 -- r3 ) 'e1 @ 0< IF FNEGATE THEN ;
: FALIGN ( -- ) ALIGN ;
: FALIGNED ( addr -- addr ) ALIGNED ;
: FSWAP ( fs: r1 r2 -- r2 r1 )
'm1 ftemp fmove 'm2 m1get ftemp 'm2 fmove ;
: FOVER ( fs: r1 r2 -- r1 r2 r3 )
c/float fsp +! 'm3 m1get ;
: FPICK ( n -- ) ( fs: -- r )
c/float fsp +! 2 + -c/float * fsp @ + m1get ;
: FNIP ( fs: r1 r2 -- r2 ) FSWAP FDROP ;
: FROT ( fs: r1 r2 r3 -- r2 r3 r1 )
'm3 ftemp fmove
'm2 'm3 c/float 2* 0 DO OVER I + @ OVER I + ! LOOP 2DROP
ftemp m1get ;
: F+ ( fs: r1 r2 -- r3 ) \ Add two floats
FDUP F0= IF FDROP EXIT THEN \ r2 = 0
FOVER F0= IF FNIP EXIT THEN \ r1 = 0
expx1 >R expx2 >R - DUP 0<
IF NEGATE 'm1 frshift 0 \ align exponents
ELSE DUP 'm2 frshift
THEN 'e2 @ +
'm2 2@ R> tneg
'm1 2@ R> tneg
t+ DUP >R ( exp L M H | sign )
DUP IF t2/ IF DNEGATE THEN 'm2 2! 1+
ELSE DROP 'm2 2!
THEN R> &sign AND sign>exp 'e2 !
FDROP normalize ;
: F- ( fs: r1 r2 -- r3 ) \ Subtract two floats
FNEGATE F+ ;
: F< ( -- flag ) ( F: r1 r2 -- ) \ Compare two floats
F- F0< ;
: F= ( -- flag ) ( F: r1 r2 -- ) \ Compare two floats
F- F0= ;
: FLOOR FDUP F0< FDUP F>D D>F FOVER F- F0= 0= AND F>D ROT
IF 1 U>D DNEGATE D+ THEN D>F ;
: FROUND ( fs: r1 -- r2 ) \ Round to the nearest integer
expx1 >R NEGATE 1- 'm1 frshift \ convert to half steps
'm1 2@ 1 U>D D+ SWAP -2 AND SWAP \ round
'm1 2! -1 R> >exp1 normalize ; \ re-float
: FMIN ( F: r1 r2 -- rmin ) FOVER FOVER F<
IF FDROP ELSE FNIP THEN ;
: FMAX ( F: r1 r2 -- rmax ) FOVER FOVER F<
IF FNIP ELSE FDROP THEN ;
: F* ( fs: r1 r2 -- r3 ) \ Multiply two floats
'm1 2@ 'm2 2@
OVER >R ftemp' 2!
OVER >R ftemp 2!
R> R> OR
\ need long multiply?
IF FTEMP CELL+ @ FTEMP' CELL+ @ UM* &sign 0 D+ NIP \ round
FTEMP @ FTEMP' @ UM*
FTEMP CELL+ @ FTEMP' @ UM* 0 t+
FTEMP @ FTEMP' CELL+ @ UM* 0 t+
ELSE 0 ftemp @ ftemp' @ UM* \ lower parts are 0
THEN 2DUP OR >R >R OVER R> R> ROT OR \ zero?
IF expx1 >R expx2 >R + bits/cell 2* + *norm
R> R> XOR sign>exp 'e2 !
ELSE DROP \ zero result
THEN 'm2 2! FDROP ;
: F/ ( fs: r1 r2 -- r3 ) \ Divide r1 by r2
FDUP F0= \ div by 0, man= umaxint
IF FDROP -1 -1 'm1 2! FE.DIVBYZERO ferror !
'e1 @ &sign AND &esign 1- OR 'e1 ! \ exponent = maxint/2
ELSE 'm2 2@ 'm1 2@
DU< 0= IF 1 'm2 frshift F2/ THEN \ divisor<=dividend
'm1 CELL+ @
IF 'm2 2@ ud2/ 'm1 2@ ud2/ longdiv \ max divisor = dmaxint
ELSE 0 'm2 2@ 'm1 @ DUP >R UM/MOD \ 0 rem quot | divisor
ROT ROT R@ UM/MOD ROT ROT R> 1 RSHIFT &unsign AND U>
IF 1 U>D D+ THEN \ round
THEN expx2 >R expx1 >R - bits/cell 2* -
>R 'm2 2! R> R> R> XOR sign>exp 'E2 !
FDROP
THEN ;
: F~ ( f: r1 r2 r3 -- ) ( -- flag ) \ f-proximate
FDUP F0< \ Win32forth version
IF FABS FOVER FABS 3 FPICK FABS F+ F* \ r1 r2 r3*(r1+r2)
FROT FROT F- FABS FSWAP F<
ELSE FDUP F0=
IF FDROP 'e1 @ 0< 'e2 @ 0< = F- F0= AND
ELSE FROT FROT F- FABS FSWAP F<
THEN
THEN ;
: FSQRT ( fs: r -- r' ) \ Square root
expx1 IF drop FE.INVALID ferror ! EXIT THEN \ sqrt of negative
DUP 1 AND DUP >R + \ round up exponent/2
2/ bits/cell - 0 >exp1
ftemp 6 BOUNDS DO 0 I ! LOOP
'm1 2@ \ x
R> IF ud2/ THEN \ exponent was rounded
ftemp 2 CELLS + 2! \ x*2^(2*bits/cell)
0 0 bits/cell 2 [*] \ p = 0
0 ?DO lstemp lstemp \ shift left x into a 2 places
D2* \ shift left p one place
2DUP D2* ftemp 2@ D<
IF ftemp 2@ 2OVER D2* 1 U>D D+ D-
ftemp 2! \ a:=a-(2*p+1)
1 U>D D+ \ p:=p+1
THEN
LOOP 'm1 2! normalize ;
\ For fixed-width fields, it makes sense to use fsplit and (F.) for output.
\ fsplit converts a float to integers suitable for pictured numeric output.
\ Fracdigits is the number of digits to the right of the decimal.
\ Left here as a public word because it requires carnal knowledge of the
\ ANS float library.
: fsplit ( F: r -- ) ( fracdigits -- sign Dint Dfrac )
\ Split float into integer component parts.
>R expx1 NIP FABS \ int part must fit in a double
FDUP F>D 2DUP D>F F- \ get int, leave frac
2 0 R> 0 ?DO D2* 2DUP D2* D2* D+ LOOP \ 2 * 10^precision
D>F F* F>D 1 U>D D+ ud2/ ; \ round
[DEFINED] 4TH# [IF]
hide c/float
hide -c/float
hide sigdigits
hide fstak
hide ftemp
hide bits/cell
hide &unsign
hide &sign
hide &esign
hide ud2/
hide frshift
hide exp>sign
hide sign>exp
hide t2*
hide t2/
hide t+
hide tneg
hide lstemp
hide *norm
hide longdiv
hide 'm1
hide 'm2
hide 'm3
hide 'e1
hide 'e2
hide fmove
hide m1get
hide m1sto
hide expx1
hide expx2
hide ftemp'
hide >exp1
hide fshift
hide normalize
[THEN]
[THEN]
|
