Struct Tower4 

pub struct Tower4<const K: usize> {
    pub v: f64,
    pub g: [f64; K],
    pub h: [[f64; K]; K],
    pub t3: [[[f64; K]; K]; K],
    pub t4: [[[[f64; K]; K]; K]; K],
}

Truncated fourth-order multivariate Taylor scalar in K variables.

See the module documentation for semantics and conventions. Copy is intentional despite the size (2 KiB at K=4): towers are per-row temporaries that live entirely in registers/stack during a row program, and value semantics keep program code readable (a * b + c).

Fields

v: f64

Value ℓ.

g: [f64; K]

Gradient ∂ℓ/∂p_a.

h: [[f64; K]; K]

Hessian ∂²ℓ/∂p_a∂p_b (symmetric).

t3: [[[f64; K]; K]; K]

Third derivatives ∂³ℓ/∂p_a∂p_b∂p_c (fully symmetric).

t4: [[[[f64; K]; K]; K]; K]

Fourth derivatives ∂⁴ℓ/∂p_a∂p_b∂p_c∂p_d (fully symmetric).

Implementations

impl<const K: usize> Tower4<K>

pub fn zero() -> Self

The additive identity.

pub fn constant(c: f64) -> Self

A constant: value c, all derivatives zero.

pub fn variable(value: f64, idx: usize) -> Self

The seeded variable p_idx with current value value: unit first derivative in slot idx, zero elsewhere and above.

pub fn mul(&self, o: &Self) -> Self

Exact truncated Leibniz product D_S(ab) = Σ_{T ⊆ S} D_T(a) · D_{S∖T}(b).

Codegen

Each output entry’s 2^m subset sum is written as a compact straight-line expression instead of the shared [jet_algebra::leibniz_product] subset walker (which, per entry, builds SlotBufs and match-dispatches the deriv closure across all 2^m subsets). The loop nest over (i,j,k,l) is unchanged — only the inner per-entry sum is unrolled — so this does NOT unroll over K and does NOT bloat code: on a Tower4<9> mul-and-read consumer the new form is faster AND smaller (asm: 34 outlined walker bl calls → 0, 21.1 KiB → 14.3 KiB, +100 NEON .2d ops).

BIT-IDENTICAL to the walker: each entry’s terms are in the walker’s exact subset-enumeration order (subset bit b ↔ position b, sub = 0..2^m), and the per-entry acc accumulator mirrors the walker’s total = 0.0 start so a signed-zero leading product collapses to +0.0 identically — which matters because real jets carry exact-0.0 channels (constant/variable towers). Proven to_bits-identical on v/g/h/t3/t4 across K ∈ {2,3,4,9}, 5000 inputs each with ~30 % exact-0.0 channels and signed values (a no-leading-0.0 form fails this stress — the accumulator start is load-bearing).

pub fn add(&self, o: &Self) -> Self

Ref-taking elementwise sum, the by-ref twin of the std::ops::Add operator (which consumes by value). Mirrors the inherent mul/scale API so a chain like a.mul(&b).add(&c) reads uniformly without moving out of the borrowed operands.

pub fn sub(&self, o: &Self) -> Self

Ref-taking elementwise difference, the by-ref twin of std::ops::Sub.

pub fn compose_unary(&self, d: [f64; 5]) -> Self

Exact multivariate Faà di Bruno composition f ∘ self.

d = [f(u), f′(u), f″(u), f‴(u), f⁗(u)] evaluated at u = self.v — the SAME [f64; 5] stack shape the families’ existing unary_derivatives_* helpers produce, so those special-function stacks (Φ, log-Φ, normal pdf, …) plug in directly.

The order-m output sums over the set partitions of the m indices (Bell(3) = 5 terms at order 3, Bell(4) = 15 at order 4), grouped by block count: each partition into r blocks contributes f⁽ʳ⁾ · Π_blocks D_block(u).

Codegen

Evaluated as a compact closed form (the Bell(4)=15 set-partitions of t4, Bell(3)=5 of t3, …) instead of routing through the recursive jet_algebra::faa_di_bruno walker (per-output for_each_partition recursion + per-block SlotBuf + closure dispatch). The loop nest is identical to the walker’s (for i,j,k,l); only the per-entry partition sum is straight-line, so this does NOT unroll over K and does NOT bloat code — measured on a Tower4<9> compose-and-read consumer the new form is both faster and SMALLER (asm: 94 outlined walker bl calls → 0, 47.5 KiB → 16.7 KiB, +197 NEON .2d ops).

