Multi-channel ("shared-grid") interpolation support

Motivation

Interpolators (Interp1D etc) currently map grid coordinates to a single function
output values. Some real cases share one grid across several values instead:
neopdf (quark flavors), RGB images (3 color channels). Sharing the grid skips repeated
binary searches (or LinearUniform-style location evaluation) for the same point across
channels, and stores one copy of the grid axes instead of one per channel.

Alongside the existing 1-to-1 interpolators, add 1-to-many versions for Interp1D,
Interp2D, Interp3D, and InterpND. The pattern is identical across
Interp1D/2D/3D; 2D is shown below as the representative case. InterpND differs
only where runtime rank forces it to, spelled out under InterpND below.

API changes

InterpData2DMulti (+ InterpData1DMulti/InterpData3DMulti): hand-rolled, not a
generic extension of InterpData<D, N> (see Design notes for why).

#[derive(Debug, Clone)]
pub struct InterpData2DMulti<D>
where
    D: Data + RawDataClone + Clone,
    D::Elem: PartialEq + Debug,
{
    pub grid: [ArrayBase<D, Ix1>; 2],
    /// Shape `[nx, ny, n_channels]`.
    pub values: ArrayBase<D, Ix3>,
}
/// [`InterpData2DMulti`] that views data.
pub type InterpData2DMultiViewed<T> = InterpData2DMulti<ViewRepr<T>>;
/// [`InterpData2DMulti`] that owns data.
pub type InterpData2DMultiOwned<T> = InterpData2DMulti<OwnedRepr<T>>;

/// Hand-written, not derived, mirroring `InterpData<D, N>`'s own `PartialEq` impl
/// exactly (`src/interpolator/data.rs`): required by `partialeq_impl!` below, which
/// bounds on `InterpData2DMulti<D>: PartialEq`. Without this, that bound is never
/// satisfied and the `PartialEq for Interp2DMulti<D, S>` impl it generates would be
/// present but permanently unusable.
impl<D> PartialEq for InterpData2DMulti<D>
where
    D: Data + RawDataClone + Clone,
    D::Elem: PartialEq + Debug,
    ArrayBase<D, Ix1>: PartialEq,
{
    fn eq(&self, other: &Self) -> bool {
        self.grid == other.grid && self.values == other.values
    }
}

impl<D> InterpData2DMulti<D>
where
    D: Data + RawDataClone + Clone,
    D::Elem: PartialOrd + Debug,
{
    /// Same checks as `InterpData<D, N>::validate` (grid length, monotonicity,
    /// grid/values shape agreement), plus a non-empty channel axis: a grid with no
    /// channels has nothing for `Strategy2DMulti::interpolate_into` to write.
    pub fn validate(&self) -> Result<(), ValidateError> {
        if self.n_channels() == 0 {
            return Err(ValidateError::Other(
                "InterpData2DMulti requires at least 1 channel".to_string(),
            ));
        }
        for i in 0..2 {
            let i_grid_len = self.grid[i].len();
            if i_grid_len < 2 {
                return Err(ValidateError::InsufficientGridPoints(i));
            }
            if !self.grid[i].windows(2).into_iter().all(|w| w[0] <= w[1]) {
                return Err(ValidateError::NonMonotonic(i));
            }
            if i_grid_len != self.values.shape()[i] {
                return Err(ValidateError::IncompatibleShapes(i));
            }
        }
        Ok(())
    }

    pub fn n_channels(&self) -> usize {
        self.values.shape()[2]
    }

    /// Borrow channel `k` as a standalone, viewed `InterpData2D`.
    ///
    /// # Note
    /// The returned `values` view is **strided, not contiguous**: with channels last,
    /// `index_axis(Axis(2), k)` steps by `n_channels`. A strategy that reaches for
    /// `data.values.as_slice().unwrap()` will panic on it. Index via `ArrayView`
    /// instead, the same guidance `Strategy2D::interpolate` already carries for
    /// `Interp*Viewed`.
    pub fn channel_view(&self, k: usize) -> InterpData2DViewed<&D::Elem> {
        InterpData2D {
            grid: std::array::from_fn(|i| self.grid[i].view()),
            values: self.values.index_axis(Axis(2), k),
        }
    }
}

