Correctness, seeding, benchmarks and a working CUDA backend

.github/workflows/python-distro.yml

@@ -46,8 +46,8 @@ jobs:
             -   name: Build C++ tests
                 run: cmake --build build --target fastdist_tests --parallel
 
-            # RNG tolerances are sized from the estimator standard error (>=7 sigma),
-            # so a failure here is a real regression rather than a flake.
+            # The RNG tests are seeded, so a failure here reproduces on every run
+            # with the same toolchain rather than being a flake.
             -   name: Run C++ tests
                 run: ctest --test-dir build --output-on-failure
 

CHANGELOG.md

@@ -12,16 +12,61 @@ call in `CMakeLists.txt`.
 
 ### Fixed
 
-- `binomial_cdf_scalar`, `poisson_cdf_scalar` and `negative_binomial_cdf_scalar` returned the raw sum of
-  PMF terms, which accumulates rounding error and could exceed `1.0`. For the binomial this also made the
-  CDF non-monotonic, because `x >= n` short-circuits to exactly `1.0` while the values just below it did
-  not. All three now clamp to `1.0`.
+- `beta_cdf_scalar` was wrong at every point, not only at the extremes: 0.0015 against a true 0.1143 for
+  Beta(2, 5) at x = 0.1, and values outside [0, 1] such as -147 for Beta(0.01, 0.01) at x = 0.5. It is now
+  the modified-Lentz continued fraction with the standard reflection, and agrees with SciPy to ~1e-12.
+- `gamma_cdf_scalar` and `chi_square_cdf_scalar` returned probabilities above 1.0 (Gamma(1.5, 1).cdf(2.5)
+  gave 1.000498). A unary minus applied to an unsigned loop index wrapped to 2^32 - i inside the continued
+  fraction. Separately, the series stopped at 100 iterations and silently truncated for large shapes (off
+  by 0.16 at alpha = 10000); the ceiling is now 1000.
+- `beta_pdf_scalar`, `gamma_pdf_scalar` and `chi_square_pdf_scalar` returned `nan` or `inf` once a shape
+  parameter passed ~171 (k above ~342 for chi-square), because the normalising constants overflowed
+  `std::tgamma`. All three are now evaluated in log space.
+- `negative_binomial_pmf_scalar` returned `inf` at k = 170 and `nan` beyond it, taking the CDF with it, for
+  the same reason. It is now evaluated in log space.
+- `binomial_cdf_scalar`, `poisson_cdf_scalar` and `negative_binomial_cdf_scalar` could exceed `1.0`
+  through accumulated rounding, which also made the binomial CDF non-monotonic. All three now clamp to
+  `1.0`.
+- `beta_sample` did not validate its parameters, and `std::gamma_distribution` has undefined behaviour for
+  a non-positive shape. It now returns `nan` for invalid input, like every other continuous sampler.
+- Python setters: `Beta.beta` and `Binomial.p` raised `TypeError` for every value, valid or not;
+  `DiscreteUniform.b` recursed until `RecursionError`; `DiscreteUniform.a` changed the attribute's type to
+  `float`; and the `Uniform` and `DiscreteUniform` setters accepted a bound that violated `a < b`, after
+  which every method silently returned `nan`.
+- `Utils.law_of_total_probability` rejected scalar arguments despite its signature. Scalars are now
+  treated as a one-element partition, and sequences of different lengths raise `ValueError`.
+- An `ImportError` raised while loading a distribution module, such as a missing numpy, was reported as a
+  missing C++ core. The original exception is now chained, and the message says how to build the
+  extension.
+- The C++ RNG tests failed in about 7.7% of runs because their tolerances sat near 2σ of the estimator's
+  own noise (#2).
+- The CUDA backend did not compile on Windows. `nvcc` 12.x's front end crashes on MSVC's C++20
+  standard-library headers, and every `.cu` file was compiled a second time into the Python module
+  target, which built as C++20. CUDA sources are now compiled once, as C++17.
+- The CUDA extension could not be imported on Windows, because it depended on `cudart64_*.dll` and
+  Python does not search `PATH` for extension dependencies. Once imported, its first GPU call crashed
+  with an access violation, because the wrapper released the GIL before touching the input and output
+  arrays. The CUDA runtime is now linked statically, and the GIL is released only around the device
+  work.
+- `normal_cdf` lost all relative precision in the lower tail: Phi(-8) was off by 1.8%, and Phi(-10)
+  returned exactly 0 instead of 7.6e-24. `exponential_cdf` did the same for small arguments,
+  returning 0 at x = 1e-17. Both now use `erfc` and `expm1` respectively, on the CPU and GPU paths.
 - `setup.py` no longer hardcodes the `Visual Studio 17 2022` CMake generator. CMake selects the newest
   Visual Studio present, so builds work on machines with a different version installed. Set
   `CMAKE_GENERATOR` to pin one.
 
 ### Added
 
+- `fastdist.seed(value)` and `fastdist.seed_from_entropy()`, with C++ equivalents `seed_rng` and
+  `seed_rng_from_entropy` in `fastdist/math/rng.h`. Every sampler now draws from one shared thread-local
+  Mersenne Twister, seeded through `std::seed_seq`, and any signed 64-bit value is accepted. A seed
+  reproduces a run on one platform and toolchain, and applies to the calling thread only.
+- A benchmark suite under `benchmarks/`, and `BENCHMARKS.md` as a performance log generated from its
+  recorded results. It compares against SciPy, checks numerical agreement before timing, flags regressions
+  against the run's measured noise, and times the CUDA paths when they are built.
+- Type stubs for the compiled extension (`_fastdist.pyi`), with a CI check that keeps them in step with the
+  bindings.
+- A PyPI release workflow using Trusted Publishing, with a TestPyPI dry-run option.
 - Cross-platform wheel building in CI via `cibuildwheel` — Linux x86_64, Windows AMD64, and macOS
   x86_64 + arm64, for CPython 3.10 through 3.14.
 - An install-from-sdist check in CI, exercising the source path an end user takes on any platform without
@@ -35,6 +80,22 @@ call in `CMakeLists.txt`.
 
 ### Changed
 
+- **Performance.** The Poisson, binomial and negative binomial CDFs sum their terms by recurrence instead
+  of re-deriving each one, and the batch paths hoist parameter validation and loop-invariant terms out of
+  their loops. On the reference machine in `BENCHMARKS.md`, `poisson_cdf` over 100k values went from
+  43.5 ms to 1.3 ms (from 0.15x SciPy's speed to 4.8x), `normal_logpdf` got 80% faster, and the uniform and
+  normal CDF paths got 22–29% faster.
+- `Utils.sigmoid` is scalar-only and raises a `TypeError` naming `Utils.sigmoid_cpu` when given a sequence.
+  Its annotation previously advertised sequences, which never worked.
+- Annotations use `typing.SupportsFloat` rather than `numbers.Real`, which mypy cannot check, so correct
+  code such as `Normal(0.0, 1.0)` no longer reports errors. The package now type-checks cleanly.
+- The `exponential` and `poisson` bindings name their rate keyword `lambda_`. The old name, `lambda`, is a
+  Python keyword and could never be passed by name.
+- The sampling tests are seeded and deterministic, with tolerances at about 5× the estimator's standard
+  error, and CI no longer retries failed C++ tests.
+- `requirements.txt` is now `requirements-dev.txt`. Runtime dependencies are declared only in the package
+  metadata.
+- Package metadata moved from `setup.py` into `pyproject.toml`.
 - `NDEBUG` is undefined for the `fastdist_tests` target, so its `assert()`-based checks stay live in Release
   builds. They were previously compiled away, meaning the suite reported success without testing anything.
 - Repository layout: bindings moved from `python/bindings` to `src/bindings`; tests split into `tests/cpp`
@@ -46,12 +107,16 @@ call in `CMakeLists.txt`.
 
 ### Known issues
 
-- `negative_binomial_pmf_scalar` returns `inf` for `k` around 169 and `NaN` beyond it, because the binomial
-  coefficient is computed with `std::tgamma`, which overflows above ~171. Computing it in log space via
-  `std::lgamma` is the fix. `beta.cpp` and `gamma.cpp` use `tgamma` similarly.
-- Three RNG test assertions have tolerances at roughly 2σ and fail a few percent of runs. CI absorbs this
-  with `ctest --repeat until-pass:3`.
-- The samplers seed `std::mt19937` from a single 32-bit `random_device` word.
+- The name `fastdist` belongs to an unrelated project on PyPI, so this package cannot be published under
+  it. The release workflow refuses to upload until the distribution is renamed.
+- Sampling draws one variate per call, which makes bulk generation 30–100x slower than numpy. There is no
+  batch sampling entry point yet.
+- The gamma CDF's iteration ceiling covers shape parameters up to roughly 20000. Beyond that the result
+  degrades without warning.
+- The `*_cpu` bindings do not expose the `step_size` default that the C++ headers declare, and `step_size`
+  is a `double` for the continuous distributions but an `int` for the discrete ones.
+- CI compiles the CUDA backend (on Linux, in the type-stub job) but has no GPU runner, so the CUDA
+  kernels are not exercised there.
 
 ---
 

CMakeLists.txt

@@ -21,6 +21,11 @@ option(FASTDIST_ENABLE_CUDA "Enable CUDA backend" OFF)
 
 if (FASTDIST_ENABLE_CUDA)
     enable_language(CUDA)
+    # CUDA is compiled as C++17 even though the rest of the project is C++20:
+    # nvcc 12.x's front end (cudafe++) crashes on the MSVC standard library
+    # headers in C++20 mode. Set globally so no target can pick up C++20.
+    set(CMAKE_CUDA_STANDARD 17)
+    set(CMAKE_CUDA_STANDARD_REQUIRED ON)
     message(STATUS "CUDA backend enabled")
 else ()
     message(STATUS "CUDA backend disabled")
@@ -159,39 +164,6 @@ if (FASTDIST_ENABLE_CUDA)
             src/cuda/utils/manhattan_distance.cu
             src/cuda/utils/cosine_similarity.cu
     )
-    target_sources(_fastdist PRIVATE
-            src/cuda/normal/pdf.cu
-            src/cuda/normal/logpdf.cu
-            src/cuda/normal/cdf.cu
-            src/cuda/normal/mgf.cu
-            src/cuda/normal/cgf.cu
-
-            src/cuda/poisson/cdf.cu
-            src/cuda/poisson/cgf.cu
-            src/cuda/poisson/mgf.cu
-            src/cuda/poisson/pmf.cu
-
-            src/cuda/bernoulli/cdf.cu
-            src/cuda/bernoulli/cgf.cu
-            src/cuda/bernoulli/mgf.cu
-            src/cuda/bernoulli/pmf.cu
-
-            src/cuda/exponential/cdf.cu
-            src/cuda/exponential/cgf.cu
-            src/cuda/exponential/mgf.cu
-            src/cuda/exponential/pdf.cu
-
-            src/cuda/uniform/cdf.cu
-            src/cuda/uniform/cgf.cu
-            src/cuda/uniform/mgf.cu
-            src/cuda/uniform/pdf.cu
-
-            src/cuda/utils/sigmoid.cu
-            src/cuda/utils/logit.cu
-            src/cuda/utils/euclidean_distance.cu
-            src/cuda/utils/manhattan_distance.cu
-            src/cuda/utils/cosine_similarity.cu
-    )
     set_target_properties(fastdist_core PROPERTIES
             CUDA_SEPARABLE_COMPILATION OFF
             CUDA_STANDARD 17
@@ -200,8 +172,12 @@ if (FASTDIST_ENABLE_CUDA)
     )
 
     # CUDA Runtime Linking
-    find_package(CUDAToolkit REQUIRED) # Finds Toolkit
-    target_link_libraries(fastdist_core PUBLIC CUDA::cudart) # Links Runtime library
+    # The runtime is linked statically so the extension is self-contained.
+    # Linked dynamically, _fastdist depends on cudart64_*.dll, which Python 3.8+
+    # on Windows will not find through PATH -- the module fails to import unless
+    # every caller adds the toolkit to the DLL search path first.
+    find_package(CUDAToolkit REQUIRED)
+    target_link_libraries(fastdist_core PUBLIC CUDA::cudart_static) # Links Runtime library
 
     target_include_directories(fastdist_core PRIVATE src)
     target_include_directories(_fastdist PRIVATE src)
@@ -280,5 +256,5 @@ endforeach ()
 
 # Link CUDA to Python module if enabled
 if (FASTDIST_ENABLE_CUDA)
-    target_link_libraries(_fastdist PRIVATE CUDA::cudart)
+    target_link_libraries(_fastdist PRIVATE CUDA::cudart_static)
 endif ()

README.md

@@ -86,6 +86,40 @@ logit, Euclidean/Manhattan distance, cosine similarity, coefficient of variation
 
 ---
 
+## Validation and error handling
+
+One rule at each layer, the same for every distribution.
+
+**Parameters are checked when you construct a distribution.** A value that cannot describe a
+distribution raises `ValueError`; the wrong type raises `TypeError`. Nothing is constructed, so no
+later call can quietly return nonsense. Non-finite parameters are refused too -- `nan` passes every
+range comparison, so it is rejected explicitly.
+
+```python
+Normal(0.0, -1.0)             # ValueError: sigma must be positive
+Normal(0.0, float("nan"))     # ValueError: sigma must be finite
+Normal(0.0, "1.0")            # TypeError: sigma must be a real number
+```
+
+**A non-finite input is not an error.** `x = nan` or `+/-inf` yields `nan`, matching the C++ core and
+numpy's elementwise behaviour, so one bad value in an array does not abort the whole call. A string
+is still a `TypeError`, even though numpy would happily read `"0.5"` as a number.
+
+```python
+Normal(0.0, 1.0).pdf(float("nan"))          # nan
+Normal(0.0, 1.0).pdf([0.0, float("inf")])   # array([0.3989..., nan])
+Normal(0.0, 1.0).pdf("0.5")                 # TypeError
+```
+
+**Outputs stay inside their mathematical range.** CDFs are clamped to `[0, 1]`, so accumulated
+rounding cannot hand back `1 + 1e-16` to code that treats the result as a probability.
+
+**The C++ API has no exceptions**, so it signals invalid parameters by return value: `NaN` from
+anything returning a real number, `-1` from the integer samplers, and `INT_MIN` from
+`discrete_uniform_sample`.
+
+---
+
 ## Reproducible sampling
 
 Every `*_sample()` call draws from one shared Mersenne Twister engine. Seeding it makes a run reproducible: