Cubic interpolation family: `CubicC1` and `CubicC2`

[kylecarow]

Tracks CubicC1, the remaining piece of the cubic interpolation family originally
scoped alongside CubicC2 in #24. CubicC2 has since shipped (merged to main via
#16 plus a follow-up commit, not yet released; see CHANGELOG.md's [Unreleased]);
this issue is now just CubicC1, sized against CubicC2's finished implementation
rather than a still-hypothetical one.

Status

• CubicC2: done. Global tridiagonal solve, C², dimension-generic (1D-ND),
  configurable boundary conditions (NotAKnot/Natural/Clamped/Periodic, per-axis
  or broadcast). Fully cached at every rank: Strategy1D caches its solved
  second-derivative vector directly; Strategy2D, Strategy3D, and StrategyND all
  call the same compute_corner_cache/spline_eval_corner_cached pair
  (src/strategy/cubic.rs) to precompute and evaluate the full 2^N-wide
  corner-derivative tensor: there's no dimensionality-specific implementation left,
  StrategyND isn't a fallback path, it's the identical function Strategy2D/3D
  call, generalized to runtime N from the start. Every query after init is O(1) in
  grid size (a flat 2^N-cost Hermite blend), regardless of dimensionality or grid
  resolution.
• CubicC1: not started. This issue's scope.

Motivation

CubicC1, originally from #24: a local
cubic Hermite spline (C¹, finite-difference derivative estimate). No global solve,
cheaper to build than CubicC2, and in the same family as the recipe LHAPDF uses, so
it's not interchangeable with a natural spline for consumers wanting that kind of
local, no-global-solve behavior. CubicC1 doesn't aim for bit-for-bit LHAPDF
compatibility, though: LHAPDF's own bicubic scheme is an order-dependent,
successive-1D construction (interpolate x first, derive the y-derivative from
finite differences of the already-x-interpolated cross-sections), not a symmetric
per-knot corner-derivative spline. CubicC1 stays symmetric and order-independent,
matching CubicC2's existing architecture, so consumers wanting exact LHAPDF
parity still need their own implementation (as neopdf already has).