InterpData1DMulti/InterpData3DMulti: grid: [ArrayBase<D, Ix1>; 1 or 3], values: ArrayBase<D, Ix2 or Ix4>, same three methods, Owned/Viewed aliases, and hand-written
PartialEq impl.

Channel axis is last ([nx, ny, n_channels]) so that evaluating every channel at one
point walks contiguous memory in C order. This is the same insight behind neopdf's
hand-rolled [cell][flavor][4] interleaved layout. The cost is that the single-channel
view is the strided one, which is the right trade given the type exists for the
all-channel case.

Strategy2DMulti: independent trait, no supertrait relationship to Strategy2D (see
Design notes for why not).

pub trait Strategy2DMulti<D>: Debug + DynClone
where
    D: Data + RawDataClone + Clone,
    D::Elem: PartialEq + Debug,
{
    fn validate(&self, _data: &InterpData2DMulti<D>) -> Result<(), ValidateError> {
        Ok(())
    }

    fn init(&mut self, _data: &InterpData2DMulti<D>) -> Result<(), ValidateError> {
        Ok(())
    }

    /// Interpolate every channel at `point`, writing channel `k` into `out[k]`.
    /// `out.len()` must equal `data.n_channels()`.
    ///
    /// The only required method. Per #45, the out-slice form is what strategies
    /// implement and the allocating forms below are defaulted wrappers over it: a
    /// multi-channel result is a `Vec` per *point*, so a `Vec`-returning required
    /// method would put an allocation in every single-point call, not just per batch.
    fn interpolate_into(
        &self,
        data: &InterpData2DMulti<D>,
        point: &[D::Elem; 2],
        out: &mut [D::Elem],
    ) -> Result<(), InterpolateError>;

    /// Unchecked [`Strategy2DMulti::interpolate_into`]. Default just unwraps it.
    fn interpolate_fast_into(
        &self,
        data: &InterpData2DMulti<D>,
        point: &[D::Elem; 2],
        out: &mut [D::Elem],
    ) { /* default: unwrap interpolate_into */ }

    /// Interpolate every channel at each of several points, sharing one grid across
    /// all of them. `out.len()` must equal `points.len() * data.n_channels()`; point
    /// `p`'s channels occupy `out[p * n_channels ..][.. n_channels]`.
    ///
    /// Flat, not `Vec<Vec<_>>`: nesting costs one allocation per point and forces a
    /// transpose on any caller that wants channel-major output. Default loops
    /// [`Strategy2DMulti::interpolate_into`] over the chunks. Override only if
    /// locating a point can be amortized across both channels *and* the batch at
    /// once, e.g. a Chebyshev-style strategy computing one coefficient matrix
    /// covering every point and every channel in a single pass.
    fn batch_interpolate_into(
        &self,
        data: &InterpData2DMulti<D>,
        points: &[[D::Elem; 2]],
        out: &mut [D::Elem],
    ) -> Result<(), InterpolateError> { /* default: chunk out, loop interpolate_into */ }

    fn batch_interpolate_fast_into(
        &self,
        data: &InterpData2DMulti<D>,
        points: &[[D::Elem; 2]],
        out: &mut [D::Elem],
    ) { /* default: chunk out, loop interpolate_fast_into */ }

    /// Allocating [`Strategy2DMulti::interpolate_into`]. Defaulted, do not override.
    fn interpolate(
        &self,
        data: &InterpData2DMulti<D>,
        point: &[D::Elem; 2],
    ) -> Result<Vec<D::Elem>, InterpolateError> { /* default: alloc, call _into */ }

    // ... `interpolate_fast`, `batch_interpolate` (-> Vec<D::Elem>, flat, same
    // chunking as above), `batch_interpolate_fast`, all defaulted the same way.
}

validate/init default exactly like Strategy2D's own. Every built-in strategy needs
a real, hand-written interpolate_into, not a mechanical loop over its scalar
interpolate: each one has a per-axis locate step (nearest index, step-direction index,
fractional blend position) that a naive per-channel loop would re-run once per channel for
the same point, exactly the repeated-search cost this issue exists to avoid.

Nearest locates each axis once (already how its scalar interpolate works, via
locate_lower_index plus a distance comparison to the nearest grid index [i, j]), then
loops channels for a direct data.values[[i, j, k]] lookup, no blending.

