statistics.NormalDist.cdf returns incorrect probability when x - mu overflows

[augusto-rehfeldt]

Bug report

Bug description:

Documented behaviour: Public API documentation for NormalDist.cdf: "Cumulative distribution function. P(X <= x)"; NormalDist describes mu and sigma as the mean and standard deviation.

Expected: 0.9772498680518209, within absolute tolerance 1e-15

Actual: 1.0

import math
import statistics

a = 1e308
mu, sigma, x = -a, a, a
if not (all(map(math.isfinite, (mu, sigma, x))) and 0 < a <= 1.2e308):
    print("REFUTATION REJECTED: invalid input")
else:
    actual = statistics.NormalDist(mu, sigma).cdf(x)
    z = x / sigma - mu / sigma  # Avoid overflow in x - mu.
    expected = (1 + math.erf(z / math.sqrt(2))) / 2
    if not abs(actual - expected) <= 1e-15:
        print("REFUTATION CONFIRMED:", (mu, sigma, x), actual, expected)
    else:
        print("REFUTATION REJECTED: agrees within absolute tolerance 1e-15")

Output on Python 3.14.6 (Windows-11-10.0.26220-SP0), standard library statistics:

REFUTATION CONFIRMED: (-1e+308, 1e+308, 1e+308) 1.0 0.9772498680518209

This report was found and written by an automated property-testing tool I run (bugforge). The reproducer above was executed and its output is pasted unedited; no person reviewed the report before it was filed. The search script is in https://github.com/augusto-rehfeldt/bugforge-results/tree/main/statistics-20261003-045701-c1

CPython versions tested on:

3.14

Operating systems tested on:

Windows

[ZeroIntensity]

See #158748 (comment).