Trait ScalarSpecialMath 

pub trait ScalarSpecialMath: ScalarSpecialMathWithPolicy {
    // Provided methods
    fn scalar_erf(self) -> Self { ... }
    fn scalar_erfc(self) -> Self { ... }
    fn scalar_erfcx(self) -> Self { ... }
    fn scalar_logistic_sigmoid(self) -> Self { ... }
    fn scalar_logit(self) -> Self { ... }
    fn scalar_logit_1m(self) -> Self { ... }
    fn scalar_softplus(self, k: Self, rcp_k: Self) -> Self { ... }
    fn scalar_tgamma(self) -> Self { ... }
    fn scalar_lgamma(self) -> Self { ... }
    fn scalar_poisson_pmf(self, lambda: Self) -> Self { ... }
    fn scalar_poisson_log_pmf(self, lambda: Self) -> Self { ... }
    fn scalar_digamma(self) -> Self { ... }
    fn scalar_trigamma(self) -> Self { ... }
    fn scalar_polygamma(self, n: u32) -> Self { ... }
    fn scalar_zeta(self) -> Self { ... }
    fn scalar_zetac(self) -> Self { ... }
    fn scalar_polylog(self, order: PolylogOrder<Self, Self::Signed>) -> Self { ... }
    fn scalar_bessel_n<F: BesselFamily, const N: i32>(self) -> Self { ... }
    fn scalar_bessel<F: BesselFamily>(
        self,
        order: BesselOrder<Self, Self::Signed>,
    ) -> Self { ... }
    fn scalar_sph_bessel_n<F: BesselFamily, const N: usize>(self) -> Self { ... }
    fn scalar_sph_bessel<F: BesselFamily>(self, n: u32) -> Self { ... }
    fn scalar_airy<W: AiryFn>(self) -> Self { ... }
    fn scalar_airy_all<const SCALED: bool>(self) -> (Self, Self, Self, Self) { ... }
    fn scalar_beta(self, y: Self) -> Self { ... }
    fn scalar_lbeta(self, y: Self) -> Self { ... }
    fn scalar_jacobi(self, alpha: Self, beta: Self, n: u32, m: u32) -> Self { ... }
    fn scalar_hermite_n<const N: usize>(self) -> Self { ... }
    fn scalar_hermitev(self, n: Self::Unsigned) -> Self { ... }
    fn scalar_hermite(self, n: u32) -> Self { ... }
    fn scalar_hermite_function_n<const N: usize>(self) -> Self { ... }
    fn scalar_hermite_function(self, n: u32) -> Self { ... }
    fn scalar_hermite_function_series_n<const N: usize>(
        self,
        coeffs: &[Self::Element; N],
    ) -> Self { ... }
    fn scalar_hermite_function_series(self, coeffs: &[Self::Element]) -> Self { ... }
    fn scalar_laguerre_n<const N: usize>(self, alpha: Self) -> Self { ... }
    fn scalar_laguerrev(self, alpha: Self, n: Self::Unsigned) -> Self { ... }
    fn scalar_laguerre(self, alpha: Self, n: u32) -> Self { ... }
    fn scalar_laguerre_function_n<const N: usize>(self, alpha: Self) -> Self { ... }
    fn scalar_laguerre_function(self, alpha: Self, n: u32) -> Self { ... }
    fn scalar_laguerre_function_i_n<const N: usize>(self, alpha: i32) -> Self { ... }
    fn scalar_laguerre_function_i(self, alpha: i32, n: u32) -> Self { ... }
    fn scalar_laguerre_function_series_n<const N: usize>(
        self,
        alpha: Self,
        coeffs: &[Self::Element; N],
    ) -> Self { ... }
    fn scalar_laguerre_function_series(
        self,
        alpha: Self,
        coeffs: &[Self::Element],
    ) -> Self { ... }
    fn scalar_laguerre_function_series_i_n<const N: usize>(
        self,
        alpha: i32,
        coeffs: &[Self::Element; N],
    ) -> Self { ... }
    fn scalar_laguerre_function_series_i(
        self,
        alpha: i32,
        coeffs: &[Self::Element],
    ) -> Self { ... }
    fn scalar_chebyshev_n<const K: usize, const N: usize>(
        self,
        coeffs: &[Self::Element; N],
    ) -> Self { ... }
    fn scalar_chebyshev<const K: usize>(self, coeffs: &[Self::Element]) -> Self { ... }
    fn scalar_gaussian(self, a: Self, c: Self) -> Self { ... }
    fn scalar_planck(self) -> Self { ... }
    fn scalar_legendre(self, n: u32, m: u32) -> Self { ... }
    fn scalar_legendre_series_n<const N: usize>(
        self,
        coeffs: &[Self::Element; N],
    ) -> Self { ... }
    fn scalar_legendre_series(self, coeffs: &[Self::Element]) -> Self { ... }
    fn scalar_zernike_r(self, n: u32, m: u32) -> Self { ... }
    fn scalar_zernike<const NORM: u8>(self, theta: Self, n: u32, m: i32) -> Self { ... }
    fn scalar_zernike_basis<const L: usize, const NORM: u8, const N: usize>(
        x: Self,
        y: Self,
        out: &mut [Self; N],
    ) { ... }
    fn scalar_lambert_w(self) -> (Self, Self) { ... }
    fn scalar_expint_n<const N: usize>(self) -> Self { ... }
    fn scalar_expint(self, n: u32) -> Self { ... }
    fn scalar_phi_n<const N: usize>(self) -> Self { ... }
    fn scalar_phi(self, n: u32) -> Self { ... }
    fn scalar_carlson<K: WrapTo>(kind: K) -> Self
       where K::Wrapped: CarlsonKind,
             <K::Wrapped as CarlsonKind>::Output: Unwrap<Unwrapped = Self> { ... }
    fn scalar_ellint<K: WrapTo>(kind: K) -> Self
       where K::Wrapped: EllipticKind,
             <K::Wrapped as EllipticKind>::Output: Unwrap<Unwrapped = Self> { ... }
    fn scalar_erfinv(self) -> Self { ... }
    fn scalar_probit(self) -> Self { ... }
    fn scalar_ndtr(self) -> Self { ... }
    fn scalar_log_ndtr(self) -> Self { ... }
    fn scalar_logerfc(self) -> Self { ... }
    fn scalar_fresnel(self) -> (Self, Self) { ... }
    fn scalar_fresnel_c(self) -> Self { ... }
    fn scalar_fresnel_s(self) -> Self { ... }
    fn scalar_sici(self) -> (Self, Self) { ... }
    fn scalar_sinint(self) -> Self { ... }
    fn scalar_cosint(self) -> Self { ... }
    fn scalar_inv_log_ndtr(self) -> Self { ... }
    fn scalar_inv_digamma(self) -> Self { ... }
    fn scalar_wright_omega(self) -> Self { ... }
    fn scalar_bessel_ratio<F: BesselRatioFamily>(self, nu: Self) -> Self { ... }
    fn scalar_inv_bessel_ratio<F: BesselRatioFamily>(self, nu: Self) -> Self { ... }
    fn scalar_bessel_ratio_1m<F: BesselRatioFamily>(self, nu: Self) -> Self { ... }
    fn scalar_inv_bessel_ratio_1m<F: BesselRatioFamily>(self, nu: Self) -> Self { ... }
    fn scalar_gauss_legendre(self, n: u32) -> (Self, Self) { ... }
    fn scalar_gauss_hermite(self, n: u32) -> (Self, Self) { ... }
    fn scalar_gauss_laguerre(self, alpha: Self, n: u32) -> (Self, Self) { ... }
    fn scalar_pochhammer(self, m: Self) -> Self { ... }
    fn scalar_jacobi_elliptic(self, k: Self) -> (Self, Self, Self) { ... }
    fn scalar_agm(self, other: Self) -> Self { ... }
    fn scalar_langevin(self) -> Self { ... }
    fn scalar_inv_langevin(self) -> Self { ... }
    fn scalar_langevin_1m(self) -> Self { ... }
    fn scalar_inv_langevin_1m(self) -> Self { ... }
    fn scalar_gelu(self, alpha: Self) -> Self { ... }
    fn scalar_swish(self, beta: Self) -> Self { ... }
    fn scalar_algebraic_sigmoid_n<const N: usize>(self) -> Self { ... }
    fn scalar_algebraic_sigmoid(self, n: u32) -> Self { ... }
    fn scalar_algebraic_swish(self) -> Self { ... }
    fn scalar_lgamma_r(self) -> (Self, Self) { ... }
    fn scalar_gaussian_integral(x0: Self, x1: Self, a: Self, c: Self) -> Self { ... }
    fn scalar_boxcox(self, lambda: Self) -> Self { ... }
    fn scalar_boxcox_1p(self, lambda: Self) -> Self { ... }
    fn scalar_inv_boxcox(self, lambda: Self) -> Self { ... }
    fn scalar_inv_boxcox_1p(self, lambda: Self) -> Self { ... }
    fn scalar_yeo_johnson(self, lambda: Self) -> Self { ... }
    fn scalar_inv_yeo_johnson(self, lambda: Self) -> Self { ... }
    fn scalar_spherical_harmonics<const L: usize, const N: usize, const CS: bool>(
        x: Self,
        y: Self,
        z: Self,
        out: &mut [Self; N],
    ) { ... }
    fn scalar_spherical_harmonics_table<const L: usize, const N: usize, const CS: bool>(
        table: &mut ShTable<Self, N>,
    ) { ... }
    fn scalar_spherical_harmonics_with<const L: usize, const N: usize>(
        table: &ShTable<Self, N>,
        x: Self,
        y: Self,
        z: Self,
        out: &mut [Self; N],
    ) { ... }
    fn scalar_spherical_harmonics_d<const L: usize, const N: usize, const CS: bool>(
        x: Self,
        y: Self,
        z: Self,
        out: &mut [Self; N],
        ddx: &mut [Self; N],
        ddy: &mut [Self; N],
        ddz: &mut [Self; N],
    ) { ... }
    fn scalar_zernike_basis_d<const L: usize, const NORM: u8, const N: usize>(
        x: Self,
        y: Self,
        out: &mut [Self; N],
        ddx: &mut [Self; N],
        ddy: &mut [Self; N],
    ) { ... }
    fn scalar_spherical_harmonics_d_with<const L: usize, const N: usize>(
        table: &ShTable<Self, N>,
        x: Self,
        y: Self,
        z: Self,
        out: &mut [Self; N],
        ddx: &mut [Self; N],
        ddy: &mut [Self; N],
        ddz: &mut [Self; N],
    ) { ... }
    fn scalar_softplus_d(self, k: Self, rcp_k: Self) -> (Self, Self) { ... }
    fn scalar_gelu_d(self, alpha: Self) -> (Self, Self) { ... }
    fn scalar_swish_d(self, beta: Self) -> (Self, Self) { ... }
    fn scalar_algebraic_sigmoid_d_n<const N: usize>(self) -> (Self, Self) { ... }
    fn scalar_algebraic_sigmoid_d(self, n: u32) -> (Self, Self) { ... }
    fn scalar_algebraic_swish_d(self) -> (Self, Self) { ... }
    fn scalar_langevin_d(self) -> (Self, Self) { ... }
}

Aggregate of all scalar math traits using the default policy.

This trait collects every method from the following trait families into a single trait implemented directly on f32 and f64, using DefaultPolicy for all operations:

• SpecialMath
• RealSpecialMath
• RealPrimalMath

All methods are prefixed with scalar_. See ScalarSpecialMathWithPolicy for the policy-aware variant, which additionally carries a _p suffix on each method.

Limitations

Implemented only for bare scalar types; see ScalarSpecialMathWithPolicy. All types implementing ScalarSpecialMathWithPolicy automatically implement this trait.

Provided Methods

fn scalar_erf(self) -> Self

Computes the error function.

For f32 vectors, this is still decently accurate even with the Medium and Worst precision policies, thanks to good approximations that don’t rely on the precision of exp. Subsequently, performance of the lower precision policies is excellent. Furthermore, if using on a GPU with native exp support, all precision policies will have good performance and accuracy.

Below Best, the f64 kernel forms erf as $1 - m\,e^{-x^2}$, whose error is a fixed absolute ulp of 1: erf(0) comes out 2.2e-16 and erf(1e-8) is only 2e-8 relative. From Best up, |x| < 0.84375 takes a direct $x + x\,R(x^2)/S(x^2)$ arm that is exact at zero and relatively accurate down to the denormals. The f32 kernel carries that arm from Average.

fn scalar_erfc(self) -> Self

Computes the complementary error function.

The f64 kernel is one product of six rationals times $e^{-x^2}$ over the whole line, within about 3 ulp everywhere on hardware with a fused multiply-add: the one error that grows, the rounding of $x^2$ under the exponential amplified by $x^2$, is removed with the exact residual of the product at every tier. Without a native FMA that residual is unavailable, so Best removes the growth with a bit-split of x instead, and the lower tiers keep it (47 ulp at x = 14, 237 at x = 24).

fn scalar_erfcx(self) -> Self

