(     Title:  ULP'S [Units in the Last Place]
       File:  fulps.fs
    Version:  0.5.0
    Revised:  July 28, 2009
    Authors:  David N. Williams

The date above may reflect unlogged, cosmetic changes.

Version 0.5.0
21Jul09 * Started.
22Jul09 * Experimental FNEXTUP variations.
25Jul09 * Replaced FNEXTUP with Andrew Haley's version.
        * Replaced hexfloat.fs with rawfloat.fs.
26Jul09 * Removed test for zero in F-/ULP.
28Jul09 * Removed redundancy in FNEXTUP +inf case.
        * Made +INF a conditional definition.

The words FULP and F-/ULP have been tested only by inspection.
)

s" rawfloat.fs" included
decimal

bits/float 64 <>
[IF] cr .( ***This code requires 64 bits/float.) ABORT [THEN]

[UNDEFINED] F2DUP  [IF] : F2DUP  fover fover ;  [THEN]
[UNDEFINED] F2DROP [IF] : F2DROP fdrop fdrop ;  [THEN]
[UNDEFINED] F-ROT  [IF] : F-ROT  frot frot ;    [THEN]
[UNDEFINED] +INF   [IF] 1e 0e f/ fconstant +inf [THEN]

: NAN?  ( f: r -- s: nan? )  fdup f= 0= ;

: FNEXTUP  ( f: x -- y)  \ aph
   fdup nan? fdup +inf f= or IF EXIT THEN
   0.0e f+      \ convert -0 to +0
   f>raw dup 0< IF  1 0 d-  ELSE  1 0 d+  THEN raw>f ;

: FNEXTDOWN  ( f: x -- y )  fnegate fnextup fnegate ;

: FULP  ( f: x -- ulp[x] )
(
Leave the size of one ulp for x.
)
  fdup fnextup fover f- fswap  ( f: [x-up - x] x)
  fdup fnextdown f-            ( f: [x-up - x] [x - x-down])
  fmin ;

: F-/ULP  ( f: x ref -- [x-ref]/ulp[ref] )
(
Leave the difference between x and ref in units of one ref-ulp.
If the reference value ref is accurate to say 0.5 ulp, and x is
the value being tested, then the output is the accuracy of x in
ref-ulps, up to about +/- 0.5 ulp.
)
  fdup fulp ( f: x ref ulp[ref]) f-rot f- fswap f/ ;