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

bd0

The Bregman divergence bd0(x, m) = x·ln(x/m) + m − x of the Poisson deviance (Loader’s bd0), summed as a series in (x − m)/(x + m) when the two arguments are within 20% of each other so the two ~equal logarithms never cancel. bd0(0, m) = m exactly.

bernoulli_kl_from_logits

Bernoulli KL divergence in natural coordinates, KL(σ(a) ‖ σ(b)), without subtracting an entropy from a cross entropy. For |b − a| ≤ 1/2 only second-order remainders are evaluated; the tail branches orient the event so the reference probability never rounds to one.

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.

expm1_minus_x

e^x − 1 − x, the second-order remainder of the exponential. For |x| ≤ 1/2 the Taylor tail is summed directly so the result does not cancel against x; beyond that exp_m1 is accurate on its own.

exprel

The relative exponential exprel(x) = (e^x − 1) / x, equal to 1 at x = 0 and summed as a series for |x| ≤ 1/2.

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.

gauss_lobatto

Gauss-Lobatto nodes and weights on [-1, 1], ascending, for n >= 2.

log1p_minus_x

ln(1 + x) − x, the second-order remainder of the logarithm, with the same small-argument series treatment as expm1_minus_x.

log_abs_one_minus_exp

ln|1 − e^x| for x ≠ 0, routed through crate::probability::log1mexp_positive on both sides of the origin.

log_exprel

ln(exprel(x)), with the large-|x| branches written so that no exponential of a positive argument is ever formed.

logaddexp

ln(e^a + e^b) without forming either exponential at full scale; returns -∞ when both arguments are -∞.

logistic

The logistic function σ(x) = 1 / (1 + e^{-x}), oriented so the exponential is always of a non-positive argument.

pentagamma

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

positive_frexp

Exact power-of-two decomposition x = mantissa · 2^exponent for a positive finite f64, including subnormals. The mantissa lies in [1, 2).

scale_normalized_power_of_two

mantissa · 2^exponent for a positive mantissa, renormalising the mantissa into [1, 2) first so the result overflows or underflows only when the final f64 itself is unrepresentable. Subnormal results are formed by scaling in units of the least positive subnormal, so IEEE rounds the final value once instead of underflowing an intermediate.

scaled_positive_product_quotient

a·b·c/d for positive finite inputs, carrying the binary exponent separately so an intermediate overflow or underflow cannot change a representable final result.

softplus

softplus(x) = ln(1 + e^x), evaluated as max(x, 0) + ln(1 + e^{-|x|}) so that neither tail overflows and the small-x result keeps full relative accuracy.

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.

xlogy

x·ln(y) with the convention 0·ln(0) = 0 used by every deviance kernel.