Trait RealSpecialMath 

pub trait RealSpecialMath: RealSpecialMathWithPolicy {
    // Provided methods
    fn erfinv(self) -> Self { ... }
    fn probit(self) -> Self { ... }
    fn ndtr(self) -> Self { ... }
    fn log_ndtr(self) -> Self { ... }
    fn logerfc(self) -> Self { ... }
    fn fresnel(self) -> (Self, Self) { ... }
    fn fresnel_c(self) -> Self { ... }
    fn fresnel_s(self) -> Self { ... }
    fn sici(self) -> (Self, Self) { ... }
    fn sinint(self) -> Self { ... }
    fn cosint(self) -> Self { ... }
    fn inv_log_ndtr(self) -> Self { ... }
    fn inv_digamma(self) -> Self { ... }
    fn wright_omega(self) -> Self { ... }
    fn bessel_ratio<F: BesselRatioFamily>(self, nu: Self) -> Self { ... }
    fn inv_bessel_ratio<F: BesselRatioFamily>(self, nu: Self) -> Self { ... }
    fn bessel_ratio_1m<F: BesselRatioFamily>(self, nu: Self) -> Self { ... }
    fn inv_bessel_ratio_1m<F: BesselRatioFamily>(self, nu: Self) -> Self { ... }
    fn gauss_legendre(self, n: u32) -> (Self, Self) { ... }
    fn gauss_hermite(self, n: u32) -> (Self, Self) { ... }
    fn gauss_laguerre(self, alpha: Self, n: u32) -> (Self, Self) { ... }
    fn pochhammer(self, m: Self) -> Self { ... }
    fn jacobi_elliptic(self, k: Self) -> (Self, Self, Self) { ... }
    fn agm(self, other: Self) -> Self { ... }
    fn langevin(self) -> Self { ... }
    fn inv_langevin(self) -> Self { ... }
    fn langevin_1m(self) -> Self { ... }
    fn inv_langevin_1m(self) -> Self { ... }
    fn gelu(self, alpha: Self) -> Self { ... }
    fn swish(self, beta: Self) -> Self { ... }
    fn algebraic_sigmoid_n<const N: usize>(self) -> Self { ... }
    fn algebraic_sigmoid(self, n: u32) -> Self { ... }
    fn algebraic_swish(self) -> Self { ... }
    fn lgamma_r(self) -> (Self, Self) { ... }
    fn gaussian_integral(x0: Self, x1: Self, a: Self, c: Self) -> Self { ... }
    fn boxcox(self, lambda: Self) -> Self { ... }
    fn boxcox_1p(self, lambda: Self) -> Self { ... }
    fn inv_boxcox(self, lambda: Self) -> Self { ... }
    fn inv_boxcox_1p(self, lambda: Self) -> Self { ... }
    fn yeo_johnson(self, lambda: Self) -> Self { ... }
    fn inv_yeo_johnson(self, lambda: Self) -> Self { ... }
    fn spherical_harmonics<const L: usize, const N: usize, const CS: bool>(
        x: Self,
        y: Self,
        z: Self,
        out: &mut [Self; N],
    ) { ... }
    fn spherical_harmonics_table<const L: usize, const N: usize, const CS: bool>(
        table: &mut ShTable<<Self as PrimalProjection>::Primal, N>,
    ) { ... }
    fn spherical_harmonics_with<const L: usize, const N: usize>(
        table: &ShTable<<Self as PrimalProjection>::Primal, N>,
        x: Self,
        y: Self,
        z: Self,
        out: &mut [Self; N],
    ) { ... }
}

RealSpecial math functions for floating-point vectors using the default policy.

The same operations as RealSpecialMathWithPolicy, with every method’s leading policy fixed to DefaultPolicy. Implementors of RealSpecialMathWithPolicy automatically implement this trait; each method here has a _p-suffixed counterpart there. Special math functions that are only defined for real-valued floating-point vectors.

These functions either rely on ordering/sign information that has no complex analogue (e.g. erfinv, probit, lgamma_r), or use the real absolute value in a way that makes them non-holomorphic (e.g. algebraic_sigmoid).

Provided Methods

fn erfinv(self) -> Self

Computes the inverse error function.

fn probit(self) -> Self

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

fn 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 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 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 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 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 fresnel_s(self) -> Self

$S(x)$ alone. See fresnel_c.

fn 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 sinint(self) -> Self

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

fn cosint(self) -> Self

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

fn 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 spherical_harmonics_table<const L: usize, const N: usize, const CS: bool>( table: &mut ShTable<<Self as PrimalProjection>::Primal, 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 spherical_harmonics_with<const L: usize, const N: usize>( table: &ShTable<<Self as PrimalProjection>::Primal, 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.

Dyn Compatibility

This trait is not dyn compatible.

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

Implementors

impl<M> RealSpecialMath for M

where M: RealSpecialMathWithPolicy,