Linear locates each axis once via locate_axis (#29, src/strategy/utils.rs, already
on main), then reuses that result to blend all n_channels values:

impl<D> Strategy2DMulti<D> for Linear
where
    D: Data + RawDataClone + Clone,
    D::Elem: Float + Debug,
{
    fn interpolate_into(
        &self,
        data: &InterpData2DMulti<D>,
        point: &[D::Elem; 2],
        out: &mut [D::Elem],
    ) -> Result<(), InterpolateError> {
        let locations = std::array::from_fn(|i| locate_axis(data.grid[i].view(), &point[i]));
        for (k, o) in out.iter_mut().enumerate() {
            *o = blend_from_locations(&locations, &data.values.index_axis(Axis(2), k));
        }
        Ok(())
    }
}

blend_from_locations is the same bilinear blend Strategy2D for Linear (already on
main, src/interpolator/two/strategies.rs) does today, factored out of that impl into
a pub(crate) free function in src/strategy/utils.rs so it can run once per channel
against a pre-computed locations instead of each channel re-deriving its own via
locate_axis. Signature:

pub(crate) fn blend_from_locations<T: Float>(
    locations: &[AxisLocation<T>; 2],
    values: &ArrayView2<T>,
) -> T

Body is the existing 4-way Exact/Interp match, verbatim, taking locations as an
argument rather than deriving them from point.

Strategy2D for Linear's own interpolate is left as-is, not required to route through
this helper too: nothing needs it to for correctness, since the two impls share only the
shape of the math, not state. Worth doing later as a dedupe pass, not part of this issue.

Step/StepLower/StepUpper/LinearUniform follow the same overall shape: locate once
via whichever helper they already use (locate_step_index/locate_lower_index_uniform),
then loop channels for the lookup or blend. None of them need a shared helper the way
Linear does, since their existing single-point logic is already a plain index lookup or
a single two-point blend, not a 4-way match worth factoring out.

A stateful strategy (illustrative only, not proposed here) would cache per-channel state
in init by looping data.channel_view(k) for k in 0..data.n_channels(), then still
locate once per point in interpolate_into and blend against each channel's cached
state. The load-bearing detail: init sees all channels at once (&InterpData2DMulti<D>),
which is what makes a joint layout across channels expressible, not just a Vec of
independent per-channel states. See Design notes.

Interp2DMulti wrapper (mirrors Interp2D, reuses its macros):

#[derive(Debug, Clone)]
pub struct Interp2DMulti<D, S>
where
    D: Data + RawDataClone + Clone,
    D::Elem: PartialEq + Debug,
    S: Strategy2DMulti<D> + Clone,
{
    pub data: InterpData2DMulti<D>,
    pub strategy: S,
    pub extrapolate: Extrapolate<D::Elem>,
}
pub type Interp2DMultiViewed<T, S> = Interp2DMulti<ViewRepr<T>, S>;
pub type Interp2DMultiOwned<T, S> = Interp2DMulti<OwnedRepr<T>, S>;

extrapolate_impl!(Interp2DMulti, Strategy2DMulti);
partialeq_impl!(Interp2DMulti, InterpData2DMulti, Strategy2DMulti);

new mirrors Interp2D::new exactly: data.validate()?, then check_extrapolate, then
strategy.validate, then strategy.init. (D::Elem: Float, this impl block's own bound,
already implies the PartialOrd that validate needs.)

Inherent methods, matching the trait's method set: interpolate_into,
interpolate_fast_into, batch_interpolate_into, batch_interpolate_fast_into, plus the
four allocating wrappers. No _multi suffix on any of them, the struct name already
carries it.

