fix(nd): value-matched grid promotion

src/akima/akima_oneshot.jl

@@ -24,9 +24,11 @@
     # Grid pre-normalized by the public `akima_interp` API via `_resolve_axis(x)`
     # before dispatching here; `_periodic_extend_1d` preserves the normalization.
     x_eff, y_ext, bc_eff, extrap_eff = _periodic_extend_1d(x, y, bc, extrap)
-    Tdy = _promote_eltype(_coeff_op, eltype(x_eff), Tv)
+    # Value-matched width: dy buffer + slope arithmetic run at `Tw` — see pchip_oneshot.jl.
+    Tw = _promote_grid_float(eltype(x_eff), Tv)
+    Tdy = _promote_eltype(_coeff_op, Tw, Tv)
     dy = acquire!(pool, Tdy, length(y_ext))
-    _akima_slopes!(dy, x_eff, y_ext; bc = bc_eff)
+    _akima_slopes!(dy, x_eff, y_ext, Tw; bc = bc_eff)
     searcher = _resolve_search(x_eff, xq, search, hint)
     return _hermite_eval_at_point(x_eff, y_ext, dy, xq, extrap_eff, deriv, searcher)
 end
@@ -47,9 +49,11 @@ end
     @boundscheck length(output) == length(x_query) || _throw_length_mismatch(length(x_query), length(output), "x_query", "output")
     x_eff, y_ext, bc_eff, extrap_eff = _periodic_extend_1d(x, y, bc, extrap)
 
-    Tdy = _promote_eltype(_coeff_op, eltype(x_eff), Tv)
+    # Value-matched width: dy buffer + slope arithmetic run at `Tw` — see pchip_oneshot.jl.
+    Tw = _promote_grid_float(eltype(x_eff), Tv)
+    Tdy = _promote_eltype(_coeff_op, Tw, Tv)
     dy = acquire!(pool, Tdy, length(y_ext))
-    _akima_slopes!(dy, x_eff, y_ext; bc = bc_eff)
+    _akima_slopes!(dy, x_eff, y_ext, Tw; bc = bc_eff)
     searcher = _resolve_search(x_eff, x_query, search, hint)
     return _hermite_vector_loop!(output, x_eff, y_ext, dy, x_query, extrap_eff, deriv, searcher)
 end
@@ -126,7 +130,8 @@ Outlier-robust, C\$^1\$ continuous.
         search::AbstractSearchPolicy = AutoSearch(),
         hint::Union{Nothing, Base.RefValue{Int}} = nothing
     ) where {Tg, Tv, Tq <: Real}
-    x = _resolve_axis(x)
+    # Value-matched Tg: Int/OneTo grid + Float32 data → Float32 axis.
+    x = _resolve_axis(x, _promote_grid_float(Tg, Tv))
     extrap_eff = _resolve_extrap(extrap, bc, x, y)
     resolved = _resolve_coeffs(coeffs, x, xq)
     if resolved isa OnTheFly
@@ -152,7 +157,7 @@ In-place Akima interpolation with outlier-robust slopes.
         search::AbstractSearchPolicy = AutoSearch(),
         hint::Union{Nothing, Base.RefValue{Int}} = nothing
     ) where {Tg, Tv, Tq <: Real}
-    x = _resolve_axis(x)
+    x = _resolve_axis(x, _promote_grid_float(Tg, Tv))
     extrap_eff = _resolve_extrap(extrap, bc, x, y)
     resolved = _resolve_coeffs(coeffs, x, x_query)
     if resolved isa OnTheFly

src/akima/akima_slopes.jl

