Trait SpecialMathWithPolicy 

pub trait SpecialMathWithPolicy: TranscendentalMathWithPolicy {
    // Required methods
    fn erf_p<P: Policy>(self) -> Self;
    fn erfc_p<P: Policy>(self) -> Self;
    fn erfcx_p<P: Policy>(self) -> Self;
    fn logistic_sigmoid_p<P: Policy>(self) -> Self;
    fn logit_p<P: Policy>(self) -> Self;
    fn logit_1m_p<P: Policy>(self) -> Self;
    fn softplus_p<P: Policy>(self, k: Self, rcp_k: Self) -> Self;
    fn tgamma_p<P: Policy>(self) -> Self;
    fn lgamma_p<P: Policy>(self) -> Self;
    fn poisson_pmf_p<P: Policy>(self, lambda: Self) -> Self;
    fn poisson_log_pmf_p<P: Policy>(self, lambda: Self) -> Self;
    fn digamma_p<P: Policy>(self) -> Self;
    fn trigamma_p<P: Policy>(self) -> Self;
    fn polygamma_p<P: Policy>(self, n: u32) -> Self;
    fn zeta_p<P: Policy>(self) -> Self;
    fn zetac_p<P: Policy>(self) -> Self;
    fn polylog_p<P: Policy>(
        self,
        order: PolylogOrder<Self::Element, <Self::Signed as GenericVector>::Element>,
    ) -> Self;
    fn bessel_n_p<P: Policy, F: BesselFamily, const N: i32>(self) -> Self;
    fn bessel_p<P: Policy, F: BesselFamily>(
        self,
        order: BesselOrder<Self, Self::Signed>,
    ) -> Self;
    fn sph_bessel_n_p<P: Policy, F: BesselFamily, const N: usize>(self) -> Self;
    fn sph_bessel_p<P: Policy, F: BesselFamily>(self, n: u32) -> Self;
    fn airy_p<P: Policy, W: AiryFn>(self) -> Self;
    fn airy_all_p<P: Policy, const SCALED: bool>(
        self,
    ) -> (Self, Self, Self, Self);
    fn beta_p<P: Policy>(self, y: Self) -> Self;
    fn lbeta_p<P: Policy>(self, y: Self) -> Self;
    fn jacobi_p<P: Policy>(
        self,
        alpha: Self,
        beta: Self,
        n: u32,
        m: u32,
    ) -> Self;
    fn hermite_n_p<P: Policy, const N: usize>(self) -> Self;
    fn hermitev_p<P: Policy>(self, n: Self::Unsigned) -> Self;
    fn hermite_p<P: Policy>(self, n: u32) -> Self;
    fn hermite_function_n_p<P: Policy, const N: usize>(self) -> Self;
    fn hermite_function_p<P: Policy>(self, n: u32) -> Self;
    fn hermite_function_series_n_p<P: Policy, const N: usize>(
        self,
        coeffs: &[Self::Element; N],
    ) -> Self;
    fn hermite_function_series_p<P: Policy>(
        self,
        coeffs: &[Self::Element],
    ) -> Self;
    fn laguerre_n_p<P: Policy, const N: usize>(self, alpha: Self) -> Self;
    fn laguerrev_p<P: Policy>(self, alpha: Self, n: Self::Unsigned) -> Self;
    fn laguerre_p<P: Policy>(self, alpha: Self, n: u32) -> Self;
    fn laguerre_function_n_p<P: Policy, const N: usize>(
        self,
        alpha: Self,
    ) -> Self;
    fn laguerre_function_p<P: Policy>(self, alpha: Self, n: u32) -> Self;
    fn laguerre_function_i_n_p<P: Policy, const N: usize>(
        self,
        alpha: i32,
    ) -> Self;
    fn laguerre_function_i_p<P: Policy>(self, alpha: i32, n: u32) -> Self;
    fn laguerre_function_series_n_p<P: Policy, const N: usize>(
        self,
        alpha: Self,
        coeffs: &[Self::Element; N],
    ) -> Self;
    fn laguerre_function_series_p<P: Policy>(
        self,
        alpha: Self,
        coeffs: &[Self::Element],
    ) -> Self;
    fn laguerre_function_series_i_n_p<P: Policy, const N: usize>(
        self,
        alpha: i32,
        coeffs: &[Self::Element; N],
    ) -> Self;
    fn laguerre_function_series_i_p<P: Policy>(
        self,
        alpha: i32,
        coeffs: &[Self::Element],
    ) -> Self;
    fn chebyshev_n_p<P: Policy, const K: usize, const N: usize>(
        self,
        coeffs: &[Self::Element; N],
    ) -> Self;
    fn chebyshev_p<P: Policy, const K: usize>(
        self,
        coeffs: &[Self::Element],
    ) -> Self;
    fn gaussian_p<P: Policy>(self, a: Self, c: Self) -> Self;
    fn planck_p<P: Policy>(self) -> Self;
    fn legendre_p<P: Policy>(self, n: u32, m: u32) -> Self;
    fn legendre_series_n_p<P: Policy, const N: usize>(
        self,
        coeffs: &[Self::Element; N],
    ) -> Self;
    fn legendre_series_p<P: Policy>(self, coeffs: &[Self::Element]) -> Self;
    fn zernike_r_p<P: Policy>(self, n: u32, m: u32) -> Self;
    fn zernike_p<P: Policy, const NORM: u8>(
        self,
        theta: Self,
        n: u32,
        m: i32,
    ) -> Self;
    fn zernike_basis_p<P: Policy, const L: usize, const NORM: u8, const N: usize>(
        x: Self,
        y: Self,
        out: &mut [Self; N],
    );
    fn lambert_w_p<P: Policy>(self) -> (Self, Self);
    fn expint_n_p<P: Policy, const N: usize>(self) -> Self;
    fn expint_p<P: Policy>(self, n: u32) -> Self;
    fn phi_n_p<P: Policy, const N: usize>(self) -> Self;
    fn phi_p<P: Policy>(self, n: u32) -> Self;
    fn carlson_p<P: Policy, K: CarlsonKind<Output = Self>>(kind: K) -> Self;
    fn ellint_p<P: Policy, K: EllipticKind<Output = Self>>(kind: K) -> Self;
}

