\ FEL float exponents & logarithms
needs core-ext
fvariable Ftmp
\ words for float access (thanks, Neal)
: F> ( F: r -- -- r. ) Ftmp dup f! 2@ ;
: >F ( r. -- -- F: r ) Ftmp 2! Ftmp F@ ;
: @F ( F: r -- F: r -- r. ) fdup f> ;
: F1/X ( F: r -- F: r ) 1E0 FSWAP F/ ;
: (FEXP) FDUP
-0.1413161E-3 F*
1.329882E-3 F+ FOVER F*
-0.830136E-2 F+ FOVER F*
0.41657347E-1 F+ FOVER F*
-0.1666653E0 F+ FOVER F*
0.5E0 F+ FOVER F*
-1E0 F+ F*
1E0 F+ F1/X ;
: FEXP ( F: r -- F: r )
FDUP F0< ABS FABS
FDUP 40E0 FSWAP F< ABORT" EXP overflow"
FDUP 1.44269504E0 F*
FDUP F>D DROP >R
FLOOR 0.6931472E0 F* F- (FEXP)
R> 65 + 32768 >F F*
IF F1/X THEN ;
: FLN ( F: r -- F: r ) fdup f0= fdup f0<
or abort" logarithm error"
fdup 1e0 f< if f1/x -1 else 0 then
@f drop 127 and 65 - >r
f> swap -128 and 65 or swap >f
1e0 f- fdup
-0.64535443E-02 f*
0.36088496E-01 f+ fover f*
-0.95329395E-01 f+ fover f*
0.16765408E00 f+ fover f*
-0.24073381E00 f+ fover f*
0.33179904E00 f+ fover f*
-0.49987414E00 f+ fover f*
0.99999648E00 f+ f* r> 0 d>f
0.6931472E00 f* f+
if fnegate then
;
: FX^Y ( F: x y -- F: x^y )
fswap fln f* fexp ;
\ this one is faster, more accurate
\ for integer exponents
: FX^N ( F: r -- n -- F: r^n -- )
dup 0= if drop fdrop 1e0
else dup abs fdup
1 ?do fover f* loop fswap fdrop
0< if f1/x then
then ;