pub fn implicit_solve<const K1: usize, const K: usize>(
f: &Tower4<K1>,
a0: f64,
) -> Result<Tower4<K>, String>Expand description
Solve the implicit relation F(a(θ), θ) ≡ 0 for the intercept tower
a(θ) over the K primaries θ, given the constraint tower f written
over K + 1 variables (slot 0 is the intercept a, slots 1..=K
are the primaries θ) evaluated at the SOLVED point — i.e. f.v is the
constraint residual at (a₀, θ₀) (≈ 0 from the production Newton solve)
and a0 is that solved intercept value.
Returns the Tower4<K> whose value is a0 and whose every derivative
tensor (∂a/∂θ, ∂²a/∂θ², …, ∂⁴a/∂θ⁴) is the exact implicit-function
derivative. This is the mechanical replacement for the hand-coded
a_u = -f_u/f_a, a_uv = -(f_uv + f_au·a_v + f_av·a_u + f_aa·a_u·a_v)/f_a
recursion (first_full.rs) and its third/fourth-order continuations.
Method: order-by-order substitution. We build a incrementally; at each
order m the composite G(θ) = f(a(θ), θ) has a top-order coefficient
that is linear in a’s order-m tensor with leading factor F_a
(= f.g[0]), plus terms in a’s lower orders already fixed. Setting the
order-m tensor of a to cancel the rest of G’s order-m coefficient
keeps G ≡ 0 through that order. The substitution G = f∘(a, θ) reuses
only the exact substitute_intercept chain rule, so the recursion is
auditable and exact, not a hand-expanded formula per order.
f.g[0] (= ∂F/∂a) must be non-zero — guaranteed by the production
solve’s strict monotonicity guard.
The expansion point a0 must be a genuine root F(a0, θ0) = 0: the
substitution recursion below cancels orders 1..=4 of G = F∘a but never
touches order 0, so a non-root a0 would yield the Taylor expansion of
the LEVEL SET F = F(a0) through a0, not the root curve F = 0. This
is guarded explicitly and re-verified by a composed-residual self-check.