Python Testing & Debugging Beyond unittestAdvanced Python testing playbooks
Hypothesis & Fuzzing

Testing Numeric Code with Floats and NaN Edge Cases

Numerical code is where property-based testing earns its keep fastest, and also where it produces the most confusing first failures. Hypothesis's st.floats() generates exactly the values that numerical code tends to forget: nan, inf, -inf, -0.0, the largest finite double, the smallest subnormal, and numbers that differ by one unit in the last place. A mean function tested on [1.0, 2.0, 3.0] passes; the same function under Hypothesis meets [1e308, 1e308], overflows to infinity, and fails.

Those failures are almost all useful, but they need interpretation. Some expose genuine bugs — overflow in an intermediate sum, cancellation that destroys precision, a comparison that silently lets NaN through. Others expose assertions that are wrong for floating-point: exact equality where rounding makes it impossible, or relative tolerance where the true answer is zero. Working through them turns a vague sense that "floats are tricky" into an explicit contract for which values your code supports and how accurate it promises to be.

Prerequisites

Solution

Python
# stats.py
import math

def mean(xs: list[float]) -> float:
    if not xs:
        raise ValueError("mean of empty sequence")
    if any(math.isnan(x) for x in xs):
        raise ValueError("NaN in input")
    return math.fsum(xs) / len(xs)            # fsum avoids accumulated rounding
Python
# test_stats.py
import math
import pytest
from hypothesis import given, strategies as st
from stats import mean

finite = st.floats(allow_nan=False, allow_infinity=False, width=64)
moderate = st.floats(-1e150, 1e150, allow_nan=False)

@given(st.lists(moderate, min_size=1))
def test_mean_between_min_and_max(xs):
    m = mean(xs)
    assert min(xs) <= m <= max(xs)

@given(st.lists(finite, min_size=1), st.just(float("nan")))
def test_nan_is_rejected(xs, nan):
    with pytest.raises(ValueError):
        mean(xs + [nan])

@given(moderate)
def test_mean_of_constant_list(x):
    assert mean([x] * 7) == pytest.approx(x, rel=1e-12, abs=1e-300)
The special values st.floats generates A number line of double-precision floats marks the special values Hypothesis targets: negative infinity, the most negative finite value, negative and positive zero, subnormals near zero, the largest finite value, positive infinity, and NaN shown separately because it has no place on the line. Where floating-point code breaks -inf -1.79e308 subnormals -0.0 · 0.0 1.79e308 +inf NaN off the line, != itself
Every marked point is generated deliberately; code that handles the middle of the line correctly still has to decide what happens at each of them.

Why this works

st.floats() is not a uniform sampler. It deliberately mixes special values, boundary values of any bounds you give, small integers expressed as floats, and values that are adjacent in the float representation. Shrinking then moves towards simple floats — integers, then short decimals — so a failure that first appears at 1.2345e-310 often shrinks to something like 5e-324, the smallest subnormal, which names the problem immediately.

The first property, that the mean lies between the minimum and maximum, looks mathematically trivial, and naive implementations still fail it. sum(xs) / len(xs) overflows to infinity for two values near 1e308, and accumulated rounding error can push the result just outside the range for long lists of values with mixed magnitudes. math.fsum computes the exact rounded sum and fixes the rounding case; bounding the strategy to ±1e150 keeps the overflow case out of scope. That bound is a decision: the function now documents, via its test, that it supports values up to that magnitude.

Choosing tolerances deliberately

Exact equality on floats is right in fewer places than people think, and pytest.approx or math.isclose with default tolerances is right in fewer places still. Tolerances should come from the computation.

Relative tolerance expresses "accurate to about this many significant digits". A double carries roughly 15–16 significant decimal digits; each arithmetic operation can lose up to half a unit in the last place, and some operations — subtracting nearly equal numbers — lose far more. For a short computation, rel_tol=1e-12 is a reasonable start. For iterative algorithms or matrix operations, the achievable accuracy depends on the condition number of the problem, and a looser tolerance with a comment explaining why is better than a tight one that flakes.

Absolute tolerance handles results near zero. When the true answer is 0.0, any computed 1e-17 has infinite relative error, so a purely relative comparison always fails. Set abs_tol to the size of rounding noise you expect for inputs of the magnitude being tested. Both tolerances together cover the full range: relative for large results, absolute for tiny ones.

Relative and absolute tolerance together A chart plots the allowed error against the magnitude of the expected result. The relative tolerance line grows with magnitude, while the absolute tolerance is a flat floor. Near zero the absolute floor dominates; for large results the relative line dominates. math.isclose accepts anything below the higher of the two. Allowed error across magnitudes |expected result| allowed error rel_tol x |result| abs_tol floor near zero
Without the absolute floor, every property whose true answer is zero fails on rounding noise alone.

NaN, infinities and signed zero as contract decisions

