Correctness, seeding, benchmarks and a working CUDA backend

[ghosteau]

Adds 20 commits on top of develop, which already contains repo-organization. develop is merged in, so CI here tests the actual merged result. Each commit is self-contained and independently revertible; the CUDA changes are kept in their own commits.

Intended for a careful read-through before the next release.

Correctness — 42 known bugs to 0

The suite started this work with 42 xfail(strict) markers documenting real defects. All are fixed and the markers are gone.

Three CDFs returned wrong answers, two of them probabilities above 1.0.

• beta_cdf was wrong at every point: 0.0015 against a true 0.1143 for Beta(2,5) at x=0.1, and -147 for Beta(0.01,0.01) at x=0.5. The series had an inverted coefficient ratio and normalised by Γ(a+1) instead of B(a,b). Replaced with the modified-Lentz continued fraction.
• gamma_cdf and chi_square_cdf exceeded 1.0 — Gamma(1.5,1).cdf(2.5) gave 1.000498 — because -i * (i - a) applied a unary minus to an unsigned loop index, wrapping it to 2³²−i.
• MAX_ITER was 100, silently truncating the gamma series for large shapes (off by 0.16 at α=10000).

Densities overflowed. beta_pdf, gamma_pdf and chi_square_pdf returned NaN once a shape passed ~171 (tgamma overflow); negative_binomial_pmf did the same at k=170. All now evaluated in log space.

Tail probabilities were zero. normal_cdf(-10) returned exactly 0 against a true 7.6e-24, and exponential_cdf(1e-17) returned 0. Both now use erfc/expm1 where the naive form cancels, and the cheaper form elsewhere. This costs normal_cdf 63–70%, recorded in BENCHMARKS.md.

Five Python defects, including a setter that recursed until RecursionError and setters that let Uniform reach a=10, b=-10, a state the constructor rejects, after which every method silently returned NaN.

Validated against SciPy over dense parameter grids: every CDF agrees to ~1e-12, and none is outside [0,1], non-monotonic or NaN.

Seeding

fastdist.seed(value) / seed_from_entropy(), plus seed_rng in C++. All samplers share one thread-local Mersenne Twister, seeded through std::seed_seq, accepting any signed 64-bit value. Two limits are documented at the API: seeding is per-thread, and a seed reproduces a run on one platform/toolchain, not across libstdc++/libc++/MSVC.

This also made the flaky RNG tests deterministic (#2): their tolerances sat near 2σ of the estimator's own noise and failed ~7.7% of runs. They are now seeded and sized at ~5σ. 150 consecutive full ctest sweeps: 0 failures.

Validation — one contract

tests/python/test_validation.py states the rule and holds all twelve classes to it. Fixed along the way: seven classes accepted NaN parameters (every comparison against NaN is false, so range checks waved it through), Normal(0,1).pdf("0.5") returned array([0.352]) because numpy reads the string as data, and Gamma/ChiSquare returned 0.0 for -inf where every other density returns NaN. The contract is documented in the README, with every example executed before being written down.

Benchmarks — and an honesty caveat worth reading

benchmarks/ plus BENCHMARKS.md, an evidence log generated from recorded JSON, never typed by hand. examples/release_benchmarks.ipynb renders it.

Two baselines, deliberately. Against scipy.stats we are 2.3–16x faster. Against the primitives underneath (numpy, scipy.special) we are roughly at parity — behind at 100k, ahead on six of ten at 1M. Most of the scipy.stats margin is that library's generic dispatch machinery, not faster arithmetic. Quoting only the first number would oversell the library, so both are recorded for every case.

Real wins: poisson_cdf 3.6x faster than scipy.special.pdtr (recurrence beats a general implementation); scalar calls ~60x cheaper than scipy.stats. Real losses, recorded: bulk sampling is 30–100x slower than numpy, and bernoulli_pmf is 0.16x against one np.where.

CUDA — now builds and wins

It previously did not compile on Windows, could not be imported once built, and crashed on first call once imported.

• cudafe++ crashes on MSVC's C++20 headers; CUDA is now pinned to C++17, and each kernel is compiled once instead of twice.
• The runtime is linked statically, so the extension no longer needs the toolkit on a DLL path.
• The wrapper released the GIL before touching the input and output arrays; it now releases only around the device work.

Measured warm, at n=1M: normal_cdf 4.6x, normal_pdf 2.6x, exponential_pdf 2.1x over our CPU path. normal_logpdf and uniform_pdf still lose. An earlier entry concluded the GPU lost everywhere; that was measured on an idle, downclocked card and is marked withdrawn in the log rather than deleted.

Verification

Linux / Windows / macOS wheels + pytest, sdist install	pass
Ubuntu C++ build, ctest, pytest, mypy	pass
Linux CUDA build (type-stub job)	pass
Python suite	1879 passed
xfail / known-bug markers	0 (was 42)

Known and deliberately not done

• The name fastdist is taken on PyPI by an unrelated library. The release workflow refuses to publish until it is renamed. This is the one hard blocker for a release.
• CUDA auto-dispatch thresholds default to 100k for every function. Correct for three of the five measured, wrong for normal_logpdf and uniform_pdf. config.auto_tune cannot express this: it searches from 500k and cannot answer "never".
• Batch sampling does not exist; every variate crosses the Python/C++ boundary individually.
• The CUDA comment-only commit was not compiled in CI on Windows; only clang-format checked it.

🤖 Generated with Claude Code