Module special 

Special functions underpinning the probability distributions — the canonical, full-double-precision implementations the whole engine shares.

These replace the scattered ad-hoc approximations (a coarse normal_cdf copied into financial_modeling, a |t| > 1.96 ⇒ p = 0.05 placeholder in the t-test): every distribution CDF/quantile is built from ln_gamma, the regularized incomplete gamma P(a,x), and the regularized incomplete beta I_x(a,b) here.

Algorithms are the standard, well-conditioned ones (Lanczos for ln_gamma; a series + continued-fraction split for the incomplete gamma; a Lentz continued fraction for the incomplete beta) — they are ordinary numerical mathematics, accurate to ~1e-12 and verified against known closed forms in the tests.

All scalar f64 — these are pointwise special functions, not GPU-amenable (per-call scalar work below any dispatch crossover); the CPU path is the right and only path here (CLAUDE.md §13: “where GPU does not help, say so”).

Functions

betai

Regularized incomplete beta I_x(a, b), 0 ≤ x ≤ 1, a,b > 0. Continued fraction (Lentz) with the symmetry switch for fast convergence.

erf

Error function erf(x) via the incomplete gamma: erf(x) = sign(x)·P(½, x²).

erfc

Complementary error function erfc(x) = 1 − erf(x).

gamma

Γ(x) for convenience (small/moderate x).

gammp

Regularized lower incomplete gamma P(a, x) = γ(a,x)/Γ(a), a > 0, x ≥ 0. Series for x < a+1, continued fraction (via Q) otherwise.

gammq

Regularized upper incomplete gamma Q(a, x) = 1 − P(a, x).

ln_gamma

Natural log of the Gamma function, ln Γ(x), via the Lanczos approximation (g = 7, 9 coefficients) with the reflection formula for x < 0.5. Accurate to ~15 significant figures for x > 0.