math.nextafter with steps as a pure Python prototype

nextafter.py

import math
import struct
from functools import reduce
from itertools import repeat

from hypothesis import assume, example, given, note
from hypothesis import strategies as st


def float_from_integer(integer):
    return struct.unpack("!d", struct.pack("!Q", integer))[0]


def float_to_integer(float):
    return struct.unpack("!Q", struct.pack("!d", float))[0]


def nextafter_reduce(start, goal, steps):
    return reduce(math.nextafter, repeat(goal, steps), start)


def nextafter_loop(start, goal, steps):
    for _ in range(steps):
        start = math.nextafter(start, goal)
    return start


# TODO: Think about supporting negative steps?
def my_nextafter_steps(start, goal, steps):
    # We only check for zero steps, so that
    # nextafter(start, nan, 0) returns start,
    # instead of giving nan.
    if steps == 0:
        return start
    if math.isnan(start) or math.isnan(goal):
        return math.nan

    if math.copysign(1, start) < 0:
        return -my_nextafter_steps(-start, -goal, steps)
    assert start >= 0

    if math.copysign(1, goal) < 0:
        i = float_to_integer(start)
        if steps > i:
            return -my_nextafter_steps(0.0, -goal, steps - i)
        else:
            # Setting goal to 0 ain't necessary, but makes the invariants
            # further down simpler to explain
            goal = 0.0
    assert start >= 0 and goal >= 0 and steps > 0

    if start < goal:
        steps = steps
        cmp = min
    elif start == goal:
        return start
    else:
        steps = -steps
        cmp = max
    return float_from_integer(
        cmp(float_to_integer(goal), steps + float_to_integer(start))
    )


@example(start=1.0, goal=math.nextafter(1.0, math.inf), steps=10)
@given(
    start=st.floats(),
    goal=st.floats(),
    steps=st.integers(min_value=0, max_value=1_000_000),
)
def test(start, goal, steps):
    soll = nextafter_reduce(start, goal, steps)
    ist = my_nextafter_steps(start, goal, steps)
    note(f"{soll} != {ist}")
    assert math.isnan(soll) and math.isnan(ist) or soll == ist


@given(
    start=st.floats(),
    goal=st.floats(),
    steps1=st.integers(min_value=0, max_value=1_000_000),
    steps2=st.integers(min_value=0, max_value=1_000_000),
)
def test_goal2(start, goal, steps1, steps2):
    goal1 = nextafter_reduce(start, goal, steps1)
    soll = nextafter_reduce(start, goal1, steps2)
    ist = my_nextafter_steps(start, goal1, steps2)
    note(f"{soll} != {ist}")
    assert math.isnan(soll) and math.isnan(ist) or soll == ist


if __name__ == "__main__":
    test()
    test_goal2()

Inspired by https://stackoverflow.com/questions/44991042/binary-search-bisection-for-floating-point-numbers