Special math functions for floating-point vectors with customizable policies.

Each function carries a leading P: Policy generic controlling the precision/performance trade-off. For convenience, SpecialMath provides the same set of operations under DefaultPolicy. Every floating-point vector type implementing specialized::SpecializedSpecialMath implements both automatically. Special math functions that are valid for both real and complex floating-point vectors.

Required Methods

fn erf_p<P: Policy>(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 erfc_p<P: Policy>(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 erfcx_p<P: Policy>(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 logistic_sigmoid_p<P: Policy>(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 logit_p<P: Policy>(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 logit_1m_p<P: Policy>(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 softplus_p<P: Policy>(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 tgamma_p<P: Policy>(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 lgamma_p<P: Policy>(self) -> Self

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

fn poisson_pmf_p<P: Policy>(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 poisson_log_pmf_p<P: Policy>(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 digamma_p<P: Policy>(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 trigamma_p<P: Policy>(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 polygamma_p<P: Policy>(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 zeta_p<P: Policy>(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 zetac_p<P: Policy>(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 polylog_p<P: Policy>( self, order: PolylogOrder<Self::Element, <Self::Signed as GenericVector>::Element>, ) -> 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 bessel_n_p<P: Policy, 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 bessel_p<P: Policy, 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 sph_bessel_n_p<P: Policy, 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 sph_bessel_p<P: Policy, F: BesselFamily>(self, n: u32) -> Self

sph_bessel_n for an order known only at runtime.

fn airy_p<P: Policy, 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 airy_all_p<P: Policy, 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 beta_p<P: Policy>(self, y: Self) -> Self

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

fn lbeta_p<P: Policy>(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 jacobi_p<P: Policy>(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 hermite_n_p<P: Policy, 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 hermitev_p<P: Policy>(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 hermite_p<P: Policy>(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 hermite_function_n_p<P: Policy, 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 hermite_function_p<P: Policy>(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 hermite_function_series_n_p<P: Policy, 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 hermite_function_series_p<P: Policy>(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 laguerre_n_p<P: Policy, 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 laguerrev_p<P: Policy>(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 laguerre_p<P: Policy>(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 laguerre_function_n_p<P: Policy, 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 laguerre_function_p<P: Policy>(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 laguerre_function_i_n_p<P: Policy, 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 laguerre_function_i_p<P: Policy>(self, alpha: i32, n: u32) -> Self

laguerre_function_i_n for a degree known only at runtime.

fn laguerre_function_series_n_p<P: Policy, 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 laguerre_function_series_p<P: Policy>( 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 laguerre_function_series_i_n_p<P: Policy, 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 laguerre_function_series_i_p<P: Policy>( 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 chebyshev_n_p<P: Policy, 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 chebyshev_p<P: Policy, 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 gaussian_p<P: Policy>(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 planck_p<P: Policy>(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 legendre_p<P: Policy>(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 legendre_series_n_p<P: Policy, 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 legendre_series_p<P: Policy>(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 zernike_r_p<P: Policy>(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 zernike_p<P: Policy, 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 zernike_basis_p<P: Policy, 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 lambert_w_p<P: Policy>(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 expint_n_p<P: Policy, const N: usize>(self) -> Self

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

fn expint_p<P: Policy>(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 phi_n_p<P: Policy, 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 phi_p<P: Policy>(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 carlson_p<P: Policy, K: CarlsonKind<Output = Self>>(kind: K) -> 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 ellint_p<P: Policy, K: EllipticKind<Output = Self>>(kind: K) -> 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)

Dyn Compatibility

This trait is not dyn compatible.

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

Implementors

impl<E, V> SpecialMathWithPolicy for V

where V: SpecializedSpecialMath<E> + FloatVector<Element = E> + TranscendentalMathWithPolicy,