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