BIT-IDENTICAL to the walker: each channel’s terms are emitted in the walker’s exact partition-enumeration order, each term’s block products are left-associated exactly as the walker’s prod *= block, and the per-channel acc accumulator mirrors the walker’s total = 0.0 start (so signed-zero products collapse to +0.0 identically). The order-4 term sequence was generated from the walker’s own enumeration. Proven to_bits-identical on v/g/h/t3/t4 across K ∈ {2,3,4,9}, 5000 random inputs each (zeroed / sign-varied stacks included).

pub fn compose_unary_with(&self, stack_fn: impl Fn(f64) -> [f64; 5]) -> Self

Compose with a unary special-function whose [f64; 5] derivative STACK is built from the base value through stack_fn — the scalar arm of the generic-over-Lane compose seam (see Tower4Lane::compose_unary_with). Evaluates stack_fn(self.v) ONCE and forwards to Self::compose_unary, so it is BIT-IDENTICAL to the explicit self.compose_unary(stack_fn(self.v)). Writing a program against this seam lets it re-instantiate, unchanged, at Tower4Lane (where each of the four lanes carries a DISTINCT base value and stack_fn is re-run per lane).

pub fn compose_unary_single_slot(&self, d: [f64; 5], slot: usize) -> Self

Single-active-slot fast path for Self::compose_unary.

When the inner jet self has derivative support ONLY on the all-slot diagonal channels — i.e. it is a univariate jet in primary slot scattered into the K-wide layout (g[a] = 0, h[a][b] = 0, t3 = 0, t4 = 0 for any axis ≠ slot) — the multivariate Faà di Bruno walk collapses. Every output channel whose axis tuple contains an axis ≠ slot is structurally 0: each set-partition has a block covering that axis, that block reads an off-slot derivative of self (which is 0), so the block product and the whole partition vanish, and the channel sums to the walker’s total = 0.0 start, i.e. +0.0. Only the five diagonal channels (v, g[slot], h[slot][slot], t3[slot]³, t4[slot]⁴) survive.

This computes exactly those five as STRAIGHT-LINE accumulations, each in the EXACT term order of Self::compose_unary’s diagonal (i = j = k = l = slot) case — so they are BIT-IDENTICAL to Self::compose_unary on the diagonal — and leaves every other channel at the zero-init +0.0, which the full walk also produces (the off-slot collapse is to_bits-+0.0, signed-zero products included; proven across K ∈ {2,3,4,9}, 5000 single-slot inputs each). At any K ≥ 2 this is far fewer floating-point operations than materialising the full 1 + K + K² + K³ + K⁴ channel set whose off-diagonal entries are all zero, and far cheaper than the recursive set-partition walker the diagonal channels previously routed through (a measured ~9.5× speedup vs the full compose_unary, recovering a 5.9× walker regression at the K ∈ {2,3} BMS tower widths).

#[inline] so an adopting consumer pays no bl call (uninlined, the five-channel build does not amortise the call/spill overhead).

Precondition

The caller guarantees the single-active-slot structure. If it does not hold, the off-slot channels would be wrongly zeroed; use the full Self::compose_unary in that case.

pub fn scale(&self, s: f64) -> Self

Multiply every channel by a plain scalar.

pub fn exp(&self) -> Self

e^self.

pub fn ln(&self) -> Self

ln(self). Caller guarantees positivity (likelihood programs do).

pub fn recip(&self) -> Self

1/self.

pub fn sqrt(&self) -> Self

√self. Caller guarantees positivity.

pub fn powf(&self, a: f64) -> Self

self^a for real exponent a. Caller guarantees a positive base.

pub fn ln_gamma(&self) -> Self

ln Γ(self). Caller guarantees positivity.

pub fn digamma(&self) -> Self

ψ(self), the digamma function. Caller guarantees positivity.

pub fn trigamma(&self) -> Self

ψ′(self), the trigamma function. Caller guarantees positivity.

pub fn third_contracted(&self, dir: &[f64; K]) -> [[f64; K]; K]

Contract t3 with one primary-space direction: out[a][b] = Σ_c t3[a][b][c] · dir[c] — exactly the row_third_contracted shape.

pub fn fourth_contracted(&self, u: &[f64; K], w: &[f64; K]) -> [[f64; K]; K]

Contract t4 with two primary-space directions: out[a][b] = Σ_{c,d} t4[a][b][c][d] · u[c] · v[d] — exactly the row_fourth_contracted shape.

Trait Implementations

impl<const K: usize> Add for Tower4<K>

type Output = Tower4<K>

The resulting type after applying the + operator.

fn add(self, o: Self) -> Self

Performs the + operation.

impl<const K: usize> Add<f64> for Tower4<K>

type Output = Tower4<K>

The resulting type after applying the + operator.

fn add(self, c: f64) -> Self

Performs the + operation.

impl<const K: usize> Clone for Tower4<K>

fn clone(&self) -> Tower4<K>

Returns a duplicate of the value.

fn clone_from(&mut self, source: &Self)