For each special value, the function should do one of three things, and the tests should say which.

Reject it. Raise a clear exception, as the mean above does for NaN. Test with a strategy that always includes the value, and pytest.raises. Rejection is usually right for inputs that indicate upstream corruption — a NaN in a price list is a bug somewhere else.

Propagate it. Return NaN or infinity following IEEE semantics, as NumPy does. Test that the output is NaN with math.isnan rather than equality, and that finite inputs never produce NaN. Propagation is right for low-level numerical libraries whose callers expect IEEE behaviour.

Exclude it from scope. Configure the strategy never to generate it and document the precondition. This is the weakest choice, because production inputs are not bound by your strategy, and it is appropriate only where an earlier validation layer guarantees the value cannot arrive.

Signed zero deserves a special mention because it is invisible to equality: -0.0 == 0.0 is true, but math.copysign(1, -0.0) is -1.0, and 1 / -0.0 raises in Python while producing -inf in NumPy. Code that uses the sign of a result — angle calculations, direction of rounding — can behave differently for the two zeros. If that matters, compare with math.copysign or struct.pack rather than ==.

A worked failure: variance that goes negative

The classic numerical bug that Hypothesis finds within seconds is a one-pass variance computed as the mean of squares minus the square of the mean:

Python
def variance(xs):
    n = len(xs)
    return sum(x * x for x in xs) / n - (sum(xs) / n) ** 2

The property is simple — variance is never negative — and it fails on a shrunk example along the lines of [1e8, 1e8 + 1]. Both terms are about 1e16, their true difference is 0.25, and the subtraction of two nearly equal large numbers leaves only a few correct bits. With slightly different values the result comes out as a small negative number, which is mathematically impossible and breaks any later sqrt.

The fix is not a tolerance on the test. It is a numerically stable algorithm: Welford's online method, or a two-pass computation that subtracts the mean before squaring. After the change, the non-negativity property passes, and a second property — that the result matches statistics.pvariance (which uses exact fractions internally) to rel_tol=1e-9 — confirms the accuracy improved rather than just the sign.

That sequence is the typical shape of floating-point property testing. An obviously true mathematical property fails, the shrunk input points at a magnitude where precision runs out, and the fix is a better algorithm rather than a looser assertion. Loosening the assertion would have hidden a bug that production data — timestamps in nanoseconds, prices in large units — would eventually hit. It is also worth keeping the shrunk input as a permanent @example on the property: the values that exposed the cancellation are exactly the ones a future refactor back to the one-pass formula would fail on, and pinning them means that regression is caught on every run rather than only when the random search happens to wander back into the same region of large, nearly equal numbers.

Cancellation in a one-pass variance Two large nearly equal quantities, the mean of squares and the square of the mean, are each about ten to the sixteen. Subtracting them leaves only a few significant bits, so the computed variance can be negative. A two-pass or Welford algorithm subtracts the mean first and keeps full precision. Subtracting two huge numbers loses the answer mean(x²) ≈ 1.0000000200e16 mean(x)² ≈ 1.0000000200e16 − -2.0 (wrong) true value 0.25 Welford / two-pass subtract mean first Property "variance ≥ 0" finds it; a tolerance would hide it.
The shrunk counterexample names the magnitude where precision runs out, which points straight at the unstable subtraction.

Edge cases and failure modes

  • float32 code tested with float64 values. NumPy code that stores float32 rounds every generated double. Use st.floats(width=32) so generated values are exactly representable.
  • Overflow in intermediate values. Squares, sums and products overflow long before the inputs reach the float limit. Bound the strategy by the square root of the limit for anything that squares.
  • Catastrophic cancellation. a - b for nearly equal a and b keeps few significant digits. Properties comparing such results need tolerances relative to the inputs, not the result.
  • hypothesis.extra.numpy arrays. arrays(np.float64, shape, elements=st.floats(...)) passes the element strategy through; forgetting elements generates the full float range including NaN.
  • Sorting with NaN. sorted with NaN produces order that depends on input position. Exclude NaN from any property involving ordering unless NaN ordering is the point.

Frequently Asked Questions

Why does my float test fail with NaN != NaN? IEEE 754 defines NaN as unequal to everything, including itself, so any equality assertion on a NaN result fails. Either exclude NaN from the strategy, decide the function should reject NaN and test that, or compare with math.isnan on both sides.

What tolerance should I use with math.isclose? Use a relative tolerance based on how many operations the computation performs — 1e-9 is a common start for short double-precision computations — plus an absolute tolerance for results near zero, where relative error is meaningless.

Should I generate subnormal floats? Yes, unless the code explicitly does not support them. Hypothesis generates subnormals by default because they are a common source of precision loss and performance surprises. Restrict with allow_subnormal=False only when that is a documented limitation.

← Back to Advanced Property-Based Testing