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

floating-point, math, python

Solution

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 `float`s, 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.

Problem

Does Python provide a function to get the floating-point value that results from incrementing the least significant bit of an existing floating-point value? I'm looking for something similar to the `std::nextafter` function that was added in C++11.

Original source

Related problems