Enum on Float/Double

[oleg at pobox.com]

> 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.