How do you get the next value in the floating-point sequence?

UPDATE:

Turns out this is a duplicate question (which comes up in google as result #2 for the search "c++ nextafter python"): Increment a python floating point value by the smallest possible amount

The accepted answer provides some solid solutions.

ORIGINAL ANSWER:

Certainly this isn't the perfect solution but using cython just a few lines will allow you to wrap the existing C++ function and use it in Python. I've compiled the below code and it works on my ubuntu 11.10 box.

First, a .pyx file (I called mine nextafter.pyx) defines your interface to the C++:

cdef extern from "cmath":
    float nextafter(float start, float to)

def pynextafter(start, to):
    cdef float float_start = float(start)
    cdef float float_to = float(to)
    result = nextafter(start, to)
    return result

Then a setup.py defines how to build the extension:

from distutils.core import setup
from distutils.extension import Extension
from Cython.Distutils import build_ext 

ext_modules=[
    Extension("nextafter",
        ["nextafter.pyx"],
        libraries=[],
        library_dirs=[],
        include_dirs=[],
        language="c++",
    )
]

setup(
    name = "nextafter",
    cmdclass = {"build_ext": build_ext},
    ext_modules = ext_modules
)

Make sure those are in the same directory then build with python setup.py build_ext --inplace. I hope you can see how you would add the other variations of nextafter to the extension (for doubles, etc...). Once built, you should have a nextafter.so. Fire up python in the same directory (or put nextafter.so on your path somewhere) and you should be able to call from nextafter import pynextafter.

Enjoy!

Here are five (really four-and-a-half) possible solutions.

Solution 1: use Python 3.9 or later

Python 3.9, released in October 2020, includes a new standard library function math.nextafter which provides this functionality directly: use math.nextafter(x, math.inf) to get the next floating-point number towards positive infinity. For example:

>>> from math import nextafter, inf
>>> nextafter(100.0, inf)
100.00000000000001

It's a bit easier to verify that this function really is producing the next float up if you look at the hexadecimal representation, provided by the float.hex method:

>>> 100.0.hex()
'0x1.9000000000000p+6'
>>> nextafter(100.0, inf).hex()
'0x1.9000000000001p+6'

Python 3.9 also introduces a closely related and frequently useful companion function math.ulp which gives the difference between a value and the next value away from zero:

>>> from math import ulp
>>> nextafter(100.0, inf) - 100.0
1.4210854715202004e-14
>>> ulp(100.0)
1.4210854715202004e-14

Solution 2: use NumPy

If you don't have Python 3.9 or later, but you do have access to NumPy, then you can use numpy.nextafter. For regular Python floats, the semantics match those of math.nextafter (though it would be fairer to say that Python's semantics match NumPy's, since NumPy had this functionality available long before Python did).

>>> from numpy import nextafter, inf
>>> nextafter(100.0, inf)
100.00000000000001

Solution 3: wrap C's nextafter yourself

C specifies a nextafter function in math.h (see for example section 7.12.11.3 of C99); this is exactly the function that Python >= 3.9 wraps and exposes in its math module. If you don't have Python 3.9 or later, you can either use ctypes or cffi to dynamically call C's nextafter, or alternatively write a simple Cython wrapper or Python C extension that exposes C's nextafter. The details of how to do this are already well-explained elsewhere: in @Endophage's answer to this question, and in this answer to a similar StackOverflow question (the one that this question is closed as a duplicate of).

Solution 4: bit manipulation via the struct module

If you're willing to make the (almost always safe in practice) assumption that Python is using IEEE 754 floating-point, it's quite easy to write a Python function to provide nextafter. A little bit of care is needed to get all the corner cases right.

The IEEE 754 binary floating-point formats are cleverly designed so that moving from one floating-point number to the 'next' one is as simple as incrementing the bit representation. This works for any number in the range [0, infinity), right across exponent boundaries and subnormals. To produce a version of nextUp that covers the complete floating-point range, you also need to deal with negative numbers, infinities, nans, and one special case involving negative zero. Below is a standards compliant version of IEEE 754's nextUp function in Python. It covers all the corner cases.

import math
import struct

def nextup(x):
    # NaNs and positive infinity map to themselves.
    if math.isnan(x) or (math.isinf(x) and x > 0):
        return x

    # 0.0 and -0.0 both map to the smallest +ve float.
    if x == 0.0:
        x = 0.0

    n = struct.unpack('<q', struct.pack('<d', x))[0]
    if n >= 0:
        n += 1
    else:
        n -= 1
    return struct.unpack('<d', struct.pack('<q', n))[0]

The implementations of nextDown and nextAfter then look like this. (Note that nextAfter is not a function specified by IEEE 754, so there's a little bit of guesswork as to what should happen with IEEE special values. Here I'm following the IBM Decimal Arithmetic standard that Python's decimal.Decimal class is based on.)

def nextdown(x):
    return -nextup(-x)

def nextafter(x, y):
    # If either argument is a NaN, return that argument.
    # This matches the implementation in decimal.Decimal
    if math.isnan(x):
        return x
    if math.isnan(y):
        return y

    if y == x:
        return y
    elif y > x:
        return nextup(x)
    else:
        return nextdown(x)

(Partial) solution 5: floating-point operations

If x is a positive not-too-tiny float and you're willing to assume IEEE 754 binary64 format and semantics, there's a surprisingly simple solution: the next float up from x is x / (1 - 2**-53), and the next float down from x is x * (1 - 2**-53).

In more detail, suppose that all of the following are true:

• You don't care about IEEE 754 corner cases (zeros, infinities, subnormals, nans)
• You can assume not only IEEE 754 binary64 floating-point format, but also IEEE 754 binary64 semantics: namely that all basic arithmetic operations are correctly rounded according to the current rounding mode
• You can further assume that the current rounding mode is the IEEE 754 default round-ties-to-even mode.

Then the quantity 1 - 2**-53 is exactly representable as a float, and given a positive non-subnormal Python float x, x / (1 - 2**-53) will match nextafter(x, inf). Similarly, x * (1 - 2**-53) will match nextafter(x, -inf), except in the corner case where x is the smallest positive normal value, 2**-1022.

There's one thing to be careful of when using this: the expression 2**-53 will invoke your pow from your system's math library, and it's generally not safe to expect pow to be correctly rounded. There are many safer ways to compute this constant, one of which is to use float.fromhex. Here's an example:

>>> d = float.fromhex('0x1.fffffffffffffp-1')  # 1 - 2**-53, safely
>>> d
0.9999999999999999
>>> x = 100.0
>>> x / d  # nextup(x), or nextafter(x, inf)
100.00000000000001
>>> x * d  # nextdown(x), or nextafter(x, -inf)
99.99999999999999

These tricks work right across the normal range of floats, including for awkward cases like exact powers of two.

For a sketch of a proof: to show that x / d matches nextafter(x, inf) for positive normal x, we can scale by a power of two without affecting correctness, so in the proof we can assume without loss of generality that 0.5 <= x < 1.0. If we write z for the exact mathematical value of x / d (thought of as a real number, not a floating-point number), then z - x is equal to x * 2**-53 / (1 - 2**-53). Combining with the inequality 0.5 <= x <= 1 - 2**-53, we can conclude that 2**-54 < z - x <= 2**-53, which since floats are spaced exactly 2**-53 apart in the interval [0.5, 1.0], is enough to guaranteed that the closest float to z is nextafter(x, inf). The proof for x * d is similar.