@@ -40,12 +40,18 @@ secant sequence.
 # Complexity
 O(n), single pass, zero allocation (writes into `dy`).
 """
+# Width-less form: delegate with the axis's own eltype — bit-identical to the
+# historic raw behavior. One-shot PreCompute backends pass the value-matched `Tw`.
+_akima_slopes!(dy::AbstractVector, x::AbstractVector, y::AbstractVector; bc::AbstractBC = NoBC()) =
+    _akima_slopes!(dy, x, y, eltype(x); bc)
+
 function _akima_slopes!(
         dy::AbstractVector,
-        x::AbstractVector{Tg},
-        y::AbstractVector;
+        x::AbstractVector,
+        y::AbstractVector,
+        ::Type{Tw};
         bc::AbstractBC = NoBC()
-    ) where {Tg}
+    ) where {Tw}
     n = length(x)
     @assert n >= 2 "Akima requires at least 2 points"
     @assert length(y) == n "y length must match x"
@@ -55,12 +61,12 @@ function _akima_slopes!(
     # 4-secant helper (cycle=2 for `:exclusive` yields a 2-secant alternation).
     if n == 2
         if bc isa PeriodicBC
-            @inbounds dy[1] = _akima_local_4secant_periodic(x, y, 1, n, bc)
-            @inbounds dy[2] = _akima_local_4secant_periodic(x, y, 2, n, bc)
+            @inbounds dy[1] = _akima_local_4secant_periodic(Tw, x, y, 1, n, bc)
+            @inbounds dy[2] = _akima_local_4secant_periodic(Tw, x, y, 2, n, bc)
             return dy
         end
         @inbounds begin
-            δ = _forward_secant(x, y, 1)
+            δ = _forward_secant(Tw, x, y, 1)
             dy[1] = δ
             dy[2] = δ
         end
@@ -73,13 +79,13 @@ function _akima_slopes!(
             # Wrap-aware path: every index is within K=5 stencil reach of the
             # join, so use the closed-cycle 4-secant formula at all 3 points.
             @inbounds for i in 1:3
-                dy[i] = _akima_local_4secant_periodic(x, y, i, n, bc)
+                dy[i] = _akima_local_4secant_periodic(Tw, x, y, i, n, bc)
             end
             return dy
         end
         @inbounds begin
-            m1 = _forward_secant(x, y, 1)
-            m2 = _forward_secant(x, y, 2)
+            m1 = _forward_secant(Tw, x, y, 1)
+            m2 = _forward_secant(Tw, x, y, 2)
             dy[1] = m1
             dy[2] = (m1 + m2) / 2
             dy[3] = m2
@@ -101,9 +107,9 @@ function _akima_slopes!(
     # Simplification: extrapolate m sequence linearly.
 
     # Compute all n-1 secant slopes
-    @inbounds m1 = _forward_secant(x, y, 1)
-    @inbounds m2 = _forward_secant(x, y, 2)
-    @inbounds m3 = _forward_secant(x, y, 3)
+    @inbounds m1 = _forward_secant(Tw, x, y, 1)
+    @inbounds m2 = _forward_secant(Tw, x, y, 2)
+    @inbounds m3 = _forward_secant(Tw, x, y, 3)
 
     # Boundary virtual / wrapped secants for the LEFT side of the domain.
     # NoBC:                          linear extrapolation (Akima's original).
@@ -112,8 +118,8 @@ function _akima_slopes!(
     #                                (closed cycle on n cells, with virtual seam cell).
     # The `_periodic_secant` abstraction absorbs both PeriodicBC variants.
     if bc isa PeriodicBC
-        @inbounds m_0 = _periodic_secant(x, y, 0, n, bc)     # m[0]
-        @inbounds m_neg1 = _periodic_secant(x, y, -1, n, bc)    # m[-1]
+        @inbounds m_0 = _periodic_secant(Tw, x, y, 0, n, bc)     # m[0]
+        @inbounds m_neg1 = _periodic_secant(Tw, x, y, -1, n, bc)    # m[-1]
     else
         @inbounds m_0 = 2 * m1 - m2                          # virtual (NoBC)
         @inbounds m_neg1 = 3 * m1 - 2 * m2                      # virtual (NoBC)
@@ -133,7 +139,7 @@ function _akima_slopes!(
     m_k = m3
 
     @inbounds for k in 3:(n - 2)
-        m_kp1 = _forward_secant(x, y, k + 1)
+        m_kp1 = _forward_secant(Tw, x, y, k + 1)
         dy[k] = _akima_weighted_slope(m_km2, m_km1, m_k, m_kp1)
         m_km2 = m_km1
         m_km1 = m_k
@@ -147,8 +153,8 @@ function _akima_slopes!(
     # PeriodicBC{:inclusive}:        m[n]=m[1], m[n+1]=m[2] (closed cycle on n-1 cells).
     # PeriodicBC{:exclusive}:        m[n]=seam, m[n+1]=m[1] (closed cycle on n cells).
     if bc isa PeriodicBC
-        @inbounds m_np1 = _periodic_secant(x, y, n, n, bc)         # m[n]
-        @inbounds m_np2 = _periodic_secant(x, y, n + 1, n, bc)     # m[n+1]
+        @inbounds m_np1 = _periodic_secant(Tw, x, y, n, n, bc)         # m[n]
+        @inbounds m_np2 = _periodic_secant(Tw, x, y, n + 1, n, bc)     # m[n+1]
     else
         m_np1 = 2 * m_k - m_km1                                    # virtual (NoBC)
         m_np2 = 3 * m_k - 2 * m_km1                                # virtual (NoBC)

src/cardinal/cardinal_oneshot.jl

@@ -30,9 +30,12 @@
     # `_CachedRange`/Vector) — the public `cardinal_interp` API pre-resolved
     # via `_resolve_axis(x)` before dispatching here, so no extra prep needed.
     x_eff, y_ext, bc_eff, extrap_eff = _periodic_extend_1d(x, y, bc, extrap)
-    Tdy = _promote_eltype(_coeff_op, eltype(x_eff), Tv)
+    # Value-matched width: dy buffer + slope arithmetic (incl. the `1 - tension`
+    # scale) run at `Tw` — see pchip_oneshot.jl.
+    Tw = _promote_grid_float(eltype(x_eff), Tv)
+    Tdy = _promote_eltype(_coeff_op, Tw, Tv)
     dy = acquire!(pool, Tdy, length(y_ext))
-    _cardinal_slopes!(dy, x_eff, y_ext, tension; bc = bc_eff)
+    _cardinal_slopes!(dy, x_eff, y_ext, tension, Tw; bc = bc_eff)
     searcher = _resolve_search(x_eff, xq, search, hint)
     return _hermite_eval_at_point(x_eff, y_ext, dy, xq, extrap_eff, deriv, searcher)
 end
@@ -54,9 +57,12 @@ end
     @boundscheck length(output) == length(x_query) || _throw_length_mismatch(length(x_query), length(output), "x_query", "output")
     x_eff, y_ext, bc_eff, extrap_eff = _periodic_extend_1d(x, y, bc, extrap)
 
-    Tdy = _promote_eltype(_coeff_op, eltype(x_eff), Tv)
+    # Value-matched width: dy buffer + slope arithmetic (incl. the `1 - tension`
+    # scale) run at `Tw` — see pchip_oneshot.jl.
+    Tw = _promote_grid_float(eltype(x_eff), Tv)
+    Tdy = _promote_eltype(_coeff_op, Tw, Tv)
     dy = acquire!(pool, Tdy, length(y_ext))
-    _cardinal_slopes!(dy, x_eff, y_ext, tension; bc = bc_eff)
+    _cardinal_slopes!(dy, x_eff, y_ext, tension, Tw; bc = bc_eff)
     searcher = _resolve_search(x_eff, x_query, search, hint)
     return _hermite_vector_loop!(output, x_eff, y_ext, dy, x_query, extrap_eff, deriv, searcher)
 end
@@ -144,7 +150,8 @@ Default `tension=0` is Catmull-Rom. C\$^1\$ continuous.
         search::AbstractSearchPolicy = AutoSearch(),
         hint::Union{Nothing, Base.RefValue{Int}} = nothing
     ) where {Tg, Tv, Tq <: Real}
-    x = _resolve_axis(x)
+    # Value-matched Tg: Int/OneTo grid + Float32 data → Float32 axis (tension follows).
+    x = _resolve_axis(x, _promote_grid_float(Tg, Tv))
     tension_f = float(eltype(x))(tension)
     extrap_eff = _resolve_extrap(extrap, bc, x, y)
     resolved = _resolve_coeffs(coeffs, x, xq)
@@ -172,7 +179,7 @@ In-place cardinal spline interpolation.
         search::AbstractSearchPolicy = AutoSearch(),
         hint::Union{Nothing, Base.RefValue{Int}} = nothing
     ) where {Tg, Tv, Tq <: Real}
-    x = _resolve_axis(x)
+    x = _resolve_axis(x, _promote_grid_float(Tg, Tv))
     tension_f = float(eltype(x))(tension)
     extrap_eff = _resolve_extrap(extrap, bc, x, y)
     resolved = _resolve_coeffs(coeffs, x, x_query)

src/cardinal/cardinal_slopes.jl

@@ -24,49 +24,60 @@ Compute cardinal spline slopes in-place.
 # Complexity
 O(n), single pass, zero allocation (writes into `dy`).
 """
