\ fel 05.3.10 NAB
\ float exponents & logarithms

needs core-ext

module fel

\ words for float access
fvariable Ftmp

public:

: 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/ ;

private:

: (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 ;

public:

: 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 ;

end-module