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. 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)) form.

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.

The output is symmetric in (a, b): t3 is fully index-symmetric, so t3[a][b][c] == t3[b][a][c] and the Σ_c contraction gives out[a][b] == out[b][a] term-for-term, in the same c order. We compute only the upper triangle a ≤ b (the inner contraction is unchanged and stays contiguous/vectorisable) and mirror into the lower triangle — this is BIT-IDENTICAL to the full a, b ∈ 0..K nest while doing ~2× fewer inner contractions, with no dense scatter (the mirror is a K × K copy).

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.

As in Self::third_contracted, the output is symmetric in (i, j) (t4[j][i][k][l] == t4[i][j][k][l], contracted in the same (k, l) order), so the upper triangle i ≤ j is computed and mirrored — BIT-IDENTICAL to the full nest, ~2× fewer inner Σ_{k,l} contractions, and the inner double loop stays the original contiguous/vectorisable form.

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> JetField for Tower4<K>

fn value(&self) -> f64

The real value channel (recurses through any nesting to the f64 leaf).

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

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

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

fn neg(&self) -> Self

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

Multiply every channel by a plain f64.

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

Faà di Bruno composition f ∘ self given the OUTER real function’s derivative stack d = [f(u), f′(u), f″(u), f‴(u), f⁗(u)] evaluated at u = self.value() — the identical [f64; 5] stack shape crate::jet_tower::Tower4::compose_unary consumes.

fn constant_like(&self, v: f64) -> Self

A constant carrying THIS element’s shape: real value v, every derivative channel zero.

fn with_value(&self, v: f64) -> Self

self with its real value channel replaced by v, every derivative channel untouched.

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 symmetric_quadratic_form<C: SymmetricQuadraticCoefficients>( inputs: &[Self], coefficients: &C, ) -> Self

Evaluate inputs' A inputs from one universal semantic primitive. Order-specific scalars may lower the mechanically derived channels directly; the default is the exact scalar program over mul/add/scale.

fn linear_combination(inputs: &[Self], weights: &[f64]) -> Self

Evaluate sum_i weights[i] * inputs[i] in one semantic primitive.

fn add_constant(&self, constant: f64) -> Self

Add a primal constant without changing derivative channels.

fn multiply_add(&self, right: &Self, addend: &Self) -> Self

Evaluate self * right + addend in one semantic primitive.

fn composed_sum(inputs: &[Self], derivative_stacks: &[[f64; 5]]) -> Self

Sum unary compositions directly from certified derivative stacks.

fn product(&self, right: &Self) -> Self

Exact product as an explicit compiled graph node.

fn affine_compose( &self, input_scale: f64, input_shift: f64, derivative_stack: [f64; 5], ) -> Self

Compose a certified outer stack after the affine map u = input_scale * self + input_shift.

fn affine_composed_sum( inputs: &[Self], input_scales: &[f64], derivative_stacks: &[[f64; 5]], ) -> Self

Sum unary compositions whose inputs each carry an affine scale.

fn shared_multiply_add_affine_composed_sum<const N: usize>( lefts: &[&Self; N], right: &Self, addend: &Self, addend_scales: &[f64; N], input_scales: &[f64; N], derivative_stacks: &[[f64; 5]; N], ) -> Self

Evaluate Σ_i f_i(input_scale_i · (left_i · right + addend_scale_i · addend)) from the certified derivative stack of each f_i. The shared operands make expression-level common subexpressions explicit, so optimized backends apply their inherited derivative channels once. Expression arity is part of the type, and borrowed operands never copy a full tower. An exact-zero addend scale (either sign of IEEE zero) removes that addend from the corresponding term entirely.

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. This evaluates stack_fn(self.value()) once and forwards to compose_unary, so it is bit-identical to the explicit self.compose_unary(stack_fn(self.value())) form.

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.