Performs copy-assignment from source.

impl<const K: usize> Copy for Tower4<K>

impl<const K: usize> Debug for Tower4<K>

fn fmt(&self, f: &mut Formatter<'_>) -> Result

Formats the value using the given formatter.

impl<const K: usize> Div for Tower4<K>

type Output = Tower4<K>

The resulting type after applying the / operator.

fn div(self, o: Self) -> Self

Performs the / operation.

impl<const K: usize> JetScalar<K> for Tower4<K>

The full dense crate::jet_tower::Tower4 is itself a JetScalar: it carries EVERY channel, so a row expression written ONCE against JetScalar can be evaluated at Tower4 to obtain the full (v, g, H, t3, t4) in one pass. This is BOTH the #932 oracle ground truth the packed Order2 / OneSeed / TwoSeed scalars are pinned against, AND a production scalar: a family whose uncontracted third / fourth derivative tensors are needed (the BMS rigid third_full / fourth_full caches) evaluates the SAME generic row-NLL expression at Tower4 and reads .t3 / .t4 off the result — so the dense tensors come from the single source of truth, not a separately hand-written jet. The packed scalars serve the consumers that need only (v, g, H) (Order2) or one / two contractions (OneSeed / TwoSeed) without paying for the dense tensors.

fn constant(c: f64) -> Self

A constant: value c, every derivative channel zero.

fn variable(x: f64, axis: usize) -> Self

The seeded variable p_axis at value x: unit first derivative in slot axis, all higher channels zero. (The nilpotent / cross channels of the directional scalars are seeded zero — callers set ε/δ directions through the scalar-specific OneSeed::seed_direction / TwoSeed::seed.)

fn value(&self) -> f64

The value channel ℓ(p).

fn add(&self, o: &Self) -> Self

Exact truncated Leibniz sum self + o.

fn sub(&self, o: &Self) -> Self

Exact truncated Leibniz difference self − o.

fn mul(&self, o: &Self) -> Self

Exact truncated Leibniz product self · o.

fn neg(&self) -> Self

Negate every channel.

fn scale(&self, s: f64) -> Self

Multiply every channel by a plain scalar s.

fn compose_unary(&self, d: [f64; 5]) -> Self

Exact multivariate Faà di Bruno composition f ∘ self, given the outer derivative stack d = [f(u), f′(u), f″(u), f‴(u), f⁗(u)] at u = self.value().

fn compose_unary_with(&self, stack_fn: impl Fn(f64) -> [f64; 5]) -> Self

Compose with a unary special-function whose derivative STACK is built from the scalar base value through stack_fn — the generic-over-Lane seam that lets a single-sourced row program instantiate at BOTH the scalar f64 jets and the SIMD f64x4 batch towers from ONE expression.

fn exp(&self) -> Self

e^self. Convenience for tame arguments (see module stability note).

fn sqrt(&self) -> Self

√self. Caller guarantees positivity.

fn ln(&self) -> Self

ln(self). Caller guarantees positivity. Same derivative stack crate::jet_tower::Tower4::ln uses, so any program written over both matches term-for-term.

fn recip(&self) -> Self

1/self.

fn powf(&self, a: f64) -> Self

self^a for real exponent a. Caller guarantees a positive base. Mirrors crate::jet_tower::Tower4::powf (falling-factorial stack).

fn ln_gamma(&self) -> Self

ln Γ(self). Caller guarantees a positive argument. Uses the SAME hand-certified derivative stack crate::jet_tower::Tower4::ln_gamma consumes (crate::jet_tower::ln_gamma_derivative_stack), so any program written over both matches term-for-term.

fn digamma(&self) -> Self

ψ(self) = d/dx ln Γ(x) (digamma). Caller guarantees a positive argument. Same hand-certified stack crate::jet_tower::digamma_derivative_stack.

impl<const K: usize> Mul for Tower4<K>

type Output = Tower4<K>

The resulting type after applying the * operator.

fn mul(self, o: Self) -> Self

Performs the * operation.

impl<const K: usize> Mul<f64> for Tower4<K>

type Output = Tower4<K>

The resulting type after applying the * operator.

fn mul(self, c: f64) -> Self

Performs the * operation.

impl<const K: usize> Neg for Tower4<K>

type Output = Tower4<K>

The resulting type after applying the - operator.

fn neg(self) -> Self

Performs the unary - operation.

impl<const K: usize> Sub for Tower4<K>

type Output = Tower4<K>

The resulting type after applying the - operator.

fn sub(self, o: Self) -> Self

Performs the - operation.

impl<const K: usize> Sub<f64> for Tower4<K>

type Output = Tower4<K>

The resulting type after applying the - operator.

fn sub(self, c: f64) -> Self

Performs the - operation.