I found a case where I really need: f :: Float -> Float where f x is the least y such that x < y
This seems to be the problem of finding the unnormalized epsilon: the smallest positive number one can meaningfully add to the given number x. If that x is 1.0, we're talking about the epsilon. There are several approaches: - Consult the file /usr/include/float.h and find FLT_EPSILON and DBL_EPSILON there. One can write a Haskell program that parses that file, if one so wishes. - apply the formula epsilon = base^(1-p) - write a program everybody has written in good old Fortran days
nextafter1:: Float -> Float nextafter1 0 = 0 nextafter1 x | GHC.Float.isNaN x = x nextafter1 x | GHC.Float.isInfinite x = x nextafter1 x = try (abs x) where try d = let d1 = d/2 in if x + d1 == x then improve d1 d else try d1 -- improve a b -- we know that the number we seek is between a and b -- try to see if there is a number in between that -- still can be meaningfully added to x improve a b = let middle = (a+b)/2 in if middle == b || middle == a then x + b else if x + middle > x then improve a middle else improve middle b
Standard caveat applies: some overly smart compilers arrange for the whole iteration to be done with CPU registers, which often have extra precision. Therefore, one might need to disable optimizations or play other tricks to force storing/reloading the number to/from memory.
test_delta:: Float -> Float -> Bool test_delta x y = y > x && (let x'=x+(0.9*(y-x)) in x'==y || x'==x)
*Main> test_delta (1.0::Float) (nextafter1 (1.0::Float)) True *Main> test_delta (0.1::Float) (nextafter1 (0.1::Float)) True *Main> test_delta (-42.3e20::Float) (nextafter1 (-42.3e20::Float)) True *Main> (nextafter1 (-42.3e50::Float)) -Infinity We can also deal with denormalized numbers: *Main> GHC.Float.isDenormalized (12.0e-40::Float) True *Main> test_delta (12.0e-40::Float) (nextafter1 (12.0e-40::Float)) True
On Tue, 21 Oct 2003, Lennart Augustsson wrote:
So this has been a while, but i think that decodeFloat, incrementing the mantissa, encodeFloat might work.
But what if the mantissa is 0xffffff? We are at mercy of encodeFloat doing the right thing. Furthermore, this technique doesn't seem to work with denormalized numbers *Main> GHC.Float.decodeFloat (1.0::Float) (8388608,-23) *Main> GHC.Float.decodeFloat (nextafter1 (1.0::Float)) (8388609,-23) *Main> GHC.Float.decodeFloat (0.1::Float) (13421773,-27) *Main> GHC.Float.decodeFloat (nextafter1 (0.1::Float)) (13421774,-27) So far, so good. However, *Main> GHC.Float.decodeFloat (12.0e-40::Float) (13701584,-153) *Main> GHC.Float.decodeFloat (nextafter1 (12.0e-40::Float)) (13701600,-153) Indeed, *Main> GHC.Float.encodeFloat 13701584 (-153) == (12.0e-40::Float) True *Main> GHC.Float.encodeFloat 13701585 (-153) == (12.0e-40::Float) True so incrementing the mantissa by one doesn't actually work in all circumstances.