Next floating point number

[Steven D'Aprano]

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)



Thanks,



-- 
Steven.