Trait SpecializedCompensatedSpecialMath 

pub trait SpecializedCompensatedSpecialMath<E>: Sized {
    const INV_LANGEVIN_STEPS: usize;

    // Required method
    fn dd_const(hi: f64, lo: f64) -> Compensated<Self>;

    // Provided methods
    fn compensated_tgamma<P: Policy>(x: Compensated<Self>) -> Compensated<Self>
       where Compensated<Self>: CompensatedGammaOps { ... }
    fn compensated_lgamma_r<P: Policy>(
        x: Compensated<Self>,
    ) -> (Compensated<Self>, Compensated<Self>)
       where Compensated<Self>: CompensatedGammaOps { ... }
    fn compensated_digamma<P: Policy>(x: Compensated<Self>) -> Compensated<Self>
       where Compensated<Self>: CompensatedGammaOps { ... }
    fn compensated_trigamma<P: Policy>(
        x: Compensated<Self>,
    ) -> Compensated<Self>
       where Compensated<Self>: CompensatedGammaOps { ... }
    fn compensated_beta<P: Policy>(
        a: Compensated<Self>,
        b: Compensated<Self>,
    ) -> Compensated<Self>
       where Compensated<Self>: CompensatedGammaOps { ... }
    fn compensated_langevin_d<P: Policy, const ONE_MINUS: bool>(
        x: Compensated<Self>,
    ) -> (Compensated<Self>, Compensated<Self>)
       where Compensated<Self>: CompensatedGammaOps { ... }
    fn compensated_inv_langevin_newton<P: Policy>(
        x: Compensated<Self>,
        y: Compensated<Self>,
        t: Compensated<Self>,
    ) -> Compensated<Self>
       where Compensated<Self>: CompensatedGammaOps { ... }
    fn langevin_small_poly_n(t: Compensated<Self>) -> Compensated<Self>
       where Compensated<Self>: CompensatedGammaOps { ... }
    fn langevin_large_parts<P: Policy>(
        x: Compensated<Self>,
    ) -> (Compensated<Self>, Compensated<Self>, Compensated<Self>)
       where Compensated<Self>: CompensatedGammaOps { ... }
}

Backend for the coefficient-bearing part of Compensated’s special math.

Implemented on the inner vector, not on Compensated

This is implemented for V, with E = V::Element (f32 or f64) as the dispatch tag, and its methods take Compensated<Self> by argument rather than by self. That looks backwards for a math trait and is load-bearing.

The natural spelling - implement it for Compensated<V>, take self, and give it FloatVector<Element = E> as a supertrait so that default bodies can do arithmetic - does not work. Naming that supertrait asserts the projection <Compensated<V> as GenericVector>::Element == Compensated<V::Element> at every use of the bound, and crate::special’s seam then normalizes through that rather than through CompensatedFloatVector, losing the Mask: CastMask<..> obligations its erf / erfinv / lambert_w bodies depend on. It surfaces a hundred lines away as unrelated mismatched types errors in code that was never touched.

Hanging the trait off V avoids that entirely: the seam then constrains V, which cannot say anything about Compensated<V>’s projections. What a default body needs is requested per method via CompensatedGammaOps, scoped to that method alone - which is what makes real default bodies possible here at all.

The element parameter is also what keeps the two per-width impls from colliding: without it both would be impl<V> .. for V, differing only in V::Element, which coherence does not accept as disjoint.

Required Associated Constants

const INV_LANGEVIN_STEPS: usize

Newton steps compensated_inv_langevin_newton needs from the inner vector’s own inv_langevin (~u of that width) to reach this width: one for double-double, two for double-single.

Required Methods

fn dd_const(hi: f64, lo: f64) -> Compensated<Self>

A double-double literal (hi, lo) splat at this width. The Langevin table below is fitted at double-double, so this is a plain splat of both limbs for f64 and a re-split of hi for f32.

Provided Methods

fn compensated_tgamma<P: Policy>(x: Compensated<Self>) -> Compensated<Self>

where Compensated<Self>: CompensatedGammaOps,

$\Gamma(x)$.

Exponentiates compensated_lgamma_r rather than running its own reduction, which costs a few bits and saves a second copy of the reflection: an absolute error d in $\ln\Gamma$ is a relative error d in $\Gamma$, so the loss is $\log_2|\ln\Gamma(x)|$ bits - about 6 near x = 30 and 10 at the overflow edge, out of 106. Avoiding it entirely means a direct Stirling for $\Gamma$, which is only worth writing if those bits are ever missed.