+# Width-less form: delegate with the axis's own eltype (historic behavior; the
+# `float(one(Tw))` scale slot reproduces the old `one(Int) - tension::Float64`
+# promotion for Int axes). One-shot PreCompute backends pass the value-matched
+# `Tw` so slopes — including the `(1 - tension)` scale — are born at the value width.
+_cardinal_slopes!(dy::AbstractVector, x::AbstractVector, y::AbstractVector, tension; bc::AbstractBC = NoBC()) =
+    _cardinal_slopes!(dy, x, y, tension, eltype(x); bc)
+
 function _cardinal_slopes!(
         dy::AbstractVector,
-        x::AbstractVector{Tg},
+        x::AbstractVector,
         y::AbstractVector,
-        tension;
+        tension,
+        ::Type{Tw};
         bc::AbstractBC = NoBC()
-    ) where {Tg}
+    ) where {Tw}
     n = length(x)
     @assert n >= 2 "Cardinal spline requires at least 2 points"
     @assert length(y) == n "y length must match x"
     @assert length(dy) == n "dy length must match x"
 
-    scale = one(Tg) - tension
+    # Dimensionless scale at the value-matched width (`one(quantity)` strips
+    # units; `float` keeps Int axes on the historic Float64 arithmetic).
+    Tsc = typeof(float(one(Tw)))
+    scale = one(Tsc) - convert(Tsc, tension)
 
     # Special case: 2 points. PeriodicBC routes through the wrap-aware central
     # FD helper (see PCHIP n=2 note for rationale).
     if n == 2
         if bc isa PeriodicBC