neopdf reimplements the LHAPDF recipe five times: LogBicubicInterpolation (2D),
LogTricubicInterpolation (3D), LogFourCubicInterpolation (4D),
LogFiveCubicInterpolation (5D), plus a standalone consolidated InterleavedHermite
whose own doc comment notes it uses "the same algorithm" as the others. Notably, only
the 2D case and InterleavedHermite cache anything (x-direction coefficients only,
src/strategy.rs/src/interleaved.rs); LogTricubicInterpolation and up cache
nothing at all, init() only validates, and every query recomputes finite differences
straight from the raw grid (src/strategy.rs's own comment: "avoids the complex 64x64
matrix"). So an uncached local-Hermite mode isn't hypothetical, it's what the actual
reference implementation ships as its default at 3D and up.

Scope

One dimension-generic strategy (impl Strategy1D/2D/3D/ND in
interpolator/{one,two,three,n}/strategies.rs), not four dimensionality-specific
types, matching CubicC2's existing shape. Caching regime is fixed at compile time
by which Interp1D/2D/3D/ND wrapper is used, never a runtime branch on grid
rank:

	Strategy1D	Strategy2D/Strategy3D	StrategyND
CubicC2 (shipped)	full cache	full corner-derivative tensor	full corner-derivative tensor
CubicC1 (this issue)	full cache	CubicC1CacheMode-dependent	CubicC1CacheMode-dependent

CubicC1 defaults to a full cache at every rank, same as CubicC2, but at
Strategy2D/3D/ND it also exposes a CacheMode choice (below) to opt out of it;
CubicC2 doesn't get an equivalent option, see the cache-mode section for why the
two aren't symmetric here.

• CubicC1 derivatives: local finite differences (central interior, one-sided
  boundary), no boundary_conditions parameter, which is what keeps it cheap to
  build (init is O(grid size), no solve), independent of the caching question
  above.
• CubicC1/CubicC2 stay separate types, not one type with a mode enum: they differ
  in init cost class (O(1) vs. a real linear solve) and continuity guarantee (C¹ vs.
  C²), unlike the boundary-condition choice within CubicC2, which is a
  same-algorithm endpoint tweak and stays a CubicC2BoundaryConditions field. Their
  respective cache-mode enums (below) stay separate types for the same reason: what
  each mode actually costs and buys differs enough between the two that a shared enum
  would either lie about one of them or force a variant neither wants.

CubicC1 derivative mode

CubicC1 carries its derivative-estimate method as a #[non_exhaustive] enum field
from v1, not bolted on later:

pub struct CubicC1<T> {
    pub derivative_mode: CubicC1DerivativeMode,
    pub cache_mode: CubicC1CacheMode,
    pub(crate) cache: ArrayD<T>,
}

#[non_exhaustive]
pub enum CubicC1DerivativeMode {
    FiniteDifference,
}

(cache_mode's type, CubicC1CacheMode, is introduced in the next section.)

Only FiniteDifference (unclipped finite differences) ships now; the field exists from the
start so Monotonic (Fritsch-Carlson clipping) can be added later as an additive enum
variant with no wire-format migration: the field, and its serialization, already
exist, so no #[serde(default)] is needed when that lands. Monotonic itself stays
out of scope until 2D/3D mixed-partial monotonicity is verified (1D Fritsch-Carlson is
well-established literature, the mixed-partial case isn't): gate it out of
Strategy2D/Strategy3D::validate until then; loosening a validation error later is
non-breaking, shipping an unverified guarantee isn't as easily undone.

CubicC1 cache mode

pub enum CubicC1CacheMode {
    Full,
    None,
}

Full (default) precomputes the corner-derivative tensor at init(), same as
CubicC2. None skips that: init() only validates, and interpolate() derives
each corner's finite differences fresh on every call. No wrapper changes are needed
for this: Interp1D::new/set_strategy keep calling init() unconditionally exactly
as they do today, CacheMode only changes what CubicC1's own init/interpolate
do internally with that call, not whether it happens. interpolate() dispatches on
self.cache_mode explicitly rather than inferring the mode from an empty cache, so
None's legitimate empty cache can't be confused with the (already-documented,
separate) misuse case of calling interpolate after deserializing without running
init_strategy.

Why caching helps CubicC1 at all: building the corner-derivative tensor is N
sequential differencing passes, mask 0 (raw values) differenced along axis 0 gives
mask 1, that field differenced along axis 1 gives masks 2 and 3, and so on, so a
k-th order mixed partial is a finite difference of a finite difference, k levels
deep. Full builds that once over the whole grid at init(), reusing each pass's
full field for the next; every query after that is just cache reads. None rebuilds
the same nested-differencing structure from scratch, scoped to the local neighborhood,
on every single call. The important property either way: that cost depends only on
N (dimension count), never on grid resolution n, so None is a bounded, safe
fallback, Full's speedup over it is real but capped by N, not the unbounded
blowup CubicC2 would see from skipping its own cache.

Why CubicC2 doesn't get an equivalent CacheMode: its only possible uncached
fallback is the pre-#16 recursive-collapse solve (there's no local, no-solve
alternative for a global spline), and that construction has a natural innermost axis
that gets revisited many times within a single query. Caching just that axis is
strictly better than not, for a ~2x memory cost, so there's no honest "less than
that" option to expose: a bare None would just be a strictly worse version of what
caching the innermost axis already gives for free. That innermost-axis-cached tier
(Partial, in the shape neopdf's LogBicubicInterpolation and the pre-#16 CubicC2
implementation both already take) is exactly what CubicC2's StrategyND measured
as, before the corner-derivative tensor was generalized to it:

D	n (per axis)	grid points	time/query
2	16	256	1.3 µs
3	16	4,096	22.8 µs
4	16	65,536	345 µs
5	16	1,048,576	5,616 µs

That's ~O(n^(D-1)) per query: every axis past the cached innermost one still has to
be re-solved in full, and correctness requires visiting every point along it, not just
local neighbors, at each recursion level. Compare a flat ~1.4 µs/query (any grid size,
n = 10 to 200 tested) once every axis is cached via the same
compute_corner_cache/spline_eval_corner_cached machinery Interp2D/Interp3D
already used (it already took grids: &[ArrayView1<T>] and derived n_axes = grids.len() at runtime, it just wasn't wired up for StrategyND yet at the time).
2^N is a dimension-count multiplier, not a grid-resolution one, it doesn't grow if
you add more grid points along an axis, only if you add axes, so full caching's
memory cost was bounded in a way Partial's per-query time cost wasn't. That's the
concrete evidence behind replacing Partial with Full for CubicC2 in the first
place; reviving Partial as an opt-in now is a separate call from having shipped it
by default, worth reconsidering only against a narrow
high-N-and-memory-constrained-and-latency-tolerant use case with no concrete request
behind it yet. Left out of this issue; CubicC2 stays Full-only for now.

Strategy1D ignores cache_mode for both types: at N=1 the cache is already just
the solved derivative vector, the same order as the values array itself, so there's
nothing meaningful to opt out of.

Cache internals: corner-derivative tensor

CubicC2 already built the shape CubicC1 reuses: ArrayD<T> sized like the value
grid, plus a trailing axis of length 2^N. Index k in that axis is a bitmask over
the N grid axes: bit i set means "∂/∂x_i of the value", so k=0 is the raw value
and k=2^N-1 is the full mixed partial. CubicC1 reuses CubicC2's evaluator
(spline_eval_corner_cached/spline_eval_corner_cached_local) unchanged; only the
population step differs:

• CubicC2 (existing, compute_corner_cache): first derivatives at each knot from a
  closed-form M → S'(x_i) formula off an already-solved moment vector (no extra
  solve); mixed partials from splining that derivative field along the other axis
  (compute_m/thomas/eval_spline_from_m, applied to a derived array instead of
  the original data). Boundary condition for that second pass: reuse bc_for_dim(j)
  for axis j, except Clamped, which falls back to Natural (the supplied value is
  a first derivative, not what's needed to spline the derivative field; Natural has
  no minimum-point requirement, unlike NotAKnot, so it can't newly reject a small
  Clamped axis that validated fine before this cache existed).
• CubicC1 (new): local finite differences at each knot and, for mixed partials, of
  the derivative field along the other axis, with no boundary condition to choose
  since there's no solve to seed.

Because both types populate the same tensor shape and share one evaluator, CubicC1
needs no new query-side code, only a new population function.

Batch interpolation caching (future, out of scope here)

batch_interpolate_into (strategy/traits.rs) defaults to a naive per-point loop.
For CubicC1 under CacheMode::None, that means every point in a batch redoes the
local finite-difference derivation from scratch, even if two points in the same batch
land in the same grid cell. A future override could memoize corner derivatives per
grid cell within one batch call (HashMap<CellIndex, [T; 2^N]> scoped to that call),
bounded by distinct cells touched in that batch rather than the whole grid, so it
never costs more than None already does and is strictly better whenever a batch
revisits cells. This is a refinement on top of CacheMode::None specifically, not a
substitute for it or for Full, and not needed for CubicC1's v1.

Naming

Decided and shipped for CubicC2 (renamed from CubicSpline on main, no
compatibility cost since it hadn't released); CubicC1 follows the same rationale.

Candidate	Problem
CubicHermite / CubicSpline	"Hermite" isn't the differentiator: a natural/not-a-knot spline is also Hermite-representable per-interval; only how per-knot derivatives are sourced differs.
Bare Cubic / CubicSpline	Reads as "normal" vs "different" cubic, not two co-equal methods.
CubicC1 / CubicC2 (chosen)	Names the actual differentiator (C¹ vs C²); doesn't collide with terms other libraries overload differently (MATLAB interp2(...,'cubic') = pchip/local, scipy interp1d(kind='cubic') = spline-based).

Rustdoc for both should still lead with "Hermite"/"spline" and literature terms
(Numerical Recipes bicubic, Lekien-Marsden tricubic, LHAPDF) so the crate stays
findable by search.

Strategy enum implications

Done, for CubicC2, in #16. strategy_enum_impl! (src/strategy/enums/mod.rs)
generates generic enums (Strategy1DEnum<T> … StrategyNDEnum<T>), threading <T>
through the enum, From impls, and trait impl block (D::Elem = T bound).
#[serde(untagged)] round-tripping is confirmed generic (tests/serde_strategies.rs).
CubicC2 is wired into all four enum files.

Remaining: wire CubicC1 into all four enum files once it exists; no further enum
rework needed, the genericity is already in place.

Explicitly out of scope

Item	Reason
Tg/Tv grid/value type split	separate issue (#57)
Per-axis heterogeneous grid types	conflicts with StrategyND's runtime-N design
Non-float/integer grid coordinates	separate issue
Global Chebyshev polynomial strategy	different interpolation family, separate issue
Generic LogSpace/GridTransformed coordinate wrapper	separate issue (#56)
Domain-specific extrapolation (LHAPDF-style)	left to downstream crates; neopdf already does this
Monotonicity-preserving (Pchip) clipping for CubicC1	plain (non-monotonic) is the minimum bar, matching neopdf's own unclipped version; the #[non_exhaustive] CubicC1DerivativeMode enum ships in v1 so this can land non-breaking later
Interpolating B-spline basis	same interpolant CubicC2 already produces, different internal algorithm (De Boor vs moment/Thomas); not a new capability
Approximating/smoothing B-spline	curve fitting, not interpolation; same reason RBF/scattered-data is already out of scope (#24)
CubicC2CacheMode / reviving the innermost-axis-cached (Partial) path for CubicC2	resurrects the pre-#16 algorithm this codebase already replaced for being too slow; no concrete request behind it yet, see the CubicC1 cache mode section above

Testing

Match CubicC2's existing coverage (src/strategy/cubic.rs unit tests,
tests/serde_strategies.rs, plus its Strategy*D integration tests): grid-point
exactness, interior accuracy against known polynomials, and serde round-trip for
CubicC1. Add cross-rank consistency coverage: CubicC1 queries via InterpND
should agree with Interp2D/Interp3D for the same rank-2/3 grid, since all three
call the same population/evaluator pair; this is as much a regression guard on that
sharing as it is new-strategy correctness testing, and mirrors what CubicC2 already
gets implicitly by using identical code across ranks. Additionally, for
CacheMode: None and Full must agree exactly for the same query (same
grid-point-exactness and known-polynomial tests, run under both modes), and serde
round-trip covers CacheMode alongside DerivativeMode.