Next floating point number

[Bengt Richter]

On Sat, 17 Dec 2005 14:23:38 +1100, Steven D'Aprano <steve at REMOVETHIScyber.com.au> wrote:

>I'm looking for some way to get the next floating point number after any
>particular float.
>
>(Any mathematicians out there: I am aware that there is no "next real
>number". But floats are not real numbers, they only have a finite
>precision.)
>
>According to the IEEE standard, there should be a routine like next(x,y)
>which returns the next float starting from x in the direction of y.
>
>Unless I have missed something, Python doesn't appear to give an interface
>to the C library's next float function. I assume that most C floating
>point libraries will have such a function.
>
>So I came up with a pure Python implementation, and hoped that somebody
>who had experience with numerical programming in Python would comment.
>
>
>
>def nextfloat(x, y=None):
>    """Return the next floating point number starting at x in the
>    direction of y. If y is not given, the direction is away from zero.
>    """
>    from operator import add, sub
>    from math import frexp # Returns a tuple (mantissa, exponent)
>    x = float(x)
>    if y is None:
>        if x < 0.0:
>            op = sub
>        else:
>            op = add
>    elif y > x:
>        op = add
>    elif y < x:
>        op = sub
>    else:
>        raise ValueError("direction float is equal to starting float")
>    # We find the next float by finding the smallest float epsilon which 
>    # is distinguishable from x. We know epsilon will be a power of two, 
>    # because floats are implemented internally in binary.
>    # We get our initial estimate for epsilon from the floating point 
>    # exponent of x.
>    epsilon = 2.0**(frexp(x)[1]) # epsilon > abs(x)
>    lasteps = epsilon # Save a copy of our last epsilon for later use.
>    while op(x, epsilon) != x:
>        lasteps = epsilon
>        epsilon /= 2.0
>    # When we get here, we've halved epsilon one time too many.
>    # We can't just double it, because that fails for the case where x is
>    # zero - epsilon will underflow to zero. So we need to save and use
>    # the previous iteration of epsilon.
>    return op(x, lasteps)
>
I wonder if this won't work (for IEEE 754 double that is)

from math import frexp
def nextf(x, y):
    f,e = frexp(x)
    if (f==0.5 or f==-0.5) and x>=y: eps = 2.0**-54
    else: eps = 2.0**-53
    if x<y: return (f+eps)*2.0**e
    else: return (f-eps)*2.0**e


Regards,
Bengt Richter