Module special 

Scalar special-function primitives shared across the workspace.

These are pure (std/libm-only) numeric kernels with no upward crate dependencies, so they live in the lowest crate (gam-math) and can be consumed by any term/basis/inference code without inducing an SCC edge.

Functions

bessel_i0_centered_second_log_derivative_from_log_abs

Second log-scale derivative of the centered Bessel primitive: d²/d(log η)²[log I0(η) − η], i.e. the derivative of the third term d1 returned by bessel_i0_centered_terms (d1 = η d/dη[log I0(η) − η]).

bessel_i0_centered_terms

Overflow-free centered Bessel value, ratio, and log-scale derivative.

bessel_i0_centered_terms_from_log_abs

Stable centered Bessel terms when only log(|eta|) is representable.

bessel_i0_log_and_ratio

Overflow-free (log I0(eta), I1(|eta|) / I0(|eta|)).

bessel_i0_log_minus_abs_and_ratio

Overflow-free (log I0(eta) - |eta|, I1(|eta|) / I0(|eta|)).

binomial_coefficient_f64

Numerically stable C(n,k) = n! / (k!·(n−k)!) as f64. Uses the symmetry C(n,k) = C(n, n−k) to keep the loop count min(k, n−k) and the multiplicative recurrence C(n,j+1) = C(n,j)·(n−j)/(j+1), avoiding the overflow of separate factorial evaluations. Returns 0.0 for k > n and exact integer results within 2^53.

digamma

Digamma ψ(x) = d/dx ln Γ(x), for x > 0; NaN otherwise.

gauss_legendre

Gauss-Legendre nodes and weights on [-1, 1] for n points, computed via Newton iteration on the Legendre-polynomial roots (Bonnet’s three-term recurrence, cosine initial guess). Returns (nodes, weights) with nodes ascending; for odd n the central node is exactly 0.0.

pentagamma

Pentagamma ψ₃(x) = d⁴/dx⁴ ln Γ(x), for x > 0; NaN otherwise.

stable_polynomial_times_exp_neg

Evaluate (Σ_k coeffs[k]·x^k) · exp(−x) without overflow. For moderate x ≤ 600 uses Horner + exp(−x) directly; for very large x rewrites xᵈ · exp(−x) = exp(d·ln x − x) and runs Horner in 1/x, which keeps both the polynomial sum and its multiplier inside double range. Returns 0.0 for non-finite x or empty coeffs.

tetragamma

Tetragamma ψ₂(x) = d³/dx³ ln Γ(x), for x > 0; NaN otherwise.

trigamma

Trigamma ψ₁(x) = d²/dx² ln Γ(x), for x > 0; NaN otherwise.