fn compensated_lgamma_r<P: Policy>( x: Compensated<Self>, ) -> (Compensated<Self>, Compensated<Self>)

where Compensated<Self>: CompensatedGammaOps,

$(\ln|\Gamma(x)|, \operatorname{sign}\Gamma(x))$.

The sign is carried separately because lgamma discards it and beta needs it.

Shift-and-Stirling, with the reflection below 1/2. See the module docs for the expansion and for why the shift target is 30.

fn compensated_digamma<P: Policy>(x: Compensated<Self>) -> Compensated<Self>

where Compensated<Self>: CompensatedGammaOps,

$\psi(x)$, the digamma function.

Same shape as compensated_lgamma_r - shift up, then the asymptotic series - but the recurrence $\psi(x) = \psi(x+1) - 1/x$ accumulates a sum of reciprocals rather than a product, so it cannot be deferred to a single log at the end.

fn compensated_trigamma<P: Policy>(x: Compensated<Self>) -> Compensated<Self>

where Compensated<Self>: CompensatedGammaOps,

$\psi_1(x)$, the trigamma function.

As compensated_digamma, with the recurrence $\psi_1(x) = \psi_1(x+1) + 1/x^2$ and the reflection $\psi_1(x) + \psi_1(1-x) = \pi^2/\sin^2(\pi x)$. Note the reflection adds rather than subtracting, unlike digamma’s.

fn compensated_beta<P: Policy>( a: Compensated<Self>, b: Compensated<Self>, ) -> Compensated<Self>

where Compensated<Self>: CompensatedGammaOps,

$B(a, b) = \Gamma(a)\Gamma(b)/\Gamma(a+b)$.

Through logs rather than as a ratio of gammas, which overflows for arguments the beta function itself handles perfectly well. Rides entirely on compensated_lgamma_r, so a width that overrides that one gets this for free and should never need to touch this.

fn compensated_langevin_d<P: Policy, const ONE_MINUS: bool>( x: Compensated<Self>, ) -> (Compensated<Self>, Compensated<Self>)

where Compensated<Self>: CompensatedGammaOps,

L(x) (or 1 - L(x) with ONE_MINUS) and L'(x). Same structure as thermite-special’s kernel with the crossover at |x| = 1 (3u/x^2 of cancellation is 3 ulp there, and the double-double table is 22 terms already). The complement on the large branch is 1/x - 2q/(1-q), which at worst (x = 1) cancels to 0.69 of 1/x.

fn compensated_inv_langevin_newton<P: Policy>( x: Compensated<Self>, y: Compensated<Self>, t: Compensated<Self>, ) -> Compensated<Self>

where Compensated<Self>: CompensatedGammaOps,

One Newton step of L(x) = y from x, with t = 1 - y supplied exactly (the residual is ((1-y) - 1/x) + 2q/(1-q) on the large branch, which is what keeps the step accurate where L sits within an ulp of 1).

fn langevin_small_poly_n(t: Compensated<Self>) -> Compensated<Self>

where Compensated<Self>: CompensatedGammaOps,

Minimax fit of L(x)/x in x^2 on [0, 1] at double-double (relative error 3e-36, crates/thermite-special/scripts/langevin_coeffs_dd.py), Horner.

fn langevin_large_parts<P: Policy>( x: Compensated<Self>, ) -> (Compensated<Self>, Compensated<Self>, Compensated<Self>)

where Compensated<Self>: CompensatedGammaOps,

1/x, 2q/(1-q) and csch^2(x) for x >= 1, one division between them (r = 1/(x(1-q)), 1/x = (1-q) r, 1/(1-q) = x r). At x = inf the products are inf * 0, so those lanes are set to their limits explicitly.

Dyn Compatibility

This trait is not dyn compatible.

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

Implementors

impl<V: FloatVector<Element = f32>> SpecializedCompensatedSpecialMath<f32> for V

const INV_LANGEVIN_STEPS: usize = 2

impl<V: FloatVector<Element = f64>> SpecializedCompensatedSpecialMath<f64> for V

const INV_LANGEVIN_STEPS: usize = 1