All notable changes to RandLib are documented here. The project follows Semantic Versioning.

[2.0.0] - Unreleased

Mathematical reference

Added a current mathematical reference for every public distribution, with TeX definitions, moments, transforms, sampling and numerical-method notes rendered through MathJax in Doxygen. Seventy focused distribution pages replace the monolithic category chapters, family landing pages provide card navigation, and 54 density, mass and CDF plots from the verified legacy set are published with responsive styling. The strict docs target validates page structure, card links, plot assets, TeX delimiters and rendered titles.

Stable distributions

Extended StableRand's valid exponent range from [0.1, 2] to (0, 2].
Added Nolan's S0 parameterization alongside the source-compatible S1 default. Constructors, location accessors, setters, names, support endpoints, probability functions, characteristic functions and sampling consistently use the selected convention, including the special logarithmic location conversion at alpha = 1.
Implemented stable densities, distribution functions and quantiles throughout the former skewed near-alpha-one gap. A cancellation-free S0 characteristic function and alpha-one sampling limit keep adjacent representable exponents continuous; 40-digit mpmath values cover both sides.
Added McCulloch quantile fitting as a fast robust estimator and initializer for constrained stable maximum likelihood. Fits are reported in the continuous S0 parameterization.
Added alpha-one density and probability tail expansions through eighth order. Geometric-stable mixture evaluation now covers its full historical tail grid in seconds instead of minutes, with focused 40-digit mpmath references on both sides.
Added the log-stable distribution as the exponential transform of a stable law, with transformed density, direct tails, sampling, quantiles, median, correct boundary/interior modes, S0/S1 parameterization and stable-family fitting. Exact Laplace-exponent moments cover totally left-skewed laws and the Gaussian case; other non-Gaussian skewnesses report divergent moments.
Implemented the log-stable characteristic function (#5) by half-period Fourier integration bounded by the survival function, with the mass near zero entering through the distribution function. Independent references cover the Gaussian case, log-Cauchy and log-Levy laws, and the convergent moment series of totally left-skewed laws.
Removed the last stable tolerance override (#1): far upper-tail Landau quantiles near 1e10 now meet the default reference tolerance.

Distributions

Added the Tier 2 StableCountRand<Real> location-scale distribution for 0 < alpha < 1, with Wright/one-sided-stable density evaluation, direct reciprocal-stable probability tails, TOMS 748 quantiles, log-gamma moments, standard-URBG inverse sampling, an exact half-stable Gamma path, constrained maximum-likelihood fitting, and independent 40-digit mpmath reference data.
Added the Tier 1 noncentral t distribution with centered Benton–Krishnamoorthy probability series, direct survival probabilities, recurrence and integral densities, transformed sampling, quantiles, moments, a conditional-normal characteristic function and constrained maximum- likelihood fitting. Saddle-split integration and a small-denominator expansion retain accuracy for large noncentrality and representable extreme tails; a 243-value 40-digit mpmath table and focused exact-series regressions provide evidence.
Made its centered series and moments cancellation-resistant on platforms where long double is binary64, using Loader's Poisson saddle-point form, stable log-beta ratios, compensated summation and grouped central-moment identities. Exact probability levels at zero now invert to zero.
Evaluated its moderate-argument density directly from the hypergeometric representation, avoiding cancellation between two probabilities on binary64-only platforms.

Accuracy and robustness

Evaluated the log-normal characteristic function (#6) by the exact saddle-point representation of Asmussen, Jensen and Rojas-Nandayapa with a complex Lambert W function, keeping relative accuracy where it is tiny, for example 8.3e-29 at t = 10 for LogNormal(3, 0.01).
Made the incomplete beta function recover 1 - x from its logarithm, so the beta survival function no longer returns one below x = 1e-16 and the beta-prime distribution function no longer rounds to one for large x. Shapes above 4096 use a modified Lentz continued fraction: binomial laws with n = 2e9 take microseconds instead of seconds, with relative errors of 1e-13 instead of 5e-5; for shapes up to 1e12 the measured worst case is 1e-8, just above the convergence boundary of extremely skewed laws. Arguments near one are never evaluated by slowly converging series, which bounds the cost of noncentral t fitting at extreme trial parameters.
Computed binomial distribution functions from a cancellation-free log-beta value; the MSVC bulk tolerance override is removed (#12).
Started the characteristic-function walk where the mass begins and bounded it by the survival function: Frechet and Rice laws concentrated far from the end of their support no longer return zero. Expectations without closed forms integrate on both sides of the median with a bulk-scaled map, fixing zero characteristic functions of narrow Gumbel and log-normal laws.
Added QUADPACK roundoff detection to adaptive quadrature, making Stable Count probabilities near alpha = 1/2 about a thousand times faster.
Evaluated Marcum P for large noncentrality by a saddle-centered lower-tail mixture, fixing NaN Rice distribution functions, and made both Marcum functions complementary to rounding.
Added Hankel expansions to the logarithms of overflowing or underflowing Bessel functions and a large-order series for underflowing I, fixing Rice moments and Skellam tail probabilities; logBesselK now uses the symmetry in the order instead of returning NaN.
Evaluated the exponentially modified Gaussian density and tails through erfcx, fixing NaN lower tails for large rate * sigma; Irwin-Hall functions on the nearer half of the symmetric support, fixing negative survival probabilities; Planck and zeta survival functions without subtracting from one; von Mises tails below the series threshold by direct integration; and generalized extreme-value moments near the Gumbel limit without cancellation.
Fixed negative-binomial sampling for small shapes, where the mixing gamma variate underflows.
Made the Brent minimizer behind the generic Mode() use a relative tolerance: an absolute 1e-8 never ended the iteration once doubles near the minimum were spaced farther apart. Stable Count modes are searched in log(x), where small exponents are no longer flat over many orders.
Fixed GetSampleStatistics for more than 10 000 observations, the sort check, cell set and degrees of freedom of PearsonChiSquaredTest, and empty-sample handling in the member goodness-of-fit tests. Bulk evaluation helpers no longer write past a shorter output span, and QuantileFunction and QuantileFunction1m are const.
fit now rejects Method::Bayes for normal and Poisson models instead of returning a likelihood estimate; conjugate updates use fit_bayes.
Added a statistical check that unified maximum-likelihood fits are local maxima of the likelihood, completing the unit-testing checklist of #12, and dev-unit and dev-statistical test presets.
Documented why binomial sampling stays exact instead of switching to a Poisson approximation (#11).

Numerical kernel

Added type-generic tanh-sinh quadrature for finite intervals with integrable endpoint singularities. Endpoint-aware callables receive the accurately computed distance to the nearer boundary, avoiding cancellation when the abscissa itself rounds to that boundary.
Added QAWF-style semi-infinite Fourier quadrature with half-period subdivision and Wynn epsilon acceleration. Heavy-tailed Pareto, Frechet, inverse-gamma, Fisher F and beta-prime characteristic functions, plus the large-variance log-normal case, now meet the standard 1e-8 absolute target; all twelve former tolerance overrides were removed.
Added direct Marcum-Q evaluation through a saddle-centered Poisson mixture. Noncentral chi-squared survival probabilities now retain relative accuracy through representable tails instead of subtracting the CDF from one; all three former tolerance overrides were removed.
Added the Tier 2 Rice distribution with direct density, CDF and survival functions, transformed- normal sampling, closed-form mean and variance, and 40-digit mpmath reference data.

Build and tooling

Avoided a GCC 12 false positive in BFGS identity-matrix construction and replaced std::jthread in the parallel-stream example for Apple libc++ versions that do not provide it.
Reduced the stable-law tangent argument, made quarter-turn tangent and inverse-tangent values exact, and evaluated the finite totally skewed integration-endpoint limit analytically. Short-tail quantiles now agree on platforms where long double is binary64.

Documentation

Rewrote the README around distribution fitting and heavy-tail work, with a quick start and an explicit comparison with the standard library, Boost.Math and StatsLib.
Added guides for fitting, parallel stream ownership, stable-law parameterization and limitations, the numerical evidence policy, packaging, and every v1-to-v2 breaking change.
Added a strict Doxygen build with the pinned Doxygen Awesome theme and a GitHub Pages deployment workflow. Redundant conversion-era @fn directives and incomplete public parameter comments were cleaned up so documentation warnings fail the build.
Standardized statistical formulas on square-bracket operator notation, with expectations written as \mathbb{E}[\,\cdot\,], and made the documentation build reject regressions.
Retired the generated legacy PDF. Its LaTeX source and figures remain under docs/legacy/ as a clearly marked 1.x mathematical reference.

Sampling efficiency

Replaced samplers whose expected work grew with a parameter by generators with bounded expected work (Devroye, 1986). Sampled streams change for the affected parameter regions; categorical streams do not.
- CategoricalRand inverts through a guide table built with the probabilities: below two expected comparisons per draw for any K, with values identical to sequential search.
- Integer negative binomial (Pascal) sampling sums geometric variates only for r <= 4. Larger shapes, which previously took r geometric draws whenever the mean exceeded 10, use the gamma–Poisson mixture with a PTRS Poisson variate; r = 1000 takes about 85 ns instead of 7–17 µs.
- HyperGeometricRand reduces draws and successes to at most N / 2 by symmetry, keeps the urn for fewer than 10 draws and otherwise uses Stadlober's HRUA ratio-of-uniforms method. Large populations take under 0.6 µs instead of time linear in the sample size.
- Static BinomialRand::Variate, used by the beta-binomial, negative-hypergeometric and trinomial samplers with a new probability per draw, uses sequential inversion for n min(p, q) < 10 and BTRS transformed rejection otherwise, instead of n min(p, q) + 1 geometric waiting times.
- NegativeHyperGeometricRand samples its beta-binomial representation unless the urn needs fewer than 16 expected draws, removing the cost linear in the population size.
The new rejection samplers evaluate their acceptance tests through Stirling's series, using the new math::logFactorialRatio, instead of differences of log-factorials near 10^10. Those differences kept only about six significant digits in double precision, and fewer where long double is binary64.
Kept the polar Cauchy sampler, which measured faster than a ratio of two normal variates, and the exact Irwin–Hall sum, whose sublinear exact alternatives are substantially more complex.

Examples and packaging

Added warning-clean, CI-compiled examples for standard-engine sampling, parallel streams, fitted model selection, and conjugate Bayesian updates. The stable-fitting example uses the stable estimator.
Added a locally validated vcpkg overlay port and Conan 2 header-only recipe with a compiled test_package consumer.
Added deterministic single-header generation and a CI smoke test. Tagged v2 releases extract their matching changelog section and attach randlib_single.hpp with its SHA-256 checksum.
Made release-candidate tags reuse the corresponding final-version changelog section and mark every SemVer prerelease GitHub release as a prerelease.
Removed a Bernoulli/binomial include cycle exposed by standalone public-header consumption without changing the exact sampling algorithm or engine advancement.

Correctness

Corrected Rayleigh UMVU scale fitting at the small-sample and asymptotic branch boundaries by consistently evaluating the exact gamma-ratio correction in the log domain (B40).
Corrected asymmetric-Laplace differential entropy to 1 + log(scale * (asymmetry + 1 / asymmetry)) (B41).
Supported zero noncentrality in noncentral chi-squared distributions, using exact central chi-squared formulas for density, tails, sampling and the mode (B42).
Defined Bernoulli entropy at deterministic endpoints by its zero-log-zero limit (B43).
Made binomial construction, scalar and bulk sampling, and modes exact at probability endpoints; near-one probabilities are no longer rounded to deterministic bulk samples (B44).
Stabilized discrete-uniform characteristic functions at integer-lattice periods with a reduced-angle sinc form (B45).
Corrected Skellam modes by searching adjacent probabilities from the median (B46).
Unified static and constructed Poisson sampling on the same exact inversion/rejection algorithms, removing the static large-rate normal approximation (B47).
Accepted the deterministic p = 1 endpoint for geometric and negative-binomial laws throughout their probability, moment and sampling APIs, while rejecting the improper static p = 0 case (B48).
Reported geometric entropy in nats consistently with the rest of the API, including its finite zero-entropy limit at p = 1 (B49).
Corrected bivariate covariance conversion from correlation, including the trinomial correlation formula and deterministic-marginal behavior (B50).
Corrected normal-inverse-gamma marginal scale, conditional sampling scale, joint CDF semantics, and the domain on which correlation is defined (B51).
Corrected negative-correlation bivariate-normal CDFs for nonzero second location or nonunit scale, including the exact countermonotonic limit (B52).
Made trinomial probability, CDF, sampling and mode evaluation exact at zero-probability and deterministic category boundaries (B53).
Hardened hypergeometric CDF traversal, sampling counters, large-population means and modes, and standardized moments for deterministic and two- or three-element populations (B54).
Prevented overflow in negative-hypergeometric moments for large integer populations (B55).

[2.0.0-beta.1] - 2026-10-17

Fitting

Began the unified non-mutating fitting API with FitResult, FitOptions, likelihood-based AIC/BIC, observed-information covariance and standard errors, and fixed-parameter MLE. Existing estimators are exposed for normal, log-normal, exponential, gamma, beta, arcsine, beta-prime, Pareto, Laplace, Lévy, Rayleigh, uniform, degenerate, Poisson, binomial, negative-binomial and discrete-uniform distributions. The generic constrained optimizer supplies full Weibull and logistic MLEs, plus MLEs for Student t, Cauchy, Gumbel, Fréchet, inverse-gamma, inverse-Gaussian and noncentral chi-squared laws. Bernoulli and geometric fits are closed form; chi-squared, Erlang and Fisher F use bounded integer likelihood searches.
Added named method-of-moments fitting, conjugate Bayesian posterior updates for the existing normal, log-normal, exponential, Poisson, binomial, negative-binomial, Weibull, Pareto and uniform models, and AIC/BIC model selection across candidate distributions.
Added continuous goodness-of-fit results for Kolmogorov–Smirnov, Anderson–Darling and Cramer–von Mises statistics with p-values, plus a discrete chi-squared test that combines sparse bins and accounts for fitted parameters.
Added seeded parameter-recovery grids over multiple sample sizes and committed fixed datasets with independently generated SciPy MLE results for nine continuous distributions. Reference-tool versions are pinned in tools/requirements.txt.
Added generalized Pareto and generalized extreme-value distributions, with high-precision mpmath reference data, generic MLEs, and stable direct survival-function implementations.
Added Hill and Pickands tail-index estimators, plot-ready Hill diagnostics, and peaks-over-threshold generalized-Pareto fitting.
Added maximum-likelihood fitting for exact, left-, right- and interval-censored observations and for one- or two-sided truncated samples. These adapters share the constrained optimizer, information criteria and observed-information covariance used by ordinary MLE.
Censored and truncated likelihood fits use a simplex search from the original initializer, avoiding gradient-based false convergence on flat transformed-likelihood regions.

API compatibility

Marked the legacy mutating member Fit* methods deprecated in favor of fit, fit_moments and fit_bayes. They remain source-compatible until v3; projects temporarily unable to migrate can define RANDLIB_DISABLE_DEPRECATED_FIT_WARNINGS before including RandLib.
Removed the unused private POISSON binomial generator selector. The exact Bernoulli-sum, waiting-time and rejection samplers remain; no approximation is introduced.

[2.0.0-alpha.1] - 2026-10-17

Architecture

Converted RandLib to a header-only C++20 library under include/randlib/ and moved the public API into namespace randlib.
Made every sampler accept any standard uniform random bit generator. Added PCG64, Xoshiro256++, SplitMix64 and deterministic stream construction while retaining the corrected legacy engines.
Removed hidden mutable sampling state. Bulk sampling uses std::span, support categories are a scoped enum, parameter factories are named, and distribution classes expose const-correct, [[nodiscard]] APIs.

Build and tooling

Added a pinned Clang 18 formatting policy and pre-commit hook, and applied it to the public headers and test sources without reordering dependency-sensitive includes.
Clang-tidy now analyzes the explicit-instantiation translation unit and public headers as a blocking CI job, with every diagnostic promoted to an error.
Added opt-in Google Benchmark sampling comparisons for normal, exponential, gamma and Poisson distributions against the standard library, plus stable-law sampling.
Engine-adapter interoperability tests now verify identical engine advancement and allow up to eight ULPs of platform-dependent variation between direct and type-erased sampling paths.
Negative-binomial maximum-likelihood fitting now optimizes its positive shape parameter in log space. Newton iterations reject non-finite states after bounded backtracking instead of stalling on a platform-dependent special-function pole (B39).

API cleanup

Kept the type-erased EngineRef sampling adapters private and removed the access-specifier churn left by the automated header conversion.
Moved implementation-only numerical helpers and tolerances into randlib::math::detail. The public Pair and Triplet aliases now live with the other core declarations, while the unconstrained arithmetic overloads for their underlying std::pair types have been removed.

Numerical kernel

Templated the elementary, gamma, beta, Bessel, zeta, Lambert W and Marcum functions on their floating-point type. Existing mixed integer/double calls remain source-compatible, while long double now stays in a separate typed path covered by closed-form and recurrence tests.
Templated factorial and combinatorial return types; calls without an explicit type retain their previous long double result.

[1.2.1] - Unreleased

Fitting and numerical methods

B39 Negative-binomial maximum-likelihood fitting optimizes the positive shape parameter in log space, preventing Newton steps from crossing a special-function pole. The Newton helper now rejects non-finite states after bounded backtracking instead of performing up to 100,000 futile iterations. A three-observation regression uses the 50-digit mpmath score root 5.2122598849254577174.
B40 Rayleigh UMVU scale fitting consistently evaluates the gamma-ratio correction in the log domain. The small-sample branch previously passed a logarithm as the scale, while the large-sample branch exponentiated half the scale squared and treated a log-ratio expansion as an ordinary correction factor.
B41 Asymmetric-Laplace entropy now returns 1 + log(scale * (asymmetry + 1 / asymmetry)); the previous log1p expression omitted the unit entropy contribution and put the scale term inside an unrelated sum.
B42 Noncentral chi-squared distributions now accept zero noncentrality, as required by their default constructor and central chi-squared limit. Density, direct tails, sampling and the mode use exact gamma-limit formulas instead of singular Bessel/Marcum expressions.
B43 Bernoulli entropy is now linkable from downstream translation units and returns zero at the accepted deterministic endpoints instead of evaluating the indeterminate floating-point products 0 * log(0).
B44 Binomial construction and cached sampling now handle the accepted deterministic probabilities without configuring an invalid geometric helper. Binomial and Bernoulli modes at p = 1 are clamped to the upper support endpoint; the ordinary floor((n + 1) * p) expression otherwise returned n + 1. Bulk sampling no longer rounds probabilities within 1e-6 of one to a deterministic law.
B45 Discrete-uniform characteristic functions now use a reduced-angle sinc form and return exactly one at integer-lattice periods instead of dividing zero by zero.
B46 Skellam modes now search adjacent probabilities from the median instead of truncating the mean, which need not maximize the PMF. A three-point regression uses 60-digit Bessel references.
B47 Static Poisson sampling now shares the exact inversion/rejection algorithms used by constructed distributions. This replaces its separate exponential-arrival algorithm through rate 1000 and, more importantly, the inexact normal approximation above that boundary.
B48 Geometric and negative-binomial distributions now accept the mathematically valid p = 1 endpoint as a deterministic zero law across probability functions, moments and sampling. The improper p = 0 geometric static sampler is rejected instead of returning the integer limit.
B49 Geometric entropy is now reported in nats, consistently with the rest of the API, and uses its finite zero-entropy limit at p = 1.
B50 Bivariate covariance matrices now convert correlation with marginal standard deviations rather than variances. Trinomial correlation uses the corresponding square-root formula and is undefined when either marginal is deterministic.
B51 Normal-inverse-gamma marginals and sampling now use the standard conditional variance sigma_squared / lambda. Its F member now evaluates the joint CDF instead of a density-like expression, and correlation is undefined when its required variances do not exist.
B52 Bivariate-normal CDFs with negative correlation now standardize the reflected second marginal before applying the normal CDF. This also restores the exact countermonotonic limit at correlation -1 when the second marginal has a nonzero location or nonunit scale.
B53 Trinomial probability, CDF and sampling methods now handle zero-probability and deterministic categories without 0 * log(0) or 0 / 0. Modes are found by a simplex-aware neighbor search and can no longer return an impossible pair whose counts exceed n.
B54 Hypergeometric means and modes no longer multiply in the integer parameter type, and standardized moments now handle deterministic and populations of size two or three. CDF summation starts at the actual support boundary, and sampling uses the distribution's integer type.
B55 Negative-hypergeometric means and variances now promote integer parameters before adding or multiplying them, preventing overflow for large populations.

[1.2.0] - 2026-10-17

A stabilization release: no API is removed and the file layout is unchanged, but sampled values change for the same seed (see B1–B6, B12, B26 below), and many probability functions now return different, correct values. Code that stored seeds to reproduce earlier runs will get different numbers.

The new test suite compares the library with 40-digit mpmath reference values, with the values recorded from v1.1, and with its own internal coherence and sampling statistics. Every fix below has a regression test in tests/unit/test_regressions.cpp under the same ID.

Build and portability

CMake build with RandLib::RandLib target, install rules and find_package(RandLib) support; presets dev, release, asan, ubsan. RandLib.pro (qmake) is deprecated and will be removed in 2.0.
Builds without warnings on GCC 11–14, Clang 15–18, AppleClang and MSVC 2022 (#18); MSVC's narrowing-conversion warnings, which GCC and Clang do not emit, are disabled. Duplicate explicit instantiations that MSVC 14.44 diagnoses as C4661 are removed. Reference tests document the double-only MSVC incomplete-beta tolerance separately. Ziggurat tables and a few constants no longer use constexpr initializers that call non-constexpr math functions. Test constants also avoid names that can trigger MSVC's template-instantiation shadow warnings.
Discrete sampling checks are discovered as one bounded CTest case per distribution. This keeps a platform-specific slow or nonterminating sampler from hiding inside a single six-hour test job.
The library no longer requires the optional C++17 special math functions, which libc++ does not provide: Bessel I and K, beta and Riemann zeta have portable implementations, used automatically when the standard library lacks them (RANDLIB_PORTABLE_SPECIAL_FUNCTIONS=ON forces them).
Constants moved to randlib::consts; the M_* names remain available when the C library does not define them. Static builds define RANDLIB_STATIC so that Windows consumers do not see dllimport.
RANDLIB_GENERATOR (CMake) or RANDLIB_JLKISS64 / RANDLIB_PCG32 (preprocessor) select the engine behind RandGenerator; JLKISS64RAND is still accepted.
Repository hygiene: removed RandLib.pro.user, LaTeX build files and the dead Travis configuration; NoncentralTRand moved to attic/ (it never compiled).

Random number generation

B1 Normal samples no longer depend on other instances: the ziggurat kept a leftover exponential in a static thread_local variable.
B2 JKiss and JLKiss64 seeding expands the seed with SplitMix64 into the whole state; seeds that differed only above bit 32 gave identical streams, and the streams of nearby seeds were strongly correlated.
B3 PCG seeding follows the reference pcg32_srandom_r; PCGRandEngine::Reseed(seed, sequence) selects the stream.
B4 Default seeds combine std::random_device, a high-resolution clock, the thread id and a per-process counter.
B5 The 64-bit engine switch works: one uniform generator derives its bit count from the engine.
B6 Uniform variates use the full mantissa (52 random bits for double, 63 for x87 long double) instead of 32 bits, so inverse-transform tails reach the full range; float uniforms can no longer round to 1.
B9 The normal sign comes from an explicit bit.
B12 BernoulliRand::StandardVariate no longer caches bits in static variables shared by all generators and threads.
B26 Composite distributions reseed every generator they use (Laplace and geometric stable families, Planck), and member generators get seeds derived with RandEngine::DeriveSeed instead of seed + k, which collided with the seed of a neighbouring object.

Correctness

B7 Stable laws with exponent within 0.01 of 1 and non-zero skewness can be constructed and sampled. Their density and distribution function throw the new randlib::not_implemented until a dedicated algorithm exists (#1).
B8 The distribution function of totally skewed stable laws no longer has a floor of 3.8·10⁻⁵ (the adaptive Simpson integrator missed narrow peaks).
B11 Generic quantiles are deterministic and no longer consume the instance's random stream: the initial guess came from sampling into a shared static vector (a data race between threads). Quantiles are now found by bracketing and the Illinois method to full double precision, in the tail with the smaller probability; the old Newton iteration stopped at tolerances of 1e-6–1e-10.
B13 M_LNPI (log π) was wrong in the 11th digit, which affected every density that uses it.
B14 NormalRand::ThirdMoment and FourthMoment returned wrong values.
B15 Classes that could not be used: BaldingNicholsRand::F was hidden by a data member, CenteredAsymmetricLaplaceRand was abstract, WrappedExponentialRand and ShiftedAsymmetricLaplaceDistribution were never instantiated, and the abstract bases had protected destructors, so distributions could not be owned through a base pointer.
B16 WrappedExponentialRand: circular variance, characteristic function at non-integer t and rate validation. VonMisesRand characteristic function at non-integer t.
B17 Regularized incomplete gamma functions are accurate to about 1e-13 everywhere; they were off by up to 1e-3 for large shapes (Gamma, chi-squared, Poisson, Erlang).
B18 logBeta no longer underflows to −∞ for large arguments (e.g. Beta(600, 600)).
B19 Triangular distributions accept a mode equal to a limit.
B20 Modified Bessel function of negative and half-integer order: the sinh/cosh cases were swapped and the reflection had the wrong sign, which broke the noncentral chi-squared density for degrees below 2 and equal to 3.
B21 Stable laws: the characteristic function had the skewness term with the wrong sign for α ≠ 1, and for α = 1 the density, distribution function and sampler used the opposite sign of the S1 location shift (2/π)βγ log γ.
B22 Stable densities and distribution functions in the tails: the integrals are split around the peak of the integrand, the heavy tail of totally skewed laws with α > 1 no longer returns 0, and survival functions are computed directly instead of as 1 − F.
B23 Survival functions and tails computed without cancellation or overflow: Cauchy, Student t (also near the median), beta prime, Fisher F, inverse Gaussian, logistic, Wigner semicircle.
B24 Discrete distributions: logarithmic distribution function, binomial median, categorical quantiles, negative binomial survival below the support, Bernoulli log-probability outside the support.
B25 Invalid parameters throw std::invalid_argument instead of being replaced silently (log-normal variance, Balding-Nichols, negative binomial), and NaN parameters and non-finite locations are rejected throughout.
B27 Characteristic functions of heavy-tailed discrete laws (Zeta, Yule) terminate; their summation used to need about 10¹⁰ terms. The tail beyond the summed terms is estimated with Euler's transformation, so a value takes milliseconds and is accurate to about 1e-16.
B28 The Cantor characteristic function ignored t.
B29 The geometric stable density and distribution function for α = 1 missed the S1 location shift of the mixing representation.
B30 The Marchenko–Pastur distribution function next to the lower edge of the support.
B31 RademacherRand::Mode returned a random value; it is now −1, like the median.
B32 Poisson and binomial probabilities use Loader's saddle-point form: the sum of logarithms lost up to four digits for large arguments, which depended on the platform's lgamma.
B33 The raised cosine (and Raab–Green) characteristic function at st = π evaluated 0 / 0; it is now computed as π² sinc(π − st) / (st (π + st)).
B34 Discrete quantiles at a probability level that F reaches exactly no longer depend on the last bits of the platform's special functions.
B35 Zeta distribution and survival functions use an Euler–Maclaurin tail for large arguments instead of summing billions of terms.
B36 Continuous and discrete uniform likelihoods negate the sample-count exponent after converting it to a signed floating-point value; negating size_t previously wrapped and produced an enormous positive exponent.
B37 Double factorial now evaluates the defining product n(n − 2)…; the previous factorial identities shifted by the full input and returned incorrect values for every n > 1.
B38 Geometric and negative-binomial sampling saturate results that exceed the result integer type instead of performing undefined conversions or additions. Binomial waiting-time sampling compares against the remaining trials before addition, preventing overflow and a beta-binomial hang.
Removed the unused private POISSON binomial generator selector. Binomial sampling continues to use its exact Bernoulli-sum, waiting-time and rejection algorithms rather than an approximation.
Noncentral chi-squared mode for 2 degrees of freedom; DegenerateRand::Fit accepted any sample.
Numerical integration uses adaptive Gauss–Kronrod (21-point) quadrature with a global error estimate instead of adaptive Simpson (#10).
Quantile searches of stable and geometric stable laws start the far tails from their power-law asymptote.

Known limitations

Tracked by the tolerance overrides in tests/data/reference_tolerances.csv:

characteristic functions of heavy-tailed laws computed by numerical integration are accurate to about 1e-3 (#10), the log-normal one to about 1e-4 for large variance (#6);
noncentral chi-squared upper tails come from 1 − Marcum P (#4);
stable densities near α = 1 with β ≠ 0 are not implemented (#1);
geometric stable densities and distribution functions take about 0.1 s per evaluation (up to a few seconds in the far tails for α = 1);
high-precision reference values for the skewed stable and geometric stable laws with α = 1 are still being computed; these cases are covered by spot checks against mpmath in the regression tests.

[1.1.0]

Last release before the modernization; tagged in the repository as v1.1.0-baseline (6536dad).