-            @inbounds dy[1] = _cardinal_boundary_slope(x, y, 1, n, scale, bc)
-            @inbounds dy[2] = _cardinal_boundary_slope(x, y, 2, n, scale, bc)
+            @inbounds dy[1] = _cardinal_boundary_slope(Tw, x, y, 1, n, scale, bc)
+            @inbounds dy[2] = _cardinal_boundary_slope(Tw, x, y, 2, n, scale, bc)
             return dy
         end
         @inbounds begin
-            δ = _forward_secant(x, y, 1)
+            δ = _forward_secant(Tw, x, y, 1)
             dy[1] = scale * δ
             dy[2] = scale * δ
         end
         return dy
     end
 
     # Left endpoint: bc-dispatched helper.
-    @inbounds dy[1] = _cardinal_boundary_slope(x, y, 1, n, scale, bc)
+    @inbounds dy[1] = _cardinal_boundary_slope(Tw, x, y, 1, n, scale, bc)
 
     # Interior: central finite difference (K=3, no wrap needed).
     @inbounds for k in 2:(n - 1)
-        dy[k] = scale * _centered_secant(x, y, k)
+        dy[k] = scale * _centered_secant(Tw, x, y, k)
     end
 
     # Right endpoint: same bc-dispatched helper. PeriodicBC{:inclusive} yields
     # dy[n] == dy[1] automatically (closed-cycle symmetry); :exclusive yields
     # a different value using the seam secant — the helper handles both via
     # `_periodic_secant`/`_periodic_cell_width`.
-    @inbounds dy[n] = _cardinal_boundary_slope(x, y, n, n, scale, bc)
+    @inbounds dy[n] = _cardinal_boundary_slope(Tw, x, y, n, n, scale, bc)
 
     return dy
 end

src/constant/constant_oneshot.jl