Computes the scaled complementary error function, $\operatorname{erfcx}(x) = e^{x^2}\operatorname{erfc}(x)$.

erfc underflows to zero at x ~ 27 in f64 and x ~ 9 in f32, where the true value is $e^{-x^2}/(x\sqrt{\pi})$, nonzero and merely too small to represent. Anything reading a Gaussian tail past that point silently gets zero: importance weights, log-likelihoods, censored-data models, the Voigt profile. erfcx removes the exponential and decays only as $1/(x\sqrt{\pi})$, so it is representable for every finite argument and keeps full relative accuracy.

Computed on the real backends as the Faddeeva function restricted to the imaginary axis, $w(ix) = \operatorname{erfcx}(x)$, where Weideman’s rational approximation degenerates to real arithmetic: one reciprocal and one Horner, no transcendental at all for x >= 0. That makes it cheaper than the erfc it complements, and measures 1.22 ulp worst over $x \in [0, 10^{15}]$ at the Best tier and above.

Negative arguments use $\operatorname{erfcx}(-x) = 2e^{x^2} - \operatorname{erfcx}(x)$ and legitimately overflow below about -26.6 (f64), the function itself growing like $e^{x^2}$ in that direction.

The two are related by $\operatorname{erfc}(x) = e^{-x^2}\operatorname{erfcx}(x)$, which is the numerically sound way to recover a tail value that erfc alone cannot hold. Keep the $-x^2$ in the log domain rather than exponentiating it.

fn scalar_logistic_sigmoid(self) -> Self

Computes the Logistic sigmoid function, defined as $\sigma(x) = \frac{1}{1 + e^{-x}}$.

It’s worth mentioning that the derivative of the logistic sigmoid can be computed very cheaply from the output of the logistic sigmoid itself, in the form of:

let s = x.logistic_sigmoid();
let derivative = s * (1.0 - s); // or s.nmul_adde(s, s), which may be slightly faster

Notably, for f32 and f64 this implementation still has good precision for the Worst precision policy, and for the Best precision policies handles very large positive and negative inputs without overflow or underflow issues.

fn scalar_logit(self) -> Self

Computes the logit $\ln\!\frac{p}{1-p}$, the inverse of logistic_sigmoid.

Evaluated as $\ln(p) - \ln_{1p}(-p)$, which is accurate for small p where the direct quotient is not. For p approaching 1 no evaluation order helps. $1 - p$ has already lost its low digits inside the input itself, and the information is not recoverable from p. A caller who knows $q = 1 - p$ should pass it to logit_1m instead, which is exact at the far end of the range.

p = 0 gives -∞, p = 1 gives +∞, and p outside [0, 1] is out of domain.

fn scalar_logit_1m(self) -> Self

Computes $\mathrm{logit}(1 - q) = \ln\!\frac{1-q}{q}$ from the complement q directly.

The companion entry point to logit, in the same relationship as langevin_1m has to langevin. The logit diverges as its argument approaches 1, and near that end $1 - p$ cannot be formed from p without losing every digit that matters. Working in q throughout sidesteps that: evaluated as $\ln_{1p}(-q) - \ln(q)$, accurate to a few ulp however small q is.

Note the sign convention follows the substitution, so logit_1m(q) == -logit(q) as functions of the same number. The two differ in which probability the argument names.

fn scalar_softplus(self, k: Self, rcp_k: Self) -> Self

Computes the softplus function, defined as $\frac{1}{k}\ln(1 + e^{kx})$.

This is a smooth approximation to the ReLU function that is more numerically stable for large inputs.

The parameter k controls the steepness of the curve, with larger values approaching ReLU more closely. Pass k = 1 and rcp_k = 1 for the standard softplus with no steepness scaling.

rcp_k must equal 1/k. It is passed explicitly so callers that invoke softplus repeatedly with the same k can pre-compute the reciprocal once rather than recomputing it per call.

To also obtain the derivative with respect to x, use softplus_d.

fn scalar_tgamma(self) -> Self

Computes the Gamma function ($\Gamma(z)$) for any real input, for each value in a vector.

This implementation uses a few different behaviors to ensure the greatest precision where possible.

• For non-integer positive inputs, it uses the Lanczos approximation.
• For small non-integer negative inputs, it uses the recursive identity $\Gamma(z) = \Gamma(z+1)/z$ until z is positive.
• For large non-integer negative inputs, it uses the reflection formula $-\pi / (\Gamma(z)\sin(\pi z)\,z)$.
• For positive integers, it simply computes the factorial in a tight loop to ensure precision. Lookup tables could not be used with SIMD.
• At zero, the result will be positive or negative infinity based on the input sign (signed zero is a thing).

NOTE: The Gamma function is not defined for negative integers.

fn scalar_lgamma(self) -> Self

Computes the natural log of the Gamma function ($\ln|\Gamma(x)|$) for any real input, for each value in a vector.

fn scalar_poisson_pmf(self, lambda: Self) -> Self

The Poisson probability mass $P(k; \lambda) = e^{-\lambda}\lambda^k / k!$ at k = self, for real $k \ge 0$ and mean $\lambda \ge 0$.

Not exp(k ln lambda - lambda - lgamma(k+1)): that forms an $O(1)$ answer as the exponential of a difference of large numbers, and half an ulp of $\ln\Gamma(k+1) = O(k \ln k)$ becomes that many ulp of the mass. For $k \ge 9$ this uses Loader’s saddle-point form (the one R’s dpois uses),

P(k; \lambda) = \frac{e^{-\mathrm{stirlerr}(k) - \mathrm{bd0}(k, \lambda)}}{\sqrt{2\pi k}}

with stirlerr the Stirling remainder (a short $1/k^2$ series) and bd0 the deviance $k \ln(k/\lambda) + \lambda - k$ (a series in $(k-\lambda)/(k+\lambda)$ near the peak, where the direct form cancels): both are small where the mass is not negligible, so the exponential amplifies nothing, and there is no lgamma and no ln at all near the peak. Below $k = 9$ the same machinery is used after shifting k up by an integer, with the exact product $(k+1)\cdots(k+m)$ taken back out, so there is no lgamma anywhere, and mixed vectors share one ln, one stirlerr and one exp. Real k is allowed because the Gamma density is the same function: $f(x; a) = P(a-1; x)$ for shape $a \ge 1$ (unit scale).

Edges: $\lambda = 0$ gives 1 at $k = 0$ and 0 above; $k = 0$ is $e^{-\lambda}$.

fn scalar_poisson_log_pmf(self, lambda: Self) -> Self

$\ln P(k; \lambda)$, the log of poisson_pmf, formed directly (no exp then ln) so it stays finite far in the tails where the mass itself underflows.

fn scalar_digamma(self) -> Self