The extrapolate handling is the same per-axis logic Interp2D::interpolate and
Interp2D::batch_interpolate (#21) already have, adapted in exactly one way:
Extrapolate::Fill(value) writes value into all n_channels slots for that point
rather than a single slot. Enable/Clamp (unconditional point transform), Wrap
(conditional on out_of_bounds, since it is not identity at the boundary), and Error
(aggregating every offending point and dimension into one ExtrapolateError) are
unchanged. The batch version's per-mode partitioning is batch_interpolate_impl! in
src/interpolator/mod.rs; parameterizing that macro over channel count is preferable to
duplicating it per Interp*Multi type.

InterpND. Same design, with three differences forced by runtime rank:

pub struct InterpDataNDMulti<D> {
    pub grid: Vec<ArrayBase<D, Ix1>>,
    /// Rank `grid.len() + 1`; the trailing axis is channels.
    pub values: ArrayBase<D, IxDyn>,
}

1. Rank check replaces the type-level guarantee. InterpDataND::validate can assume
   grid.len() == values.ndim(). Here validate must first check values.ndim() == grid.len() + 1 and return ValidateError::IncompatibleShapes otherwise, since nothing
   in the type distinguishes a rank-n multi grid from a rank-n single grid.
   n_channels() is values.shape()[self.grid.len()], and channel_view(k) is
   values.index_axis(Axis(self.grid.len()), k).
2. Slice-typed points. StrategyNDMulti::interpolate_into(&self, data, point: &[D::Elem], out: &mut [D::Elem]), matching StrategyND. InterpNDMulti's batch
   methods take &[&[D::Elem]].
3. Hand-written batch methods. There is no fixed N to route an inherent
   array-typed method through, so InterpNDMulti::batch_interpolate_into gets the same
   hand-written treatment InterpND::batch_interpolate already got in #21, rather than
   going through the shared macro.

InterpNDMulti is not optional scope. A downstream consumer with more than three grid
axes has no other route, and neopdf specifically routes 7 of its 13 grid configurations
through InterpND (see Downstream motivation).

InterpolatorMulti<T> (src/interpolator/mod.rs), the multi-channel analog of
Interpolator<T>, filling out the same structural parallel Interp*Multi/Strategy*Multi
already commit to:

pub trait InterpolatorMulti<T>: DynClone {
    fn ndim(&self) -> usize;
    fn n_channels(&self) -> usize;
    fn validate(&self) -> Result<(), ValidateError>;
    fn set_extrapolate(&mut self, extrapolate: Extrapolate<T>) -> Result<(), ValidateError>;

    /// Interpolate every channel at `point`, writing into `out`. `out.len()` must equal
    /// `self.n_channels()`. The required method, per #45.
    fn interpolate_into(&self, point: &[T], out: &mut [T]) -> Result<(), InterpolateError>;

    fn interpolate_fast_into(&self, point: &[T], out: &mut [T]) {
        self.interpolate_into(point, out)
            .expect("interpolate_fast_into: invalid point or data")
    }

    /// Default chunks `out` by `n_channels` and loops [`InterpolatorMulti::interpolate_into`].
    fn batch_interpolate_into(
        &self,
        points: &[&[T]],
        out: &mut [T],
    ) -> Result<(), InterpolateError> { /* default: chunk out, loop interpolate_into */ }

    fn batch_interpolate_fast_into(&self, points: &[&[T]], out: &mut [T]) { /* same, _fast_into */ }

    /// Allocating [`InterpolatorMulti::interpolate_into`]. Defaulted, do not override.
    fn interpolate(&self, point: &[T]) -> Result<Vec<T>, InterpolateError> { /* default: alloc, call interpolate_into */ }

    // ... `interpolate_fast`, `batch_interpolate`, `batch_interpolate_fast`, all
    // defaulted the same way: allocate, call the `_into` form.
}

clone_trait_object!(<T> InterpolatorMulti<T>);

impl<T> InterpolatorMulti<T> for Box<dyn InterpolatorMulti<T>> {
    // forwards every method to `(**self)`, same shape as `Interpolator<T>`'s own
    // `Box<dyn Interpolator<T>>` impl
}

Implemented for Interp1DMulti/2DMulti/3DMulti/NDMulti<D, S>, both Owned and
Viewed (unlike DynInterpolatorMulti below, nothing here needs 'static): each
interpolate_into/batch_interpolate_into converts the incoming slice(s) to the fixed
&[T; N]/&[[T; N]] the inherent methods take first, the same shadowing #39/#21 force on
the scalar Interpolator<T> impls, then forwards. InterpNDMulti needs no conversion,
same as InterpND. Bounds match each Interp*Multi's own Interpolator-equivalent impl
(e.g. D::Elem: Num + PartialOrd + Euclid + Copy + Debug for 2D, mirroring Interp2D's).

Every method needs an explicit override in each impl, not just interpolate_into,
for the reason #21/#46 give for Interpolator<T>: a defaulted body called through Box<dyn InterpolatorMulti<T>> dispatches through this trait's vtable once per point, instead of
reaching the concrete type's real implementation in one.

DynInterpolatorMulti, extending InterpolatorMulti<T> the same way #46's
DynInterpolator extends Interpolator<T>:

pub trait DynInterpolatorMulti<T>: InterpolatorMulti<T> + Send + Sync {
    fn as_any(&self) -> &dyn Any;
}

Blanket impls for Interp1DMultiOwned/Interp2DMultiOwned/Interp3DMultiOwned/
InterpNDMultiOwned only: as_any requires Self: 'static, which Interp*MultiViewed
can't satisfy. Each impl is just fn as_any(&self) -> &dyn Any { self }; n_channels/
interpolate_into/etc. are all inherited from InterpolatorMulti<T> for free.

Where this lives

Mirrors the existing per-dimensionality split (mod.rs/strategies.rs/tests.rs under
one/two/three/n); 2D shown, others identical.

File	Contents
src/interpolator/two/multi.rs (new)	InterpData2DMulti/Owned/Viewed, Interp2DMulti/Owned/Viewed, extrapolate_impl!/partialeq_impl! invocations, impl InterpolatorMulti<T> for Interp2DMulti<D, S>, impl DynInterpolatorMulti<T> for Interp2DMultiOwned<T, S>
src/interpolator/two/mod.rs	+ mod multi; and re-exports of its public types
src/interpolator/two/strategies.rs	+ impl Strategy2DMulti<D> for Nearest/Linear/Step/StepLower/StepUpper/LinearUniform
src/strategy/traits.rs	+ Strategy1DMulti/2DMulti/3DMulti/NDMulti definitions
src/strategy/utils.rs	+ pub(crate) blend_from_locations
src/interpolator/data.rs	+ pub use two::{InterpData2DMulti, InterpData2DMultiOwned, InterpData2DMultiViewed}; and equivalents, mirroring the existing non-multi line
src/interpolator/mod.rs	+ InterpolatorMulti<T> definition, Box<dyn InterpolatorMulti<T>> forwarding impl, and clone_trait_object! call, next to Interpolator<T>; + DynInterpolatorMulti<T> definition, next to DynInterpolator (#46); + re-exports of Interp{1,2,3}DMulti{,Owned,Viewed}, InterpNDMulti{,Owned,Viewed}, InterpolatorMulti, DynInterpolatorMulti; parameterize batch_interpolate_impl! over channel count
src/lib.rs prelude	+ the Interp*Multi types and InterpolatorMulti, joining Interp*D/Interpolator today. InterpData*Multi stays out of prelude, matching InterpData2D today: fully public via interpolator::data, just not in the curated re-export. DynInterpolatorMulti stays out too, same reasoning as #46's DynInterpolator (see Design notes). Strategy*Multi needs no line either, pub use crate::strategy; already covers it.

Strategy2DEnum (src/strategy/enums/two.rs) is untouched; no Strategy2DMultiEnum or
InterpolatorMultiEnum proposed here, see Non-goals.

Design notes

• Hand-rolled, not InterpData<D, N> generalized: computing Dim<[Ix; N + 1]> from
  const N: usize needs unstable generic_const_exprs. Not available on stable. Does
  not apply to InterpDataNDMulti, where IxDyn already carries rank at runtime, which
  is why that one is a validation check rather than a type-level guarantee.
• interpolate is a separate method, not a unified signature: forcing the scalar path
  to return Vec<D::Elem> would tax every single-channel call with a heap allocation.
• interpolate_into is the required method, not interpolate: see #45. The
  argument is sharper here than for the scalar batch case, because a multi-channel result
  is per-point rather than per-batch.
• Strategy2DMulti does not extend Strategy2D: an earlier draft did
  (Strategy2DMulti<D>: Strategy2D<D>), with interpolate's default delegating to the
  inherited scalar interpolate per channel. That breaks two ways. The scalar method has
  no channel index, so a strategy caching per-channel or joint coefficients cannot tell
  which channel it is being asked about. And the inherited validate/init take
  &InterpData2D<D>, single-channel-shaped, so there is no correct way to call them
  against InterpData2DMulti<D>: one arbitrary channel ignores the rest, and looping the
  scalar method against the same &mut self makes the last channel silently win.
  Dropping the bound fixes both, at the cost of no blanket/macro opt-in for stateless
  strategies. That cost is zero in practice: no strategy in this crate is stateless
  enough to want the mechanical per-channel loop, every one has a locate step worth
  sharing. This is not a theoretical concern either, see Downstream motivation.
• InterpolatorMulti<T> mirrors Interpolator<T>, not Strategy*Multi: its required
  method is the allocation-free interpolate_into, not an allocating interpolate, same
  reasoning as Strategy*Multi's own required method above: a multi-channel result is
  per-point, so an allocating required method taxes every call, not just batches. Every
  other method (interpolate, _fast, batch_*) is a defaulted wrapper over it, mirroring
  the shape Strategy*Multi already commits to. n_channels lives on the trait itself, not
  just DynInterpolatorMulti: even a non-erased Box<dyn InterpolatorMulti<T>> caller has
  to size out before calling, without downcasting.
• DynInterpolatorMulti extends InterpolatorMulti<T>: same reasoning #46 gives for
  DynInterpolator extending Interpolator<T>: the borrowed Interp*MultiViewed types
  simply don't implement DynInterpolatorMulti (still blocked by as_any needing Self: 'static), InterpolatorMulti<T> itself is untouched. Collapses DynInterpolatorMulti to
  one method, as_any; Send/Sync stay scoped per-impl for the same reason a custom
  Strategy2D (examples/custom_strategy.rs) may hold non-thread-safe state.
• InterpolatorMulti<T> joins the prelude, DynInterpolatorMulti stays out: mirrors
  #46's identical split for Interpolator<T>/DynInterpolator. prelude is curated for
  the common path; heterogeneous storage + downcasting (the neopdf case) is a narrower,
  advanced use case, one explicit use away for consumers who need it.
• Flat batch output, not Vec<Vec<T>>: nesting allocates per point and fixes an
  orientation on the caller. A flat out with documented n_channels chunking lets the
  caller own the layout.
• Naming: type and trait names carry a trailing Multi throughout; none of their own
  methods repeat it, Multi in the name already says so, InterpolatorMulti/
  DynInterpolatorMulti included, matching how #46 settled the same question for
  DynInterpolator (the trait, not a method-name suffix, is what disambiguates a type-erased
  call). batch_ is a prefix, not a suffix like _fast/_into, for the reason #21 gives: it
  changes what is being operated on (many points instead of one), rather than being a
  variant of the same operation.

Downstream motivation

QCDLab/neopdf is the driving consumer, and it has moved well past hand-rolling one
trait. As of aeb45a0 it maintains two complete shared-grid multi-channel evaluators
that bypass ninterp entirely, both structured exactly as this issue proposes:

neopdf	Shape
InterleavedHermite (neopdf/src/interleaved.rs, 508 lines)	locate() once per point, then eval_allpids(&loc, pid_slots, force_positive_fn, out: &mut [f64]) over flavors. Hermite x-coefficients precomputed at build time into a [cell][flavor][4] interleaved layout.
ChebyshevAllPids (neopdf/src/strategy.rs)	Same split. locate() returns barycentric coefficients per dimension, eval_allpids contracts them per flavor.

Both are reached before the ninterp-backed path in GridPDF::xfxq2_allpids, with
Vec<Vec<Box<dyn DynInterpolator>>> as the generic fallback. The ninterp path is already
the slow path for every grid type neopdf ships.

Three things follow.

The trait shape is validated by working code. InterleavedHermite's interleaved
coefficient layout is per-channel state built jointly across all channels, which is
expressible as Strategy2DMulti::init precisely because init sees
&InterpData2DMulti<D>. Under the rejected supertrait design, with init taking
&InterpData2D<D> per channel, it would not have been. The Design note above is not
hypothetical.

Memory duplication is the unglamorous win. InterpolatorFactory::create is called
once per (subgrid, flavor) and does subgrid.grid_slice(pid_index).to_owned() plus its
own subgrid.xs.mapv(f64::ln) and q2s.mapv(f64::ln). So the log-transformed axis arrays
are recomputed and stored once per flavor. On top of that, GridPDF retains the original
knot_array and the interleaved coefficients (themselves 4 floats per cell per flavor).
One InterpData*Multi per subgrid collapses the axis duplication outright and makes the
values a single array.

Dimensional coverage matters. InterpolationConfig has 13 variants routed to
Interp2D (1), Interp3D (5), and InterpND (7). Both fast paths have build arms for
2D through 5D. Shipping this issue without InterpNDMulti would leave both hand-rolled
evaluators alive for the 4D and 5D configurations, which is why ND is in scope here rather
than deferred.

Not absorbed by this issue, correctly: neopdf's cross-subgrid point routing (grouping a
batch by which of several interpolators each point falls into before calling the batch
method on each) is PDF-domain logic, and its force_positive clipping is a post-map.

Non-goals

• No unification of scalar and multi-channel trait method signatures.
• No const-generic channel-count parameter: encoding it into the array rank on stable
  needs the same unstable generic_const_exprs as the hand-rolled-struct decision above.
  Channel count is a runtime axis size instead.
• No per-channel strategy selection: sharing one locate step per point requires one
  strategy for all channels. Separate Interp2D/Interp2DViewed instances already cover
  channels that genuinely need different strategies. Same reasoning rules out a
  Vec<S>/[S; N] of per-channel strategy instances even of the same type: it
  reintroduces the naive per-channel cost.
• No channel subsetting, tracked separately in #47. interpolate_into here
  always evaluates every channel. Subsetting is additive on top of this trait (defaultable
  in terms of interpolate_into), but it should not lag far behind: without it, a
  consumer with both single-channel and all-channel access patterns has to keep
  per-channel Interp2D instances alongside Interp2DMulti, which gives back the memory
  win this issue is partly here for.
• No Box<dyn Strategy1DMulti/2DMulti/3DMulti/NDMulti<D>> support in this pass, and
  therefore no forwarding concern for one. Interp2D-style boxed-strategy runtime
  swapping was never proposed for the *Multi wrappers. Revisit if it is added later,
  following #21's Box<dyn Strategy1D/2D/3D/ND<D>> precedent exactly.
• No Strategy*MultiEnum/InterpolatorMultiEnum, the *Multi counterpart to
  Strategy*Enum/InterpolatorEnum (src/interpolator/enums.rs, src/strategy/enums/).
  InterpolatorEnum's variants are Interp1D<D, Strategy1DEnum> etc., so this needs four
  new Strategy*MultiEnum types before an InterpolatorMultiEnum wrapping them is even
  possible, doubling the existing enum-dispatch surface. That module already hand-rolls
  what enum_dispatch would give for free if it supported a generic trait on a
  non-generic enum (see the NOTE at the top of enums.rs); doubling the current
  boilerplate before addressing that seems like the wrong order. Box<dyn InterpolatorMulti<T>>/Box<dyn DynInterpolatorMulti<T>> cover the runtime-polymorphism
  need in the meantime, same as Box<dyn Interpolator<T>> did before InterpolatorEnum
  existed.

Dependencies

• #29 (closed, merged): locate_axis/AxisLocation are on main.
• #21 (closed, merged): batch_interpolate/batch_interpolate_fast and the
  batch_interpolate_impl! macro this issue parameterizes.
• #39 (closed, merged): the array-based inherent interpolate whose name shadowing
  forces the try_into conversions in the blanket impls.
• #45: establishes the _into convention. Should land first; this issue's
  signatures assume it.
• #46: defines the scalar Interpolator<T>/DynInterpolator<T> pair in
  src/interpolator/mod.rs, the precedent InterpolatorMulti<T>/DynInterpolatorMulti<T>
  follow here (subtrait extension, not a redeclared method set). Independent in both
  directions.
• Followed by #47.

Start by reading the existing data and wrapper patterns in src/interpolator/data.rs, then inspect src/interpolator/two/strategies.rs and src/strategy/utils.rs, including locate_axis and the existing Linear blend. Implement the 1D/2D/3D/ND multi-channel data, strategies, wrappers, and built-in strategy support described in the issue; done means the APIs validate, interpolate into channel-major output, and preserve the specified extrapolation behavior.