@@ -213,8 +213,10 @@ vals = constant_interp(x, y, sorted_queries; search=LinearBinarySearch(linear_wi
     @boundscheck length(y) == length(x) || throw(ArgumentError("x and y must have same length"))
 
     # Surface-level BC-aware resolvers (zero-alloc reference wrapping). BC info
-    # lives in axis type after resolution → searcher uses `NoBC()`.
-    x_eff = _resolve_axis(x, bc)
+    # lives in axis type after resolution → searcher uses `NoBC()`. Raw Tg:
+    # the selection kernel keeps natural promotion (no float forcing) — an Int
+    # range stays `_CachedRange{Int}`, mirroring the ND constant scalar rule.
+    x_eff = _resolve_axis(x, bc, Tg)
     y_eff = _resolve_data(y, bc)
     extrap_eff = _resolve_extrap(extrap, bc, x_eff, y_eff)
     searcher = _resolve_search(x_eff, xi, search, hint)
@@ -266,8 +268,9 @@ function constant_interp!(
     @assert length(y) == length(x) "x and y must have same length"
     @assert length(output) == length(x_targets) "output must match x_targets length"
 
-    # Surface-level BC-aware resolvers (same template as Linear oneshot).
-    x_eff = _resolve_axis(x, bc)
+    # Surface-level BC-aware resolvers (same template as Linear oneshot);
+    # raw Tg — selection kernel keeps natural promotion (no float forcing).
+    x_eff = _resolve_axis(x, bc, eltype(x))
     y_eff = _resolve_data(y, bc)
     extrap_eff = _resolve_extrap(extrap, bc, x_eff, y_eff)
     searcher = _resolve_search(x_eff, x_targets, search, nothing)

src/constant/nd/constant_nd_adjoint_types.jl

@@ -74,7 +74,7 @@ struct ConstantAdjointND{
             Tg, N, B <: NTuple{N, AbstractBC},
             EP <: Tuple{Vararg{AbstractExtrap, N}}, SD <: Tuple{Vararg{AbstractSide, N}}, Tq,
         }
-        grids_c = map((g, bc) -> _convert_copy(_cache_axis(g, bc, Tg), Tg), grids, bcs)
+        grids_c = map((g, bc, T) -> _convert_copy(_cache_axis(g, bc, T), T), grids, bcs, ntuple(_ -> Tg, Val(N)))
         return new{Tg, N, typeof(grids_c), B, EP, SD, Tq}(grids_c, bcs, extraps, sides, anchors, grid_size)
     end
 end

src/constant/nd/constant_nd_oneshot.jl

@@ -29,7 +29,11 @@ function _constant_interp_nd_oneshot(
         ops::NTuple{N, AbstractEvalOp},
         hints = nothing
     ) where {Tv, N}
-    grids_eff = map(_resolve_axis, grids, bcs)
+    # Selection kernel: no x·y arithmetic → RAW grid eltype (no float forcing), mirroring
+    # `_nd_promote_grids_raw`/batch/persistent. All-Int stays Int; the output follows the
+    # natural promote_type(grid, data, query) — e.g. Int grid + Float32 query → Float32.
+    Tg = _promote_grid_eltype(grids)
+    grids_eff = _resolve_axes(grids, bcs, Tg)  # @generated static-Tg unroll (no Type-captured closure)
     # Bare GridIdx(k).val is NaN → resolve to the grid coordinate for the value kernel (search still uses .idx).
     query = map(_resolve_grididx, query, grids_eff)
     # Validate AND promote per axis: an in-domain NoExtrap axis becomes InBounds for the lean

src/constant/nd/constant_nd_types.jl

@@ -72,7 +72,7 @@ struct ConstantInterpolantND{
             bcs::NTuple{N, AbstractBC} = ntuple(_ -> NoBC(), Val(N)),
             store::StorePolicy = StorePolicy()
         ) where {Tg, Tv, N}
-        grids_c = map((g, bc) -> _own_or_ref_axis(_cache_axis(g, bc, Tg), Tg, store), grids, bcs)
+        grids_c = map((g, bc, T) -> _own_or_ref_axis(_cache_axis(g, bc, T), T, store), grids, bcs, ntuple(_ -> Tg, Val(N)))
         data_c = _own_or_ref_data(data, store)
         return new{Tg, Tv, N, typeof(grids_c), typeof(extraps), typeof(sides), typeof(searches), typeof(data_c)}(
             grids_c, data_c, extraps, sides, searches