Computes the digamma function $\psi(x) = \frac{\mathrm{d}}{\mathrm{d}x}\ln\Gamma(x) = \frac{\Gamma'(x)}{\Gamma(x)}$ for any real input, for each value in a vector.

The argument is handled in three regimes:

• For x >= 10, an asymptotic expansion in $1/x^2$ is used.
• For smaller x, the recurrence $\psi(x) = \psi(x+1) - 1/x$ shifts the argument into [1, 2], where a rational minimax approximation $\psi(x) = (x - x_0)(Y + R(x-1))$ is used ($x_0$ is the positive root of $\psi$).
• For x <= -1, the reflection formula $\psi(1-x) = \psi(x) + \pi\cot(\pi x)$ is applied.

NOTE: The digamma function is not defined at zero or the negative integers. Those inputs yield NaN when overflow checking is enabled.

fn scalar_trigamma(self) -> Self

Computes the trigamma function $\psi_1(x) = \frac{\mathrm{d}}{\mathrm{d}x}\psi(x)$, the second derivative of $\ln\Gamma$.

Real vectors run a dedicated kernel (three minimax rational regions with a single recurrence step and the $\pi^2/\sin^2(\pi x)$ reflection) that is a little tighter than the general polygamma machinery at order 1. polygamma(1) routes here, so the two spellings agree exactly. Complex vectors have their own implementation, which is the reason this lives on SpecialMath while polygamma is real-only.

The poles at zero and the negative integers evaluate to +inf: $\psi_1$ has double poles, so unlike digamma the two one-sided limits agree.

fn scalar_polygamma(self, n: u32) -> Self

Computes the polygamma function $\psi_n(x) = \frac{\mathrm{d}^n}{\mathrm{d}x^n}\psi(x)$, the n-th derivative of digamma (n = 0 is digamma, n = 1 is trigamma).

The order n is a runtime scalar shared by every lane. That is a deliberate design choice: it closes the Gamma family under differentiation, since $\psi_n'(x) = \psi_{n+1}(x)$ is reachable by passing n + 1, which is what lets forward-mode AD (Dual) differentiate through any member of the family to any depth. All order-dependent coefficients are scalar work splatted once, so uniform n costs a vector nothing.

For n >= 2, real vectors run a masked recurrence up to the transition point $N = 0.4\,d_{10} + 4n$ and then the Bernoulli asymptotic series on the positive axis. Negative arguments reflect through the n-th derivative of $\cot(\pi x)$ (tabulated to n = 20, above which negative arguments return NaN). At zero and the negative integers, odd n returns +inf (the correct two-sided limit) and even n has one-sided limits of opposite sign, so it returns NaN when overflow checking is enabled.

Complex vectors run the same recurrence-plus-series in complex arithmetic, gated on $\operatorname{Re} z$, reflecting the half-plane $\operatorname{Re} z < 1/2$ through the same tabulated $\cot$ derivative (so the n <= 20 reflection reach applies there too). Only psi_n of a real variable is real, so this is the family member that makes polygamma complex-capable at all orders.

Orders where $n!$ overflows the element type (n >= 171 for f64, n >= 35 for f32) return the signed infinity carried by the leading term on the real positive axis, and NaN over C.

fn scalar_zeta(self) -> Self

Computes the Riemann zeta function $\zeta(s) = \sum_{n\ge1} n^{-s}$.

Evaluated as 1 + zetac, which is where the accuracy argument lives (see there). Worst relative error measured against mpmath at 40 digits: 4.4e-16 for s in [1.5, 5], 4.3e-16 for [5, 40], 2.3e-15 through the critical strip [0.1, 0.9], and 4.6e-16 approaching the pole at s = 1, which returns infinity.

Negative s goes through the functional equation $\zeta(s) = 2^s\pi^{s-1}\sin(\pi s/2)\,\Gamma(1-s)\,\zeta(1-s)$, landing back at $1-s > 1$ where the series is at its most accurate. That arm costs a tgamma and a sin_pi beyond the main path, so it is gated on a lane needing it.

This is the Riemann zeta of one real argument. The two-argument Hurwitz form $\zeta(s, q)$ is not provided: it generalizes the same expansion but loses the prime factorization that makes this one cheap, so it is a separate and materially more expensive function rather than a special case of this one.

fn scalar_zetac(self) -> Self

Computes $\zeta(s) - 1$, accurately where $\zeta(s)$ is within rounding of 1.

$\zeta$ approaches 1 quickly: $\zeta(40) - 1$ is about 9.1e-13, already below the mantissa of $\zeta$ itself, and $\zeta(80) - 1$ is 8.3e-25. Forming zeta and subtracting 1 therefore destroys the answer: at s = 40 it is off by 9e-8 relative, at s = 80 by 100%, and past s = 200 it returns a flat zero.

This is not a wrapper around that subtraction. The Euler-Maclaurin sum underneath opens with the $n = 1$ term, which is the 1, so the complement is obtained by omitting it, with no cancellation anywhere and still full relative accuracy at s = 700, where the value is around 1e-211. $\zeta$ is the derived form here, the same way exp relates to exp_m1.

Same accuracy and the same negative-s handling as zeta.

fn scalar_polylog(self, order: PolylogOrder<Self, Self::Signed>) -> Self

Computes the polylogarithm $\mathrm{Li}_s(z) = \sum_{k \ge 1} z^k / k^s$, continued to the whole plane, at a scalar real order given as a PolylogOrder.

The order is uniform across the packet and tagged by class, because whole-number order is a different, far cheaper algorithm than arbitrary real order and every order-dependent coefficient is a per-call scalar precompute. See the order module for why it is not a vector. Integer covers both signs: $n \le 0$ is the closed rational form (a polynomial in $z/(1-z)$), $n = 1$ is $-\ln(1-z)$, and $n \ge 2$ runs entirely on tabulated $\zeta$ values. Real is the general algorithm (Wood 1992, Roughan 2026): the defining series, the unity series about $z = 1$ with its two cancelling poles fused algebraically so orders arbitrarily close to an integer cost nothing extra, and Wood’s m-th-root identity in the far field.

On a real vector the argument is real and the result is the real part of the principal value, which for $z > 1$ (the cut) is the same from either side. Complex vectors return the full value. On the cut it follows the sign of $\mathrm{Im}\,z$’s zero, C99 style, with -0 giving mpmath’s and Wood’s convention for a bare real.

let li2 = z.polylog(PolylogOrder::Integer(2));   // the dilogarithm
let fd  = (-x.exp()).polylog(PolylogOrder::Real(1.5)); // -F_{1/2}(x)/Gamma(3/2)

The order is spelled in the vector’s own element types: Real carries Self::Element (a complex element on a complex vector, of which only a real value is implemented and anything else answers NaN, or a dual element on a dual vector, whose derivative part must be zero) and Integer carries the signed lane element (i64 on an f64 vector, i32 on an f32 one). Every order-dependent coefficient is computed once per call in that element type through the scalar math surface.

Special values: $\mathrm{Li}_s(1) = \zeta(s)$ for $s > 1$ and $+\infty$ below, $\mathrm{Li}_s(-1) = -\eta(s)$, $\mathrm{Li}_s(0) = 0$. Every arm is a fixed-length series whose length follows the policy’s precision tier. Whole-number orders past $n = 79$ (binary64) or $n = 34$ (binary32, where $n!$ overflows) return NaN in the far field ($|\ln z| > 3.2$). The series and unity arms have no such limit. Cost grows with $\ln|z|$ in the far field at real order (one unity series per root, $m \approx \ln|z| / 2.08$ roots).

Measured against mpmath on 4952 points (real and complex $z$, orders from -6 to 30 and a dozen real ones including $2 + 10^{-9}$), binary64 at Precision: whole-number orders $n \ge 0$ within 1.3e-14 relative on the real line. Negative whole orders within 1.5e-13 (the alternating defining series on the negative axis peaks at ~2500x its sum). Real orders within 3.1e-13, with the far field’s m-th-root sum cancelling by $m^{s-1}$, which is what makes binary32 real order 1.1e-4 there and 2e-5 elsewhere. On the cut the real part is accurate normwise (the imaginary part can be a millionth of it near $z = 1$ at $s = 1 + 10^{-6}$).

Autodiff closes by $\mathrm{Li}_s'(z) = \mathrm{Li}_{s-1}(z)/z$ with the order lowered by one, which is why the runtime order is what the trait carries.

fn scalar_bessel_n<F: BesselFamily, const N: i32>(self) -> Self

A cylindrical Bessel function at compile-time order, selected by family marker: J, Y, I, K, or any of them under Scaled. N is signed and the families reflect at negative order ($J_{-n} = (-1)^n J_n$, $I_{-n} = I_n$).

let j2 = x.bessel_n::<J, 2>();                 // J_2(x)
let ke = x.bessel_n::<Scaled<K>, 0>();         // e^x K_0(x)

The marker only selects: each spelling is a one-line route into the kernel for that family, scaling and order form, with nothing evaluated that was not asked for. Scaled<J> and Scaled<Y> are the SciPy jve/yve scalings by $e^{-|\mathrm{Im}\,z|}$, which is 1 on the real axis, so on a real vector they are J and Y unchanged. On a complex vector they are the scaled values.

fn scalar_bessel<F: BesselFamily>( self, order: BesselOrder<Self, Self::Signed>, ) -> Self

bessel_n with the order taken per lane, at runtime, as a BesselOrder of any class.

let iv = x.bessel::<Scaled<I>>(BesselOrder::Real(nu));   // e^{-|x|} I_nu(x)
let jh = x.bessel::<J>(BesselOrder::HalfInteger(k));     // J_{k/2}(x), elementary

fn scalar_sph_bessel_n<F: BesselFamily, const N: usize>(self) -> Self

A spherical Bessel function at compile-time order, the twin of bessel_n for $j_n$, $y_n$, $i_n$, $k_n$.

let j3 = x.sph_bessel_n::<J, 3>();             // j_3(x)
let ke = x.sph_bessel_n::<Scaled<K>, 1>();     // e^x k_1(x)

fn scalar_sph_bessel<F: BesselFamily>(self, n: u32) -> Self

sph_bessel_n for an order known only at runtime.

fn scalar_airy<W: AiryFn>(self) -> Self

One Airy function selected by marker: Ai, AiPrime, Bi, BiPrime, or any of them under Scaled.

Not a slice of airy_all: the four outputs come from two Bessel passes (order 1/3 for the values, 2/3 for the derivatives), and asking for one runs one pass (Ai skips the I half of it too, so it is roughly a quarter of the tuple). Take the tuple when you want more than one of them.

let ai = x.airy::<Ai>();
let bp = x.airy::<Scaled<BiPrime>>();      // e^{-zeta} Bi'(x) on the positive axis

fn scalar_airy_all<const SCALED: bool>(self) -> (Self, Self, Self, Self)

$(\mathrm{Ai}, \mathrm{Ai}', \mathrm{Bi}, \mathrm{Bi}')$, all four, with the exponential factored out on the positive axis when SCALED (SciPy airy / airye).

Prefer the scaled form on accuracy grounds, not only range: on the positive axis the kernel produces $e^{\zeta}K$ natively, so it evaluates no exponential anywhere and holds 1-3 eps where the unscaled one reaches 684 at x = 100 ($\zeta = \tfrac{2}{3}x^{3/2}$). Unscaled, Ai underflows past x ~ 104 and Bi overflows past x ~ 104.5. For x < 0 the functions oscillate, nothing is factored out, and the phase error grows like $|x|^{3/2}$ in every library.

fn scalar_beta(self, y: Self) -> Self

Computes the Beta function $\mathrm{B}(x, y)$

fn scalar_lbeta(self, y: Self) -> Self

Computes $\ln\left|\mathrm{B}(x, y)\right|$, the log of the absolute Beta function.

beta itself underflows to zero for quite ordinary arguments ($\mathrm{B}(200, 200)$ is about 1e-121, already gone in f32) and overflows for arguments straddling the poles. The log form has range to spare in both directions and is what the surrounding computation usually wants anyway, since Beta almost always appears inside a product of Gammas that is about to be logged.

Evaluated as $\ln\Gamma(x) + \ln\Gamma(y) - \ln\Gamma(x+y)$. The absolute value follows lgamma, so recover the sign from lgamma_r if the arguments can be negative.

This buys range at some cost in relative accuracy. The three lgamma terms cancel against each other, shedding roughly $\log_{10}\frac{\ln\Gamma(x+y)}{|\ln \mathrm{B}|}$ digits. That is under one digit at $x = y = 200$, and a little over two at $x = 200,\ y = 1$ where the terms are near 860 and the answer is near -5.3. It remains far better conditioned than beta, which simply has no value to return across most of that domain.

fn scalar_jacobi(self, alpha: Self, beta: Self, n: u32, m: u32) -> Self

Computes the m-th derivative of the n-th degree Jacobi polynomial

A the special case where α and β are both zero, the Jacobi polynomial reduces to a Legendre polynomial.

NOTE: Given constant α, β or n, LLVM will happily optimize those away and unroll loops.

fn scalar_hermite_n<const N: usize>(self) -> Self

Computes the N-th degree physicists’ Hermite polynomial $H_N(x)$ where x is self and N is the polynomial degree.

Evaluated by the three-term recurrence

H_{n+1}(x) = 2x\,H_n(x) - 2n\,H_{n-1}(x)

seeded with $H_0 = 1$ and $H_1(x) = 2x$. The trip count is N, with no data dependence, so LLVM unrolls the whole thing into straight-line FMA.

The derivative is another member of the same family, $H_n'(x) = 2n\,H_{n-1}(x)$, so a value-and-slope pair costs one extra call rather than a separate kernel. The probabilists’ polynomials are a rescaling, $He_n(x) = 2^{-n/2} H_n(x/\sqrt{2})$.

NOTE: this is the raw polynomial, which grows fast: $H_n(0) = (-2)^{n/2} (n-1)!!$ for even n, and $H_n(x) \sim (2x)^n$ in the tails. It leaves binary32 range at the origin around degree 48 and binary64 around 300, and much earlier for |x| of a few units. If what you actually want is the normalized Hermite function (the quantum harmonic oscillator eigenstate, a Hermite-Gauss beam mode, or the basis of a Hermite spectral method), use hermite_function, which folds the Gaussian weight and the normalization into the recurrence and stays $O(1)$ at every degree. The raw polynomial is the right primitive for Gauss-Hermite quadrature node-finding at modest n and for anything that genuinely wants $H_n$ itself.

fn scalar_hermitev(self, n: Self::Unsigned) -> Self

Computes the n-th degree physicists’ Hermite polynomial H_n(x) where x is self and n is a vector of unsigned integers representing the polynomial degree.

The polynomial is calculated independently per-lane with the given degree in n.

This uses the recurrence relation to compute the polynomial iteratively.

fn scalar_hermite(self, n: u32) -> Self

$H_n(x)$ for a degree known only at runtime: hermitev with the degree splatted, which is the cheapest correct spelling of a uniform degree. The runtime twin of hermite_n.

fn scalar_hermite_function_n<const N: usize>(self) -> Self

Computes the orthonormal Hermite function

\psi_N(x) = \frac{1}{\sqrt{2^N N! \sqrt{\pi}}}\, e^{-x^2/2}\, H_N(x)

where x is self. These are the eigenfunctions of the quantum harmonic oscillator and of the Fourier transform, the Hermite-Gauss modes of a paraxial beam, and the basis of Hermite spectral methods. They are orthonormal on the whole line, $\int \psi_m \psi_n\, dx = \delta_{mn}$.

Evaluated by the recurrence on the functions themselves,

\psi_{n+1}(x) = \sqrt{\tfrac{2}{n+1}}\, x\, \psi_n(x) - \sqrt{\tfrac{n}{n+1}}\, \psi_{n-1}(x)

which keeps every intermediate $O(1)$ (the polynomial’s growth and the Gaussian’s decay cancel inside each step), so unlike hermite it does not overflow at high degree. Both square roots are literals under the unrolled loop. The per-step cost is one FMA on the critical path.

Range

The only quantity that can leave the exponent range is the Gaussian seed, which is carried as $e^{-x^2/4}$ in two halves to double the reach. Full accuracy at every degree holds for $|x|$ under about 18.7 (binary32) or 53 (binary64), which covers every degree up to about 175 / 1400 everywhere on the line, since past the turning point $\sqrt{2n+1}$ the true value decays faster than the seed. Beyond that the result is still correct wherever $e^{-x^2/4}$ is representable, and zero past it.

Under a Best-or-better precision policy on true-FMA hardware, the rounding of $x^2$ (which is the entire error budget of a Gaussian at large x) is recovered exactly and corrected to first order.

fn scalar_hermite_function(self, n: u32) -> Self

$\psi_n(x)$ for a degree known only at runtime. The runtime twin of hermite_function_n: the same seed and recurrence, with the per-step constants computed rather than folded.

fn scalar_hermite_function_series_n<const N: usize>( self, coeffs: &[Self::Element; N], ) -> Self

Evaluates a finite series of Hermite functions at x = self:

\sum_{k=0}^{N-1} \mathrm{coeffs}[k] \cdot \psi_k(x)

with $\psi_k$ as in hermite_function. Evaluated by Clenshaw’s backward recurrence, which is more stable than summing the functions one at a time and never forms them individually. N is the length of the coefficient array, so the highest function is $\psi_{N-1}$; N = 0 is rejected.

Same range as hermite_function: the coefficients are pre-scaled by half of the Gaussian and the outer factor carries the other half, so the running Clenshaw values grow no faster than $e^{x^2/4}$.

fn scalar_hermite_function_series(self, coeffs: &[Self::Element]) -> Self

hermite_function_series_n over a runtime-length coefficient slice.

Same recurrence, same pre-scaling, same range. The length is the only difference, and it costs real work rather than only unrolling: the recurrence coefficients $\sqrt{2/(k+1)}$ and $\sqrt{k/(k+1)}$ fold to literals when N is a constant and become per-step square roots when it is not. Prefer the const form when the degree is known.

An empty coefficient slice is 0, where the const form rejects N = 0 at compile time.

fn scalar_laguerre_n<const N: usize>(self, alpha: Self) -> Self

Computes the generalized (associated) Laguerre polynomial $L_N^{(\alpha)}(x)$, where x is self and N is the polynomial degree.

Passing alpha = Self::ZERO gives the ordinary Laguerre polynomial $L_N(x)$; because alpha is an ordinary argument rather than a const generic, that case folds away completely when the zero is visible at the call site.

Evaluated by the three-term recurrence

(n+1)\,L_{n+1}^{(\alpha)}(x) = (2n + \alpha + 1 - x)\,L_n^{(\alpha)}(x) - (n + \alpha)\,L_{n-1}^{(\alpha)}(x)

seeded with $L_0^{(\alpha)} = 1$ and $L_1^{(\alpha)}(x) = 1 + \alpha - x$. The trip count is N, with no data dependence, so LLVM unrolls the whole thing into straight-line FMA.

The derivative is another member of the same family, $\frac{\mathrm{d}}{\mathrm{d}x} L_n^{(\alpha)}(x) = -L_{n-1}^{(\alpha+1)}(x)$, so a value-and-slope pair costs one extra call rather than a separate kernel.

NOTE: the forward recurrence is the standard evaluation route (Boost and GSL both use it) and is well behaved across the oscillatory region $0 \le x \lesssim 4n$. Past that $L_n^{(\alpha)}$ itself grows like $(-x)^n/n!$ and will overflow for large N and x on its own account.

Laguerre-Gaussian beam modes, the radial part of the hydrogen wavefunction, the quantum harmonic oscillator and coherent-state expansions, and Gauss-Laguerre quadrature.

fn scalar_laguerrev(self, alpha: Self, n: Self::Unsigned) -> Self

Computes the generalized (associated) Laguerre polynomial $L_n^{(\alpha)}(x)$ where n is a vector of unsigned integers giving the degree per lane.

The per-lane counterpart of laguerre, in the same relation to it as hermitev is to hermite. The recurrence runs to the largest n in the vector and lanes freeze at their own degree, so the cost is set by max(n) rather than by any one lane.

fn scalar_laguerre(self, alpha: Self, n: u32) -> Self

$L_n^{(\alpha)}(x)$ for a degree known only at runtime: laguerrev with the degree splatted. The runtime twin of laguerre_n.

fn scalar_laguerre_function_n<const N: usize>(self, alpha: Self) -> Self

Computes the orthonormal generalized Laguerre function

l_N^{(\alpha)}(x) = \sqrt{\frac{N!}{\Gamma(N+\alpha+1)}}\; x^{\alpha/2} e^{-x/2}\, L_N^{(\alpha)}(x)

where x is self. Orthonormal on the half-line, $\int_0^\infty l_m l_n\, dx = \delta_{mn}$. This is the radial factor of Laguerre-Gauss beam modes and (up to a power of x from the spherical measure) of the hydrogen wavefunctions. Defined for $x \ge 0$ and $\alpha > -1$, and nothing is checked outside that.

Evaluated by the recurrence on the functions themselves, with $s_k = \sqrt{(k+1)(k+\alpha+1)}$:

l_{k+1} = \frac{(2k + \alpha + 1 - x)\, l_k - s_{k-1}\, l_{k-1}}{s_k}

which keeps every intermediate $O(1)$, so unlike laguerre it does not overflow at high degree or large x. alpha is a runtime vector, so each step also carries a sqrt and a reciprocal, beside the recurrence rather than on its critical path, and folded to literals when alpha is a visible constant. The seed is skipped outright by a uniform branch when every lane has alpha = 0, which is the ordinary Laguerre function and by far the common case.

Range

The Gaussian-like seed $x^{\alpha/2} e^{-x/2}$ is carried as $e^{-x/4}$ in two halves, as in hermite_function. Full accuracy at every degree for x under about 350 (binary32) or 2800 (binary64), covering every degree up to roughly 87 / 700 everywhere on the half-line (the turning point of $l_n^{(\alpha)}$ is near 4n).

alpha is unrestricted over the same x range. The seed’s whole parameter dependence, $x^{\alpha/2}/\sqrt{\Gamma(\alpha+1)}$, is the square root of the Poisson mass $P(\alpha; x)$ and is evaluated as poisson_pmf is (Loader’s saddle-point form, one exponential of a small exponent), so neither factor materializes (separately $x^{\alpha/2}$ overflows binary64 near $\alpha = 250$ and $1/\sqrt{\Gamma(\alpha+1)}$ underflows near $\alpha = 320$, and their overlap would be inf * 0) and nothing large is exponentiated: 0-3 ulp at the peak x ~ alpha out to $\alpha = 1400$, against a 50-digit oracle.

fn scalar_laguerre_function(self, alpha: Self, n: u32) -> Self

$\ell_n^{(\alpha)}(x)$ for a degree known only at runtime. The runtime twin of laguerre_function_n: the same seed and recurrence, with the per-step scales computed rather than folded.

fn scalar_laguerre_function_i_n<const N: usize>(self, alpha: i32) -> Self

laguerre_function at an integer weight, taken as a scalar i32 rather than a vector.

Same function and same range. What changes is what the compiler can see. Every quantity the recurrence derives from the weight (the $s_k = \sqrt{(k+1)(k+\alpha+1)}$ and their reciprocals, and the $2k+\alpha+1$ offsets) becomes a scalar constant instead of a vector sqrt and reciprocal per step, and folds to a literal outright when alpha is compile-time known.

The seed changes too. Up to $\alpha = 170$ (binary64) / 29 (binary32) the normalization $x^{\alpha/2}/\sqrt{\alpha!}$ is a scalar factorial, a powi and at most one sqrt, with no ln, lgamma or second exp at all, and a few ulp more accurate than the log form, whose lgamma error is amplified by the exponential. $\alpha = 0$ is a scalar test that skips even that. Beyond the cap it takes the vector form’s saddle-point seed. Measured on AVX2 f64x4 at degree 4: about 5x faster than the vector form at a literal small weight, 2x at a runtime one.

Prefer this whenever the weight is a non-negative integer, which every classical application has: the hydrogen radial functions use $\alpha = 2\ell+1$ and the Laguerre-Gauss beam modes use $\alpha = |\ell|$. Negative values are out of domain, as $\alpha \le -1$ is for the general form.

Like the series forms this is inlined into the caller rather than given its own dispatch trampoline: the weight is a plain i32 argument, and a shared out-of-line copy would take it at runtime, which both defeats the folding above and (measured) stops LLVM overlapping consecutive evaluations, at 7x the cost. Call it from inside a #[thermite::dispatch] body.

fn scalar_laguerre_function_i(self, alpha: i32, n: u32) -> Self

laguerre_function_i_n for a degree known only at runtime.

fn scalar_laguerre_function_series_n<const N: usize>( self, alpha: Self, coeffs: &[Self::Element; N], ) -> Self

Evaluates a finite series of generalized Laguerre functions at x = self:

\sum_{k=0}^{N-1} \mathrm{coeffs}[k] \cdot l_k^{(\alpha)}(x)

with $l_k^{(\alpha)}$ as in laguerre_function. Clenshaw’s backward recurrence, same range as the single function; N is the coefficient count and N = 0 is rejected.

fn scalar_laguerre_function_series( self, alpha: Self, coeffs: &[Self::Element], ) -> Self

laguerre_function_series_n over a runtime-length coefficient slice.

Same recurrence, same pre-scaling, same range. The per-step weights are computed rather than folded, as in hermite_function_series. An empty coefficient slice is 0.

fn scalar_laguerre_function_series_i_n<const N: usize>( self, alpha: i32, coeffs: &[Self::Element; N], ) -> Self

laguerre_function_series at a scalar integer weight, in the same relation to it as laguerre_function_i is to laguerre_function. See there for what the integer form buys.

fn scalar_laguerre_function_series_i( self, alpha: i32, coeffs: &[Self::Element], ) -> Self

laguerre_function_series_i_n over a runtime-length coefficient slice.

The _n is the coefficient count and the _i is the integer weight, in that order because the length is the newer axis, and both mean what they do everywhere else. An empty coefficient slice is 0.

fn scalar_chebyshev_n<const K: usize, const N: usize>( self, coeffs: &[Self::Element; N], ) -> Self

Evaluates a finite series of Chebyshev polynomials of the K-th kind at x = self:

\sum_{k=0}^{N-1} \mathrm{coeffs}[k] \cdot P_k(x)

where P_k is T_k, U_k, V_k, or W_k depending on K. All four kinds share the recurrence $P_{k+1}(x) = 2x \cdot P_k(x) - P_{k-1}(x)$ with P_0(x) = 1, and differ only in P_1(x):

K	Kind	P_1(x)	Notes
1	First (T_k)	x	Most common, the minimax/approximation basis on [-1, 1].
2	Second (U_k)	2x	Related to $\sin((k+1)\theta)/\sin(\theta)$ under $x = \cos\theta$.
3	Third (V_k)	2x - 1	“Airfoil” polynomials; $\cos((k+\tfrac12)\theta)/\cos(\theta/2)$.
4	Fourth (W_k)	2x + 1	$\sin((k+\tfrac12)\theta)/\sin(\theta/2)$.

Any other value of K is a compile-time error.

There is deliberately no single-polynomial T_n(x) entry point beside this, unlike legendre or hermite. Chebyshev polynomials are used almost exclusively as an approximation basis, i.e. as a series; their quadrature nodes and weights are closed-form, so nothing needs to iterate on a lone $T_n$; and the one genuine single-$T_n$ application (Chebyshev filter response, Dolph-Chebyshev windows) needs $|x| > 1$, where the right evaluation is $\cosh(n \cosh^{-1} x)$ and not this recurrence at all. A unit coefficient array recovers $T_n$ if it is ever wanted.

Evaluation is done via Clenshaw’s backward recurrence with FMA, which is more numerically stable than a forward sum when the partial sums of $\sum c_k P_k$ are much smaller than $\max_k |c_k P_k|$ (e.g. fitted minimax series with alternating-sign coefficients). N is the length of the coefficient slice, so the highest polynomial term is P_{N-1}; N = 0 is rejected, N = 1 evaluates to coeffs[0].

coeffs[0] multiplies P_0 = 1, coeffs[1] multiplies P_1(x) (which depends on K), and so on. Because LLVM sees both K and N as constants, the recurrence loop and the P_1 selection are fully unrolled and specialized at monomorphization time.

Accuracy near $x = \pm 1$

The plain recurrence forms $2x b_{k+1} - b_{k+2}$ with consecutive $b_k$ of nearly equal magnitude as x approaches either endpoint, and cancels. This is a property of the recurrence, not of the series: measured against a 60-digit oracle at N = 24, it costs up to 37 ulp on sums whose own condition number is about 1, and up to 230 ulp on unstructured coefficients.

Under a Best-or-better precision policy, real vectors instead take Reinsch’s modification, which recurs on the differences (near +1) or sums (near -1) so the small quantity is never formed by subtraction. On the same grid that bounds the error envelope 2.5x to 17x tighter across all four kinds. It is an envelope improvement rather than a pointwise one (individual arguments can land worse), and costs roughly 2x on the recurrence’s dependency chain, which is why it is gated.

binary32 gains the same way, 2.6x to 13.5x on its own grid. Measuring it needs an f32-native one: 1 - 2^-j rounds to exactly 1.0 for every j >= 24, so an f64 grid piles two thirds of its points onto the endpoint itself, where the endpoint form degenerates into a plain running sum and the two policies agree, and never samples the f32 neighbourhood where the cancellation actually bites.

Coefficients from a minimax or least-squares fit decay geometrically and barely notice either way (about 3 ulp to 1). The gap opens on slowly-decaying or non-decaying spectra: truncated expansions, near-singular functions, or coefficients that came from somewhere other than a fit.

Complex and the composite arithmetics keep the plain recurrence at every policy, since Reinsch needs a real copysign and a meaningful nearest endpoint.

fn scalar_chebyshev<const K: usize>(self, coeffs: &[Self::Element]) -> Self

chebyshev_n over a runtime-length coefficient slice.

K stays a const generic, since it selects which Chebyshev kind, not how many coefficients, and there are exactly four. Only the length becomes dynamic.

Same recurrence and the same Best-precision Reinsch form near $x = \pm 1$; what the runtime length costs is the unrolling and the folded coeffs indices. An empty coefficient slice is 0.

fn scalar_gaussian(self, a: Self, c: Self) -> Self

Computes the Gaussian function with amplitude a and standard deviation c, defined as $a\, e^{-\frac{1}{2}(x/c)^2}$.

The position b is assumed to be zero. For a non-zero position, use self - b as the input.

fn scalar_planck(self) -> Self

Computes the Planck shape factor $\frac{x^3}{e^x - 1}$, finite at x = 0 where it vanishes like $x^2$.

The dimensionless kernel of Planck’s law: substituting $x = h\nu/kT$ recovers the spectral radiance up to a scale factor, so this is the part worth computing carefully and the constants are left to the caller. Radiative transfer, climate radiation budgets, and stellar atmospheres.

The denominator cancels for small x and the quotient is $0/0$ at the origin. Evaluated here as $x^2/\varphi_1(x)$ using phi_n::<1>, which is finite and equal to 1 there, so the singularity never forms rather than being patched after the fact.

fn scalar_legendre(self, n: u32, m: u32) -> Self

Computes the m-th associated n-th degree Legendre polynomial, where m=0 signifies the regular n-th degree Legendre polynomial.

If m is odd, the input is only valid between -1 and 1

NOTE: Given constant n and/or m, LLVM will happily unroll and optimize inner loops.

Internally, this is computed with jacobi when m > 0.

fn scalar_legendre_series_n<const N: usize>( self, coeffs: &[Self::Element; N], ) -> Self

Evaluates a finite Legendre series at x = self:

\sum_{k=0}^{N-1} \mathrm{coeffs}[k] \cdot P_k(x)

The form a Legendre-moment expansion takes: Mie and Henyey-Greenstein scattering phase functions tabulated by their moments, multipole expansions in $\cos\theta$, and the polar factor of a spherical-harmonic expansion at fixed order.

Evaluated by Clenshaw’s backward recurrence on the Legendre three-term relation, which is more stable than building each $P_k$ with legendre and summing, and does $O(N)$ work rather than $O(N^2)$. The recurrence ratios $(2k+1)/(k+1)$ and $k/(k+1)$ are literals under the unrolled loop, so the per-step cost matches chebyshev: one FMA on the critical path. N is the coefficient count; N = 0 is rejected, N = 1 evaluates to coeffs[0].

Plain Clenshaw at every policy: the endpoint cancellation that chebyshev treats under Best precision exists here too ($P_n(1) = 1$ for every n), but its Reinsch-style rewrite for the Legendre ratios has not been derived or measured.

fn scalar_legendre_series(self, coeffs: &[Self::Element]) -> Self

legendre_series_n over a runtime-length coefficient slice.

Plain Clenshaw here too. The recurrence ratios $(2k+1)/(k+1)$ and $k/(k+1)$ are literals only when N is a constant, so this pays a division per step where the const form pays none, the widest const-versus-slice gap of the series family. An empty coefficient slice is 0.

fn scalar_zernike_r(self, n: u32, m: u32) -> Self

Computes the Zernike radial polynomial $R_n^m(\rho)$, where rho is self.

Returns zero unless $m \le n$ with $n - m$ even, the condition for the mode to exist. m is the absolute azimuthal frequency here. The sign only affects the angular factor, which lives in zernike.

Evaluated through the shifted Jacobi identity

R_n^m(\rho) = \rho^m\, P_{(n-m)/2}^{(0,\,m)}\!\left(2\rho^2 - 1\right)

rather than the textbook sum $\sum_k (-1)^k \frac{(n-k)!}{k!\,((n+m)/2 - k)!\,((n-m)/2 - k)!} \rho^{n-2k}$, which alternates factorials of size $(n-k)!$ against an answer bounded by 1 and loses all precision somewhere around n = 10-15. That is well inside the range adaptive optics, ophthalmology and surface metrology actually use.

The $(-1)^{(n-m)/2}$ prefactor usually seen with this identity is absent because the argument is written $2\rho^2 - 1$ rather than $1 - 2\rho^2$: reflecting a Jacobi polynomial swaps its two parameters and absorbs exactly that sign.

The polynomial is only orthogonal on $\rho \in [0, 1]$ and grows quickly outside it. Nothing clamps the argument, so an unnormalized pupil coordinate stays the caller’s problem.

fn scalar_zernike<const NORM: u8>(self, theta: Self, n: u32, m: i32) -> Self

Computes the Zernike polynomial $Z_n^m(\rho, \theta)$ on the unit disc, with rho as self:

Z_n^m(\rho, \theta) = N_n^m\, R_n^{|m|}(\rho) \times
  \begin{cases} \cos(m\theta) & m \ge 0 \\ \sin(|m|\theta) & m < 0 \end{cases}

Returns zero unless $|m| \le n$ with $n - |m|$ even.

NORM selects the normalization $N_n^m$ and must be either ZERNIKE_UNIT_PEAK ($N = 1$, so $R_n^m(1) = 1$ and coefficients read as peak amplitude) or ZERNIKE_ORTHONORMAL ($N_n^m = \sqrt{2(n+1)/(1 + \delta_{m,0})}$, the ANSI Z80.28 and Noll convention, under which coefficients read as RMS contributions). Any other value is a compile-time error. There is deliberately no default: the two differ by a factor of up to $\sqrt{2(n+1)}$ per mode, and picking one silently is how coefficient sets get misinterpreted.

(n, m) is a runtime pair rather than a const generic on purpose. The workload is a basis, not a function. A wavefront fit evaluates tens to hundreds of modes over thousands of pupil samples, with the mode list coming from a config or a sensor geometry, so the degree is loop-invariant across the vector axis and const-generic specialization would buy a jump table rather than an unrolled loop.

The single-index conventions (ANSI Z80.28 / OSA, Noll, Fringe) and the conversions between them are in crate::zernike. They disagree from the second term onward, so convert at the boundary rather than assuming.

fn scalar_zernike_basis<const L: usize, const NORM: u8, const N: usize>( x: Self, y: Self, out: &mut [Self; N], )

Evaluates all Zernike modes through degree L at the Cartesian pupil point (x, y), into out[j] for the ANSI Z80.28 / OSA index $j = (n(n+2) + m)/2$.

N must equal (L+1)(L+2)/2 (compile-time checked), and NORM is ZERNIKE_UNIT_PEAK or ZERNIKE_ORTHONORMAL as on zernike.

This is the entry point a wavefront fit or reconstruction wants. It is not merely a loop over zernike. Substituting $s = x^2+y^2$ splits every mode into a polynomial in s times $\operatorname{Re}$ or $\operatorname{Im}$ of $(x+iy)^{|m|}$, which is where the $\rho^{|m|}$ and the $\cos m\theta$ both come from at once. Evaluation is then pure polynomial arithmetic: no atan2, no sqrt, no trigonometry, no division, $O(L^2)$ FMAs for the entire basis, and no singularity at the pupil centre. Calling the single-mode form per mode instead costs a sin_cos and a powi each and restarts the radial recurrence every time, for $O(L^3)$ work.

Cartesian input is part of that, not a convenience: pupil samples arrive as (x, y), and a polar entry point would charge an atan2 per sample for an angle this kernel immediately dissolves.

Fully unrolled at compile time for each L up to MAX_ZERNIKE_DEGREE; above that it takes a rolled path that is correct at any degree and substantially slower.

Nothing normalizes (x, y) onto the unit disc. Outside it the polynomials still evaluate correctly and simply are not orthogonal.

The layout is ANSI because it is the scheme whose index has a closed form and whose degree truncation is contiguous. Noll and Fringe callers gather through noll_to_ansi / fringe_to_ansi.

use thermite::prelude::*;
use thermite_special::{SpecialMath, ZERNIKE_ORTHONORMAL};
use thermite_special::zernike::noll_to_ansi;

type V = Vector<f64>;
const L: usize = 4;
const N: usize = 15; // (L+1)(L+2)/2

let mut basis = [V::ZERO; N];
V::zernike_basis::<L, ZERNIKE_ORTHONORMAL, N>(V::splat(0.3), V::splat(0.4), &mut basis);

// Noll 4 is defocus, Z_2^0 = sqrt(3) (2 rho^2 - 1) orthonormal.
let defocus = basis[noll_to_ansi(4) as usize].extract::<0>();
assert!((defocus - 3f64.sqrt() * (2.0 * 0.25 - 1.0)).abs() < 1e-14);

fn scalar_lambert_w(self) -> (Self, Self)

Computes both branches of the Lambert W function simultaneously: ($W_0(x)$, $W_{-1}(x)$).

The $W_0$ result is valid for x >= -1/e; the $W_{-1}$ result is valid for -1/e <= x < 0. Outside these domains, the respective result is NaN (when overflow checking is enabled).

fn scalar_expint_n<const N: usize>(self) -> Self

Computes the generalized exponential integral E_n(x) for integer order n.

fn scalar_expint(self, n: u32) -> Self

E_n(x) for an order known only at runtime. The runtime twin of expint_n: the same E_1 kernel, the same recurrence and the same continued-fraction handover, so the two agree to the bit.

fn scalar_phi_n<const N: usize>(self) -> Self

Returns $\varphi_N(x)$, the N-th phi-function of exponential integrators.

\varphi_0(x) = e^x, \qquad
\varphi_{k+1}(x) = \frac{\varphi_k(x) - 1/k!}{x}, \qquad
\varphi_k(x) = \sum_{n \ge 0} \frac{x^n}{(n + k)!}, \qquad
\varphi_k(0) = \frac{1}{k!}

phi_n::<0> is exp. phi_n::<1> is $(e^x - 1)/x$, which written out directly is $0/0$ at the origin and loses most of the mantissa near it, so it is evaluated as $\mathrm{expm1}(x)/x$ with the removable singularity filled in (the value is 1), which is accurate across the whole line. Outside the exponential-integrator literature phi_n::<1> goes by exprel, which is the name SciPy, Boost and the statistics literature use for it. There is no separate exprel here because this is it. Beyond that the recurrence is the wrong way to compute them: each step subtracts 1/k! from a value that is barely larger while |x| is small, so $\varphi_2 = (\mathrm{expm1}(x) - x)/x^2$ loses twice the bits phi_n::<1> would have, and gets worse with N. Below |x| = N this sums the series instead (its terms are monotone there, so nothing cancels), and above it runs the recurrence upward from expm1, where the amplification per step is bounded. Measured against mpmath, both arms sit within a few ulp for N <= 8.

The series arm’s length is bounded by the policy’s max_iterations. The primitive float types know their precision statically and use a fixed count instead. Nothing caps N, though nothing needs it large: ETDRK4 wants phi_1..phi_3, and exponential Rosenbrock methods rarely go past phi_4.

phi_n::<1> alone is the coefficient that keeps appearing wherever an exponential is integrated over a finite step:

• The in-scattering integral through a homogeneous medium, $\int_0^t e^{-\sigma s}\,ds = t\,\varphi_1(-\sigma t)$. The singular case is the empty medium, which is not an edge case in practice.
• Exact stepping of an Ornstein-Uhlenbeck process, and the Langevin thermostat’s mean-reversion factor.
• Frame-rate-independent exponential smoothing, usually written 1 - exp(-k * dt) and then divided by k.

The higher orders are the coefficients of exponential time differencing: integrating y' = Ly + N(y) exactly over a step gives $y(h) = e^{hL} y_0 + h\,\varphi_1(hL)\,N$, and expanding N in time along the step brings in $\varphi_2, \varphi_3, \ldots$ as the weights of the higher-order terms.

fn scalar_phi(self, n: u32) -> Self

$\varphi_n(x)$ for an order known only at runtime. The runtime twin of phi_n: the same series and recurrence arms, with the series length worked out from n per call rather than at compile time.

fn scalar_carlson<K: WrapTo>(kind: K) -> Self

where K::Wrapped: CarlsonKind, <K::Wrapped as CarlsonKind>::Output: Unwrap<Unwrapped = Self>,

Carlson symmetric elliptic integral, selected by a CarlsonKind request struct with named fields. The arity (and which argument is the parameter / repeated one) is fixed per kind, so the wrong shape is a compile error.

let rf = V::carlson(CarlsonRf { x, y, z });
let rj = V::carlson_p::<Precision, _>(CarlsonRj { x, y, z, p });

fn scalar_ellint<K: WrapTo>(kind: K) -> Self

where K::Wrapped: EllipticKind, <K::Wrapped as EllipticKind>::Output: Unwrap<Unwrapped = Self>,

Legendre elliptic integral, selected by an EllipticKind request struct. Each form (EllintK/EllintF/EllintE/ EllintEInc/EllintD/EllintDInc/ EllintPi/EllintPiInc) carries exactly its own arguments, and completeness is encoded by whether the struct has a phi field.

Two family members that are not Legendre integrals dispatch through here as well, because they are built from the same Carlson forms and belong beside their siblings: JacobiZeta, the oscillating part of $E(\varphi, k)$, and HeumanLambda, its complementary-modulus companion.

let k_int = V::ellint(EllintK { k });                       // K(k)
let e_inc = V::ellint_p::<Precision, _>(EllintEInc { phi, k }); // E(phi, k)
let z     = V::ellint(JacobiZeta { phi, k });               // Z(phi, k)

fn scalar_erfinv(self) -> Self

Computes the inverse error function.

fn scalar_probit(self) -> Self

Computes the Probit function, the inverse of the cumulative distribution function of the standard normal distribution.

fn scalar_ndtr(self) -> Self

Computes the cumulative distribution function of the standard normal distribution, the inverse of probit:

\Phi(x) = \frac{1}{\sqrt{2\pi}} \int_{-\infty}^{x} e^{-t^2/2}\,dt
        = \tfrac12 \operatorname{erfc}\!\left(-\frac{x}{\sqrt 2}\right)

The probability that a standard normal variable falls below x: z-scores to p-values, the N(d_1)/N(d_2) terms of Black-Scholes, the probit link, and x * ndtr(x) is GELU. The name is Cephes/SciPy’s.

Underflows to zero below about x = -38.6 (f64) and -14.4 (f32). When the tail probability itself is the quantity of interest, use log_ndtr, which is finite there.

fn scalar_log_ndtr(self) -> Self

Computes $\ln \Phi(x)$, the logarithm of the standard normal CDF, finite for every finite x.

ln(ndtr(x)) is -inf below x ~ -38.6 in f64 (-14.4 in f32), exactly where a probit or censored-regression likelihood, a truncated-normal density, or an expected-improvement acquisition needs the tail: log_ndtr(-100) is an ordinary -5004.6. The kernel keeps $-x^2/2$ in the log domain and takes the rest from erfcx, which has no underflow, so the left tail carries full relative accuracy to the largest x whose square is representable. On the right it is ln_1p of the complement, so log_ndtr(10) = -7.6e-24 rather than a rounded zero.

Costs one erfcx, one ln_1p, and an exp for the lanes with x > 0.

fn scalar_logerfc(self) -> Self

Computes $\ln \operatorname{erfc}(x)$, finite for every finite x.

erfc underflows at x ~ 27 (f64) / 9.3 (f32) and its logarithm does not: logerfc(100) = -10004.8. This is the log-domain form of a Gaussian tail wherever erfc rather than the normal CDF is the natural quantity (Ewald sums, Gaussian-smeared edges, the Mills ratio in the log domain), and it is log_ndtr with x = -\sqrt 2 x'. Built on erfcx with $-x^2$ kept in the log domain. On the left, where erfc(x) is between 1 and 2, it is ln_1p(erf(|x|)), so the result stays accurate down to logerfc(-1e-20) = 1.13e-20.

fn scalar_fresnel(self) -> (Self, Self)

The Fresnel integrals $S(x) = \int_0^x \sin(\pi t^2/2)\,dt$ and $C(x) = \int_0^x \cos(\pi t^2/2)\,dt$, together.

Returns (S, C), the same order as SciPy’s fresnel and this crate’s own sici.

Both are odd, both tend to 1/2, and both stay in [0.32, 0.72] past the first oscillation. Measured against a 45-digit oracle over x from 1e-4 to 1e15: 2.80 ulp (C) and 2.64 (S) in f64, 2.14 and 3.40 in f32 out to 1e7.

The phase $\pi x^2/2$ is carried in two words and reduced exactly, which is not a refinement but the whole of the large-argument accuracy: computed the obvious way as x*x*0.5, the phase is already 5.3e-6 wrong at x = 98765 and returns the wrong sign by $x \approx 10^9$, and since C and S are 1/2 plus a term of size $1/(\pi x)$ that error lands straight on the result. Below Average the residual is dropped and that behaviour returns.

Above x = 1.147e16 (f64) / 2.136e7 (f32) the oscillating correction is under half an ulp of 1/2, and both are exactly 1/2.

fn scalar_fresnel_c(self) -> Self

$C(x)$ alone. See fresnel.

Unlike airy::<Ai> this is not a cheaper evaluation by much: C and S share the argument reduction, the phase and both auxiliaries, so asking for one drops a single Chebyshev series and one reconstruction: roughly a third, not three quarters.

fn scalar_fresnel_s(self) -> Self

$S(x)$ alone. See fresnel_c.

fn scalar_sici(self) -> (Self, Self)

The trigonometric integrals $\mathrm{Si}(x) = \int_0^x \frac{\sin t}{t}\,dt$ and $\mathrm{Ci}(x) = \gamma + \ln x + \int_0^x \frac{\cos t - 1}{t}\,dt$, together. Returns (Si, Ci).

Si is odd. Ci is real only on the positive axis ($\mathrm{Ci}(-x) = \mathrm{Ci}(x) + i\pi$), so this returns Ci(|x|), dropping the imaginary part, which is what SciPy’s sici does. Ci(0) is $-\infty$.

Measured 2.03 ulp (Si) and 1.42 (Ci, against its envelope) in f64 over x from 1e-4 to 1e15. In f32, 1.34 and 1.99.

Two things worth knowing before relying on Ci:

• It has zeros, the first near x = 0.6165, and no algorithm is relatively accurate at one. The accuracy above is relative to $\lvert\gamma + \ln x\rvert + \lvert\mathrm{Cin}\rvert$ below the crossover and to the $1/x$ envelope above it.
• Its large-argument accuracy is sin_cos’s: for Ci the oscillation is the value, so a phase error is a relative error, and full argument reduction is a Best-tier property. Si is insulated, tending to $\pi/2$ with the oscillation only a $1/x$ correction, and is $\pi/2$ exactly above x = 1.147e16 (f64) / 2.136e7 (f32). Ci has no such cutoff: it decays like $1/x$ and stays representable for every finite x.

fn scalar_sinint(self) -> Self

$\mathrm{Si}(x)$ alone. See sici, and fresnel_c for what a single accessor saves.

fn scalar_cosint(self) -> Self

$\mathrm{Ci}(x)$ alone. See sici.

fn scalar_inv_log_ndtr(self) -> Self

Computes the inverse of log_ndtr: the x with $\ln \Phi(x) = y$, for y <= 0. The quantile of a log-probability.

probit of $e^y$ stops working once $e^y$ underflows (y < -745 in f64), which is exactly where a log-likelihood, a truncated-normal EM step or an extreme-value fit needs the quantile. This inverts log_ndtr directly, by Newton with the inverse Mills ratio as the derivative, from a probit(e^y) seed one precision tier down where that exists and from the tail asymptotic below. Within a few ulp of the true inverse of the given y over the whole domain. y = 0 gives +inf, y = -inf gives -inf, and y > 0 is NaN.

fn scalar_inv_digamma(self) -> Self

Computes the inverse of the digamma function on $(0, \infty)$: the x with $\psi(x) = y$.

The maximum-likelihood estimate of a gamma shape or a Dirichlet concentration is this function of a mean log. Newton on digamma with trigamma from Minka’s seed ($e^y + 1/2$ above y = -2.22, $-1/(y + \gamma)$ below). Above y = 6 the Stirling series is solved for x directly, since there Newton on digamma cannot see past digamma’s own rounding. +inf maps to +inf and -inf to 0.

fn scalar_wright_omega(self) -> Self

Computes the Wright omega function, the $\omega > 0$ with $\omega + \ln \omega = x$.

This is $W_0(e^x)$, the principal Lambert W of an exponential, evaluated without forming $e^x$: $W_0(e^x)$ overflows past x = 709 where $\omega(x) \approx x - \ln x$ is ordinary. Newton on $\omega + \ln \omega - x$ from a cheap seed per region. Below x = -7 the Lagrange series in $e^x$ is the answer outright.

fn scalar_bessel_ratio<F: BesselRatioFamily>(self, nu: Self) -> Self

Computes the modified Bessel ratio $A_\nu(x) = I_\nu(x) / I_{\nu-1}(x)$ for nu >= 1, odd in x.

With $p = 2\nu$ this is the mean resultant length of a von Mises-Fisher distribution on $S^{p-1}$ at concentration x. nu = 1 is the von Mises circle $I_1/I_0$, and nu = 3/2 is the langevin function. Never forms the two Bessel functions where they would underflow: a series pair for small x, the continued fraction for the ratio in the middle, and the scaled quotient only where x dominates the order. The order is a plain vector, but whole and half-integer orders reach their fast Bessel kernels through the order simplifier.

fn scalar_inv_bessel_ratio<F: BesselRatioFamily>(self, nu: Self) -> Self

Computes the inverse of bessel_ratio: the concentration $\kappa$ with $I_\nu(\kappa)/I_{\nu-1}(\kappa) = r$, for 0 <= r < 1, odd in r.

The maximum-likelihood concentration of a von Mises-Fisher distribution from its observed mean resultant length, in any dimension $p = 2\nu$. Banerjee’s $r(p - r^2)/(1 - r^2)$ seeds a Newton whose derivative is the closed form $1 - A^2 - (2\nu - 1)A/\kappa$, so each step is one ratio evaluation. r = 1 gives +inf, r > 1 NaN.

As r -> 1 the problem itself is ill-conditioned: $\kappa \sim (p-1)/(2(1-r))$, and an ulp of r is a relative $2\kappa\epsilon/(p-1)$ of $\kappa$. The result is the exact inverse of the given r to that extent.

fn scalar_bessel_ratio_1m<F: BesselRatioFamily>(self, nu: Self) -> Self

Computes $1 - A_\nu(x)$, the complement of bessel_ratio, to full relative accuracy where the ratio itself is within an ulp of 1.

1 - bessel::ratio::<I>(x) is gone once $A$ rounds to 1 (x past 1e16 (p-1)/2), and is only accurate to $\epsilon/(1 - A)$ before that. This evaluates the complement directly for x >= 8 nu, from the Hankel expansions at a reduced order and the ratio recurrence walked upward in complement form. $A$ is odd, so $1 - A(-x) = 2 - (1 - A(x))$.

fn scalar_inv_bessel_ratio_1m<F: BesselRatioFamily>(self, nu: Self) -> Self

Computes the inverse of bessel_ratio_1m: the concentration $\kappa$ with $1 - I_\nu(\kappa)/I_{\nu-1}(\kappa) = t$, for 0 < t <= 2 (t = 1 - r).

The complement form of inv_bessel_ratio for nearly concentrated data: $\kappa \sim (p-1)/(2t)$ as t -> 0. This form keeps full relative accuracy there instead of losing $2\kappa\epsilon/(p-1)$ to the rounding of r. It is the inv_langevin_1m move in every dimension. t = 0 gives +inf. t in (1, 2] is a negative r and returns the mirrored $\kappa$.

fn scalar_gauss_legendre(self, n: u32) -> (Self, Self)

Computes the k-th node and weight of the n-point Gauss-Legendre quadrature rule on $[-1, 1]$, with the root index k taken per lane.

The rule integrates every polynomial through degree $2n - 1$ exactly: $\int_{-1}^{1} f \approx \sum_k w_k f(x_k)$, $x_k$ the roots of $P_n$ in descending order (k = 0 is the largest, $x_{n-1-k} = -x_k$) and $w_k = 2 / ((1 - x_k^2) P_n'(x_k)^2)$. The packet is the rule: sweep k over 0..n in packets of consecutive indices and store the two vectors. Every lane runs the same O(n) recurrence, so a packet of roots costs one root.

Tricomi’s $\cos(\pi(k + 3/4)/(n + 1/2))$ seeds a Newton on $P_n$ from the recurrence, and nodes land within a few $\epsilon$ absolute. A non-integer or out-of-range k gives NaN in both.

let n = 16;
for base in (0..n).step_by(V::LANES) {
    let k = V::from_array(core::array::from_fn(|i| (base + i) as f64));
    let (x, w) = k.gauss_legendre(n as u32); // lanes past n - 1 are NaN
}

fn scalar_gauss_hermite(self, n: u32) -> (Self, Self)

Computes the k-th node and weight of the n-point Gauss-Hermite rule, for $\int_{-\infty}^{\infty} f(x) e^{-x^2}\,dx \approx \sum_k w_k f(x_k)$, the root index k per lane (k = 0 the largest root, $x_{n-1-k} = -x_k$).

Same shape as gauss_legendre: a packet of consecutive indices is the rule. Seeded from the WKB phase of the Hermite equation and finished by Newton on $H_n/n!$, whose recurrence stays in range where the raw $H_n$ overflows at degree 48. The weights are the unscaled ones, which reach $e^{-x_k^2}$ at the outer nodes. The scalar factor in them underflows past n = 170 in f64 and n = 40 in f32, which bounds the rule.

fn scalar_gauss_laguerre(self, alpha: Self, n: u32) -> (Self, Self)

Computes the k-th node and weight of the n-point Gauss-Laguerre rule, for $\int_0^{\infty} f(x)\, x^\alpha e^{-x}\,dx \approx \sum_k w_k f(x_k)$, the root index k and alpha > -1 per lane (k = 0 the largest root).

Same shape as gauss_legendre. Seeded from the WKB phase of the Laguerre equation, whose phase count between the turning points carries the Bessel-zero offset on the left and the Airy offset on the right, and finished by Newton on the raw $L_n^\alpha$ with Hildebrand’s weight $\Gamma(n+\alpha+1)/(n!\,x_k\,L_n^{\alpha\prime}(x_k)^2)$. Unscaled weights, which reach $e^{-x_k}$ at the outer nodes. $L_{n-1}$ at the largest root grows like $e^{x/2}$, which bounds the rule near n = 170 in f64 and n = 20 in f32.

fn scalar_pochhammer(self, m: Self) -> Self

Computes the Pochhammer symbol $(z)_m = \dfrac{\Gamma(z+m)}{\Gamma(z)}$.

Combinatorics calls this the rising factorial, and for a non-negative integer m it is exactly the ascending product $z(z+1)\cdots(z+m-1)$. The name here is the special-function one because the function is not restricted to integers: m is any real, which is what the hypergeometric series need and what “factorial” would misdescribe.

Note that the notation $(z)_m$ is ambiguous in the literature: it means the rising factorial in special functions and the falling factorial through much of combinatorics and statistics. This function is the rising one. The falling factorial is pochhammer(z - n + 1, n), and the two are related by $z^{(\bar n)} = (-1)^n (-z)^{(\underline n)}$. Neither is shipped separately, being an argument transform away.

Accuracy

The obvious spelling exp(lgamma(z+m) - lgamma(z)) cancels catastrophically whenever m is small beside z: at z = 1e8, m = 1e-4 it has no correct digits. This does not use it.

At Average precision and above (which includes the default policy), integer m up to 20 in absolute value takes an exact product, 0.00 ulp median and 4.2 worst. That path also covers negative z and returns exact zeros at the poles: $(-2)_3$ is 0.

Below Average it is compiled out and integer m goes through the Stirling difference like anything else, which measures 4.2 ulp median and 172 worst. The difference that shows is the exactness rather than the ulp count: $(3)_1$ comes back as 3.0000000000000018 there, and $(200)_2$ as 40200.00000000002.

Any other m with z and z+m both positive takes a Stirling difference arranged so nothing large is ever subtracted from anything large. Its error is the floor for anything exponentiating a logarithm, tracking $|\ln (z)_m|\cdot\epsilon$. Over 6924 measured points with z in [0.1, 8.9] that is a median of 2.6 ulp and a 99th percentile of 25. Individual points scale with the result’s own logarithm, reaching 259 ulp where the value is near 1e163, and falling to nothing as the result approaches 1.

A non-integer m with z or z+m non-positive (a ratio taken across Gamma’s poles) has no cheap rearrangement and does fall back to the logarithmic form, inheriting its cancellation.

fn scalar_jacobi_elliptic(self, k: Self) -> (Self, Self, Self)

Computes the Jacobi elliptic functions $(\mathrm{sn}, \mathrm{cn}, \mathrm{dn})$ at argument self and modulus k, all three from one evaluation.

All three are made from a single angle, the amplitude $\varphi = \mathrm{am}(u, k)$, defined by $F(\varphi, k) = u$, so this function inverts the incomplete integral of the first kind that ellint evaluates:

\mathrm{sn}(u, k) = \sin\varphi, \qquad
\mathrm{cn}(u, k) = \cos\varphi, \qquad
\mathrm{dn}(u, k) = \sqrt{1 - k^2 \sin^2\varphi}

Hence their names: sine amplitude, cosine amplitude and delta amplitude. At k = 0 the amplitude is u and they collapse to $(\sin u, \cos u, 1)$. At k = 1 they stop being periodic and become $(\tanh u, \operatorname{sech} u, \operatorname{sech} u)$.

Why one function and not three

The triple is closed under differentiation in u, each derivative a product of the other two:

\frac{d\,\mathrm{sn}}{du} = \mathrm{cn}\,\mathrm{dn}, \qquad
\frac{d\,\mathrm{cn}}{du} = -\mathrm{sn}\,\mathrm{dn}, \qquad
\frac{d\,\mathrm{dn}}{du} = -k^2\,\mathrm{sn}\,\mathrm{cn}

so they are one object the way $(\sin, \cos)$ are, and Dual differentiates them without touching the iteration underneath. It also costs nothing to return all three: they share the entire computation, and only the last few operations differ.

The other nine Jacobi functions in Glaisher’s notation (ns, nc, nd, sc, sd, cs, cd, ds, dc) are reciprocals and ratios of these three, so this gives all twelve.

Domain and accuracy

Only $k^2$ enters, so the sign of k does not matter. |k| > 1 is out of domain and gives NaN. Worst absolute error measured against mpmath at 40 digits over |u| <= 8 and k in [0, 1) is 8.3 eps for sn, 4.1 for cn and 3.8 for dn. Absolute is the meaningful metric: all three are bounded by 1 and all three have zeros, so relative accuracy at a zero depends on how well that zero’s location is known, exactly as for sin. For the same reason accuracy falls off slowly with |u|, that being the argument of the single trigonometric call inside.

fn scalar_agm(self, other: Self) -> Self

Computes the arithmetic-geometric mean $\mathrm{AGM}(a, b)$ of two non-negative arguments.

Iterating $a \mapsto (a + b)/2$ against $b \mapsto \sqrt{ab}$ drives the two sequences to a common limit, quadratically: the pair closes to within a factor of a few in a handful of passes from any starting ratio, and the correct digits then double per pass. The loop is branchless and costs one sqrt per iteration, with no transcendentals anywhere, which is why it is also the engine behind the complete elliptic integrals, $K(k) = \pi / (2\,\mathrm{AGM}(1, k'))$, reached through ellint rather than by calling this directly.

Symmetric in its arguments and homogeneous, $\mathrm{AGM}(ca, cb) = c\,\mathrm{AGM}(a, b)$. AGM(a, 0) is 0 and AGM(inf, b) is inf. A negative argument is outside the domain (the geometric mean’s sign becomes ambiguous after the first pass) and returns NaN under overflow checking, as does a zero paired with an infinity.

The geometric mean is formed as one product, so two arguments both above $\sqrt{\text{MAX}}$ (about 1.3e154 in f64, 1.8e19 in f32) overflow to infinity even where the mean is representable. Scale both by a common power of two first if that range matters. Homogeneity makes it exact.

fn scalar_langevin(self) -> Self

Computes the Langevin function $L(x) = \coth x - \frac{1}{x}$.

Odd, strictly increasing, L(0) = 0, L'(0) = 1/3, L(x) -> 1 as x -> ∞. This is the mean resultant length $A_3(\kappa)$ of a von Mises-Fisher distribution on the sphere, and the freely-jointed-chain force-extension law in polymer physics.

Evaluated as an odd minimax polynomial for |x| <= 2 (the direct form coth x - 1/x cancels catastrophically there, losing 3u/x^2), and as 1 - 1/x + 2/(e^{2x} - 1) beyond. Both branches are accurate to a few ulp at every precision policy. The policy mainly selects the exp.

To also obtain the derivative L'(x), use langevin_d.

fn scalar_inv_langevin(self) -> Self

Computes the inverse Langevin function $L^{-1}(y)$ for |y| < 1.

Odd, with a simple pole at y = 1: L^-1(y) ~ 1/(1-y). |y| = 1 returns ±∞, and |y| > 1 returns NaN under overflow checking (an unspecified value otherwise). Its condition number is 1/(1-y), so near the pole the result cannot be more accurate than that, however exact the arithmetic. A consumer that knows 1 - y should form it before rounding.

A rational seed (the same family as Cohen’s Pade approximant, which the vMF literature knows as the Banerjee et al. concentration estimator) is refined by Newton (f32) or Halley (f64) steps whose count follows the precision policy:

precision	steps	relative error
Worst	0	~2e-5
Medium, Average, Best	1	full (a few ulp)
Reference	2	full

fn scalar_langevin_1m(self) -> Self

Computes 1 - L(x), the complement of the Langevin function, accurately where L(x) is within rounding of 1.

1 - L(x) ~ 1/x, so once x > 1/u (sharpness ~1e7 in f32, ~1e16 in f64) langevin(x) rounds to exactly 1 and its complement is gone. This returns it to full relative precision at any x, from the same intermediates. Same cost as langevin. Negative x gives 1 + L(|x|).

Pairs with inv_langevin_1m: the vMF convolution kappa' = L^-1(L(k1) L(k2)) should be formed as inv_langevin_1m(a + b - a*b) with a = langevin_1m(k1), b = langevin_1m(k2), which is cancellation-free at every sharpness.

fn scalar_inv_langevin_1m(self) -> Self

Computes L^-1(1 - t) from the complement t directly.

The inverse Langevin function has a pole at y = 1 and a condition number of 1/(1-y), so a caller that knows 1 - y (see langevin_1m) should pass it here rather than form y and lose its low digits: this entry point works in t throughout and is accurate to a few ulp at any sharpness. t = 0 returns +∞, t > 1 gives the negative branch, and t < 0 is out of the domain (NaN under overflow checking). Same cost as inv_langevin.

fn scalar_gelu(self, alpha: Self) -> Self

GELU activation function, defined as $\tfrac{1}{2} x \left(1 + \operatorname{erf}\!\left(\frac{\alpha x}{\sqrt{2}}\right)\right)$, where alpha helps control the shape of the curve. The standard GELU function is recovered when alpha is 1.

For f32 vectors, this remains decently accurate even with the Medium and Worst precision policies, thanks to good erf implementations at the various precision levels. See erf for more details.

To also obtain the derivative with respect to x (which shares most of the computation), use gelu_d.

fn scalar_swish(self, beta: Self) -> Self

Swish activation function, defined as $x\,\sigma(\beta x) = \frac{x}{1 + e^{-\beta x}}$, where beta controls the sharpness of the gate. The standard Swish/SiLU function is recovered when beta is 1. As beta -> 0, the output approaches x/2 (half-identity); as beta -> inf, Swish approaches ReLU.

To also obtain the derivative with respect to x, use swish_d.

fn scalar_algebraic_sigmoid_n<const N: usize>(self) -> Self

Computes the algebraic sigmoid function, defined as $\frac{x}{(1 + |x|^N)^{1/N}}$, where N is a positive integer parameter that controls the steepness of the curve.

This also has the unique behavior where for N=0, the function is just the identity function, and for N=1 it is the softsign function.

Note: This function uses $|x|^N$ (the real absolute value), so it is non-holomorphic and only meaningful for real-valued inputs.

To also obtain the derivative with respect to x, use algebraic_sigmoid_d.

fn scalar_algebraic_sigmoid(self, n: u32) -> Self

The algebraic sigmoid for a degree known only at runtime. The runtime twin of algebraic_sigmoid_n, same arithmetic.

fn scalar_algebraic_swish(self) -> Self

Algebraic analogue of the Swish activation, defined as $x\left(\frac{1}{2} + \frac{x}{2\sqrt{1 + x^2}}\right)$. Equivalent to gating x by (1 + algebraic_sigmoid_n::<2>(x)) / 2, the [0, 1]-rescaled N=2 algebraic sigmoid.

Like standard Swish/SiLU, this is smooth and non-monotonic (it dips slightly below zero for moderately negative x before rising) and shares the same asymptotes (f(x) -> x as x -> ∞, f(x) -> 0 as x -> -∞). Unlike Swish, it requires no exp or log, which is substantially cheaper on hardware without fast transcendentals.

To also obtain the derivative with respect to x (which shares most of the underlying computation, notably $1/\sqrt{1 + x^2}$), use algebraic_swish_d.

Historical note

Algebraic gating functions of this form are effectively unknown in modern deep learning, which standardized on exp-based activations (sigmoid, Swish/SiLU, GELU) once GPUs made exp essentially free, a single-cycle special-function-unit op on most modern hardware. On CPUs the calculus is different: a vectorized exp still costs ~20+ cycles even with good polynomial approximations, while sqrt/rsqrt are cheap hardware ops (often approximated in 4-7 cycles). For CPU-side inference, training on CPU, or embedded targets without a transcendental SFU, this remains a competitive Swish-shaped activation at a fraction of the cost.

fn scalar_lgamma_r(self) -> (Self, Self)

Computes the natural log of the Gamma function ($\ln|\Gamma(x)|$) for any real input, for each value in a vector, and returns the sign of the Gamma function from before the absolute value was taken.

fn scalar_gaussian_integral(x0: Self, x1: Self, a: Self, c: Self) -> Self

Computes the definite integral of the Gaussian function from x0 to x1, with amplitude a and standard deviation c. This is more efficient than evaluating the indefinite integral at both limits and subtracting.

The position b is assumed to be zero, so offset the limits accordingly for a non-zero position.

fn scalar_boxcox(self, lambda: Self) -> Self

The Box-Cox transform of x = self with parameter lambda.

\mathrm{boxcox}(x, \lambda) = \begin{cases} \dfrac{x^\lambda - 1}{\lambda} & \lambda \ne 0 \\[6pt] \ln x & \lambda = 0\end{cases}

The variance-stabilizing power transform of applied statistics: $\lambda$ is fitted to make skewed data as close to normal as possible before a model sees it, and the family interpolates the transforms people otherwise pick by hand: $\lambda = 1$ leaves the data alone up to a shift, $1/2$ is a square root, $0$ a logarithm, $-1$ a reciprocal. A fixture of statistical software since Box and Cox introduced it in 1964.

The two cases are one function: $\ln x$ is the limit as $\lambda \to 0$, not a separate rule. Written out, $(x^\lambda - 1)/\lambda$ is $0/0$ there, and the trouble is not confined to the point. Computing $x^\lambda$ and subtracting one cancels, so the naive form is already wrong in the fifth digit at $\lambda = 10^{-12}$ and returns a flat zero by $10^{-300}$. That matters because a fitting routine searches $\lambda$ near zero, which is the usual answer for right-skewed data.

Evaluated as powf_m1(x, lambda)/lambda, which forms $x^\lambda - 1$ without ever forming $x^\lambda$, so there is nothing to cancel and no series or crossover is needed. Measured against a 60-digit oracle, it holds a few ulp from $\lambda = 10^{-300}$ to $\lambda = \pm 8$. Only the exact $\lambda = 0$ is selected apart.

Domain is $x > 0$, and a negative x gives NaN. At $x = 0$ the limits are taken: $-1/\lambda$ for $\lambda > 0$ and $-\infty$ otherwise, which is the conventional choice. That needs no special case: powf_m1(0, lambda) is $-1$ above zero and $+\infty$ below, and the division does the rest.

fn scalar_boxcox_1p(self, lambda: Self) -> Self

The Box-Cox transform of $1 + x$, where x = self.

\mathrm{boxcox1p}(x, \lambda) = \begin{cases} \dfrac{(1 + x)^\lambda - 1}{\lambda} & \lambda \ne 0 \\[6pt] \ln (1 + x) & \lambda = 0\end{cases}

The shifted form exists for the same reason ln_1p does: when x is small, $1 + x$ rounds it away, and every digit of the answer with it. Calling boxcox(1 + x, lambda) loses x entirely once $|x| < \varepsilon$, where this returns $\lambda x$ to full precision. Built on compound_m1, which forms $(1 + x)^\lambda - 1$ without forming either $1 + x$ or $(1+x)^\lambda$.

This is also the kernel underneath yeo_johnson, whose argument is data centered near zero by construction.

Domain is $x > -1$; below that the result is NaN. At $x = -1$ the limits are $-1/\lambda$ for $\lambda > 0$ and $-\infty$ otherwise.

fn scalar_inv_boxcox(self, lambda: Self) -> Self

The inverse Box-Cox transform of y = self with parameter lambda, undoing boxcox.

\mathrm{boxcox}^{-1}(y, \lambda) = \begin{cases} (\lambda y + 1)^{1/\lambda} & \lambda \ne 0 \\[6pt] e^y & \lambda = 0\end{cases}

Wanted by anyone who uses the forward transform: a model fitted on transformed data predicts in transformed units, and the prediction has to come back.

Evaluated as $\exp\!\left(\ln(1 + \lambda y)/\lambda\right)$ rather than as a literal power, which is not merely a rearrangement. The whole point of boxcox is that it stays accurate as $\lambda \to 0$, and $\lambda$ fitted near zero is the common case. There $\lambda y$ is tiny, so forming $\lambda y + 1$ and raising it to the power $1/\lambda$ throws away exactly the digits the forward transform took care to keep. Through ln_1p the exponent tends smoothly to y, so the $\lambda = 0$ case is the limit rather than a discontinuity, and only the exact zero is selected apart.

The range of the forward transform is $\lambda y + 1 > 0$. Outside it the result is NaN, and on the boundary it is $0$ for $\lambda > 0$ and $+\infty$ below.

fn scalar_inv_boxcox_1p(self, lambda: Self) -> Self

The inverse of boxcox_1p.

\mathrm{boxcox1p}^{-1}(y, \lambda) = \begin{cases} (\lambda y + 1)^{1/\lambda} - 1 & \lambda \ne 0 \\[6pt] e^y - 1 & \lambda = 0\end{cases}

The same exponent as inv_boxcox with expm1 outside it instead of exp, so a result near zero keeps its relative accuracy, which, this being the inverse of a transform applied to data centered near zero, is the ordinary case rather than an edge one. Also the kernel underneath inv_yeo_johnson.

fn scalar_yeo_johnson(self, lambda: Self) -> Self

The Yeo-Johnson transform of y = self with parameter lambda.

\psi(y, \lambda) = \begin{cases}
  \dfrac{(y + 1)^\lambda - 1}{\lambda} & y \ge 0,\ \lambda \ne 0 \\[6pt]
  \ln(y + 1) & y \ge 0,\ \lambda = 0 \\[6pt]
  -\dfrac{(1 - y)^{2 - \lambda} - 1}{2 - \lambda} & y < 0,\ \lambda \ne 2 \\[6pt]
  -\ln(1 - y) & y < 0,\ \lambda = 2
\end{cases}

Box-Cox’s sibling, and the one that gets used more, since it is defined on the whole real line rather than on $x > 0$. Same job (fit $\lambda$ by maximum likelihood to make skewed data as close to normal as a power transform can) without the “add a constant to make everything positive first” step, which is an arbitrary choice that changes the fitted $\lambda$. Introduced by Yeo and Johnson in 2000.

One kernel, not four

The four cases are one function seen twice. The $y < 0$ branch is the $y \ge 0$ branch applied to $|y|$ with $\lambda$ reflected to $2 - \lambda$ and the result negated, which is what makes $\psi$ smooth in $\lambda$ across $y = 0$ in the first place. Folding the sign out first therefore collapses the two logarithmic special cases ($\lambda = 0$ above zero, $\lambda = 2$ below) into the single seam that boxcox_1p already handles, and the whole transform is $\pm\,\mathrm{boxcox1p}(|y|, \lambda\ \mathrm{or}\ 2 - \lambda)$.

That the kernel is the 1p form and not boxcox applied to $1 + |y|$ matters here more than anywhere else. $\psi(y, \lambda) \approx y$ near the origin for every $\lambda$, and the origin is where the data is: the transform’s reason for existing is samples that straddle zero. Forming $1 + |y|$ would round away everything below $\varepsilon$ and return a flat zero there.

The value is finite for every finite y, so there is nothing to guard: the two domain edges of the kernel are at $|y| = -1$, which the fold never reaches.

fn scalar_inv_yeo_johnson(self, lambda: Self) -> Self

The inverse Yeo-Johnson transform, undoing yeo_johnson.

\psi^{-1}(z, \lambda) = \begin{cases}
  (\lambda z + 1)^{1/\lambda} - 1 & z \ge 0,\ \lambda \ne 0 \\[6pt]
  e^z - 1 & z \ge 0,\ \lambda = 0 \\[6pt]
  1 - \left((\lambda - 2) z + 1\right)^{1/(2 - \lambda)} & z < 0,\ \lambda \ne 2 \\[6pt]
  1 - e^{-z} & z < 0,\ \lambda = 2
\end{cases}

The same sign fold as the forward transform, over inv_boxcox_1p. $\psi$ is increasing and fixes the origin, so the branch on the way back is the sign of the transformed value, which is the sign of y.

Unlike the forward direction this one has a range to respect: for $\lambda > 0$ the transform’s image is bounded below by $-1/\lambda$, and a z past that came from no y. Such an input gives NaN rather than a plausible-looking number.

fn scalar_spherical_harmonics<const L: usize, const N: usize, const CS: bool>( x: Self, y: Self, z: Self, out: &mut [Self; N], )

Evaluates all real spherical harmonics through degree L at the unit direction (x, y, z), into out[l * (l + 1) + m] for m in -l..=l.

Orthonormal real harmonics. Evaluation is pure polynomial arithmetic: no trigonometry, no division, O(L^2) FMAs total, exact zeros for every m != 0 harmonic at the poles, fully unrolled at compile time for each L up to MAX_SH_DEGREE (above that it takes the rolled general path, which is correct at any degree but roughly 10x slower).

CS picks the phase convention. false gives the standard real-SH tables ($Y_{11} = \sqrt{3/4\pi}\,x$); true applies the Condon-Shortley $(-1)^{|m|}$ phase, negating every odd-|m| harmonic to match Sloan’s SHEval and the physics convention ($Y_{11} = -\sqrt{3/4\pi}\,x$). The choice is baked into a constant table, so neither costs an instruction, but mixing the two silently corrupts any projection/reconstruction round-trip, which is why it must be named.

N must equal (L + 1)^2 (compile-time checked). The direction is assumed unit-length, and nothing renormalizes. See sh_impl for the full convention, algorithm, and domain notes.

use thermite::prelude::*;
use thermite_special::RealSpecialMath;

type V = Vector<f64>;
let (x, y, z) = (V::splat(0.6), V::splat(0.0), V::splat(0.8));

let mut sh = [V::ZERO; 9];
V::spherical_harmonics::<2, 9, false>(x, y, z, &mut sh);
// Y(1,1) = sqrt(3/4pi) * x
assert!((sh[3].extract::<0>() - 0.48860251190292 * 0.6).abs() < 1e-14);

// Condon-Shortley negates odd |m|, and agrees on even |m|.
let mut cs = [V::ZERO; 9];
V::spherical_harmonics::<2, 9, true>(x, y, z, &mut cs);
assert_eq!(cs[3].extract::<0>(), -sh[3].extract::<0>());
assert_eq!(cs[8].extract::<0>(), sh[8].extract::<0>());

fn scalar_spherical_harmonics_table<const L: usize, const N: usize, const CS: bool>( table: &mut ShTable<Self, N>, )

Builds the runtime coefficient table that spherical_harmonics_with and spherical_harmonics_d_with evaluate.

The table depends only on L and CS, never on the direction, so a caller sweeping many directions should build it once rather than calling the one-shot spherical_harmonics per direction. The phase is baked in here, which is why the evaluators take no CS.

The table is typed by Self::Primal, the unaugmented value type: the recurrence coefficients are constants, so a Dual’s derivative parts and a Complex’s imaginary part would only store zeros. For plain vectors and Compensated the primal is Self and nothing changes. For Dual the table is a fraction of the size and its entries multiply as reals.

use thermite::prelude::*;
use thermite_special::{RealSpecialMath, ShTable};

type V = Vector<f64>;
const L: usize = 3;
const N: usize = (L + 1) * (L + 1);

let mut table = ShTable::<V, N>::zeroed();
V::spherical_harmonics_table::<L, N, false>(&mut table);

let mut sh = [V::ZERO; N];
for &(x, y, z) in &[(1.0, 0.0, 0.0), (0.0, 1.0, 0.0)] {
    V::spherical_harmonics_with::<L, N>(
        &table, V::splat(x), V::splat(y), V::splat(z), &mut sh,
    );
}
assert!((sh[1].extract::<0>() - 0.48860251190292).abs() < 1e-14);

fn scalar_spherical_harmonics_with<const L: usize, const N: usize>( table: &ShTable<Self, N>, x: Self, y: Self, z: Self, out: &mut [Self; N], )

Evaluates all harmonics through degree L from a prebuilt table.

The table holds Self::Primal coefficients. See spherical_harmonics_table for how to build it and why, and spherical_harmonics for the conventions and layout.

fn scalar_spherical_harmonics_d<const L: usize, const N: usize, const CS: bool>( x: Self, y: Self, z: Self, out: &mut [Self; N], ddx: &mut [Self; N], ddy: &mut [Self; N], ddz: &mut [Self; N], )

spherical_harmonics plus the ambient Cartesian gradient of every harmonic, into ddx/ddy/ddz.

Lives on RealPrimalMath rather than RealSpecialMath, so Dual does not get it, and should not want it. If you need $\partial/\partial(x,y,z)$, call this directly rather than evaluating spherical_harmonics on a Dual<V, 3> seeded with an identity Jacobian: this shares the recurrence between the value and all three gradients, whereas dual arithmetic carries a derivative through every operation and costs roughly twice as much.

Dual earns its keep on the value form instead, where (x, y, z) are themselves functions of upstream parameters and the chain rule has real work to do. Even there, going the other way (contracting these three gradients against an upstream Jacobian) loses: spherical harmonics cost about two operations per harmonic to evaluate but three per harmonic per parameter to contract, because one recurrence produces the whole basis.

The derivatives are those of the polynomial form at the given (unit) input. Project out the radial component (g - (g . n) n) for the tangential gradient. Shares all recurrence work with the value pass, since the gradients come from tabulated norm ratios, not new recurrences.

fn scalar_zernike_basis_d<const L: usize, const NORM: u8, const N: usize>( x: Self, y: Self, out: &mut [Self; N], ddx: &mut [Self; N], ddy: &mut [Self; N], )

zernike_basis plus $\partial Z_n^m/\partial x$ and $\partial Z_n^m/\partial y$ for every mode, in the same ANSI layout.

This is what a Shack-Hartmann wavefront reconstruction integrates against. The sensor measures local wavefront slopes, not the wavefront itself, so the fit matrix is built from the gradient basis and the value basis never appears in it.

Lives on RealPrimalMath rather than SpecialMath for the same reason spherical_harmonics_d does: Dual should not get it and should not want it. Seeding a Dual<V, 2> and calling the value form carries two derivative components through every operation of the whole ladder, where this differentiates only the two factors that depend on the point and shares the radial recurrence between the value and both gradients.

The gradient is finite everywhere, including the pupil centre. That is the practical dividend of the Cartesian formulation: the polar $\partial_\theta Z/\rho$ is singular there, and hand-rolled polar implementations guard the origin with a special case.

N must equal (L+1)(L+2)/2, and NORM is as on zernike_basis. All three output buffers are written in full.

fn scalar_spherical_harmonics_d_with<const L: usize, const N: usize>( table: &ShTable<Self, N>, x: Self, y: Self, z: Self, out: &mut [Self; N], ddx: &mut [Self; N], ddy: &mut [Self; N], ddz: &mut [Self; N], )

spherical_harmonics_with plus the ambient Cartesian gradients, from a prebuilt table.

fn scalar_softplus_d(self, k: Self, rcp_k: Self) -> (Self, Self)

softplus together with its derivative w.r.t. x (the logistic sigmoid $\sigma(kx)$).

fn scalar_gelu_d(self, alpha: Self) -> (Self, Self)

gelu together with its derivative w.r.t. x.

fn scalar_swish_d(self, beta: Self) -> (Self, Self)

swish together with its derivative w.r.t. x.

fn scalar_algebraic_sigmoid_d_n<const N: usize>(self) -> (Self, Self)

algebraic_sigmoid together with its derivative w.r.t. x.

fn scalar_algebraic_sigmoid_d(self, n: u32) -> (Self, Self)

algebraic_sigmoid_d_n for a degree known only at runtime.

fn scalar_algebraic_swish_d(self) -> (Self, Self)

algebraic_swish together with its derivative w.r.t. x.

fn scalar_langevin_d(self) -> (Self, Self)

langevin together with its derivative $L'(x) = \frac{1}{x^2} - \operatorname{csch}^2 x$.

The derivative shares every intermediate with the value, so this costs a handful of arithmetic ops over langevin alone.

Dyn Compatibility

This trait is not dyn compatible.

In older versions of Rust, dyn compatibility was called "object safety".

Implementors

impl<M> ScalarSpecialMath for M

where M: ScalarSpecialMathWithPolicy,