CubicSpline

[kylecarow]

Supersedes #6

Boundary conditions:

• Natural (default)
• Not a knot
• Periodic
• Clamped
• CubicExtrapolate::Spline? What does this do?

[kylecarow]

CubicSpline has two different cost models depending on dimensionality right now.

1D caches: Strategy1D::init computes m once via the tridiagonal solve, and
interpolate() is just an index lookup + O(1) blend, same shape as Linear.
Strategy2D/Strategy3D/StrategyND have no init override though, so
interpolate() calls spline_eval_nd_recursive cold every time.

For an N-D grid with n points per axis, that means the innermost dimension gets
solved from scratch once per combination of the other N-1 indices (n^(N-1)
separate O(n) solves), plus each outer recursion level does its own O(n) solve on top.
A single interpolate() call is roughly O(n^N), and it's the same O(n^N) on the next
call too, since nothing carries over. Most real use cases are build-once-query-many,
so this seems like the wrong default.

Two ways to fix it, in increasing order of precompute cost:

• Cache the innermost dimension's per-line coefficients in init: one compute_m
  per 1D pencil along the last axis, done once instead of per call. Turns the
  dominant n^(N-1) solves into n^(N-1) lookups. Outer levels still solve fresh
  since their inputs depend on the query point. Roughly what neopdf's
  LogBicubicInterpolation::init does for its x-coefficients.
• Precompute the full tensor of directional derivatives at every grid point (2^N
  arrays sized like the grid), the bicubic/tricubic-with-cross-derivatives approach.
  interpolate() drops to O(N log n) regardless of grid size, at 2^N the memory and
  a heavier init.

The uncached path shouldn't stick around as an option. Nothing beats caching it once,
it's strictly worse in every workload. CubicSpline should just cache the
innermost-dimension coefficients in init for N≥2, same as 1D already does — one ND
behavior, not a choice.

Full tensor precompute is a real tradeoff though, not a strict improvement, so it
probably shouldn't fold into CubicSpline itself. 2^N memory is fine at N=2 or 3
but blows up fast past that (128x at N=7), and the heavier init cost multiplies
across however many interpolator instances get built. A consumer building one
interpolator per flavor per subgrid pays that 2^N cost dozens or hundreds of times
over. It's the same tradeoff the crate already draws a line on with hardcoded
Interp1D/2D/3D vs general InterpND, so a tensor-precomputed variant fits that
split too: a separate strategy scoped to Interp2D/Interp3D down the line, not
something InterpND needs to carry, and not blocking this PR.

Side note on the same init gap: Periodic's y[0] != y[n] check in compute_m
only sees raw grid data for the 1D case. In the recursive ND path it's checking
derived g values at one query's coordinates, not the raw boundary hyperplane — so a
malformed dataset could pass for some queries and fail for others, and either way it
won't surface as a ValidateError until the first interpolate() call, since there's
no init to catch it earlier. Worth checking Periodic/NotAKnot against the raw
grid once wherever the caching fix above lands.