(fix+refac): Cell-local NaN + Tq carrier through deriv-zero paths

src/constant/constant_anchor.jl

@@ -274,13 +274,13 @@ end
 # Used by both interpolant anchor dispatch AND series one-shot evaluation.
 
 # Default case (extension, wrap, inbounds): kernel with right-boundary check.
-# Short-circuits use `* one(aq.xq)` to match the kernel's `* one(dL)` carrier
-# propagation — without it, Int y + Float xq returns Union{Int,Float}.
+# Right-edge short-circuits thread both Tq (`one(aq.xq)`) and Tg (`one(x_last)`)
+# carriers and use cell-local `y[aq.idxR]` so NaN at the queried cell propagates.
 @inline function _constant_eval_at_anchor(
         y::AbstractVector, x_last, aq::_ConstantAnchoredQuery,
         op::AbstractEvalOp, side_param::AbstractSide, ::AbstractExtrap
     )
-    aq.xq == x_last && return (op isa EvalValue ? (@inbounds y[aq.idxR] * one(aq.xq)) : 0 * first(y) * one(aq.xq))
+    aq.xq == x_last && return (op isa EvalValue ? (@inbounds y[aq.idxR] * one(aq.xq) * one(x_last)) : @inbounds 0 * y[aq.idxR] * one(aq.xq) * one(x_last))
     @inbounds return _constant_kernel(op, y[aq.idxL], y[aq.idxR], aq.h, aq.dL, side_param)
 end
 
@@ -290,7 +290,7 @@ end
         op::AbstractEvalOp, side_param::AbstractSide, ::NoExtrap
     )
     aq.state != IN_DOMAIN && throw(DomainError(aq.xq, "query point outside domain"))
-    aq.xq == x_last && return (op isa EvalValue ? (@inbounds y[aq.idxR] * one(aq.xq)) : 0 * first(y) * one(aq.xq))
+    aq.xq == x_last && return (op isa EvalValue ? (@inbounds y[aq.idxR] * one(aq.xq) * one(x_last)) : @inbounds 0 * y[aq.idxR] * one(aq.xq) * one(x_last))
     @inbounds return _constant_kernel(op, y[aq.idxL], y[aq.idxR], aq.h, aq.dL, side_param)
 end
 
@@ -303,7 +303,7 @@ end
         y_bnd = aq.state == OOB_LEFT ? first(y) : last(y)
         return _eval_extrapolation(op, y_bnd, extrap, aq.xq)
     end
-    aq.xq == x_last && return (op isa EvalValue ? (@inbounds y[aq.idxR] * one(aq.xq)) : 0 * first(y) * one(aq.xq))
+    aq.xq == x_last && return (op isa EvalValue ? (@inbounds y[aq.idxR] * one(aq.xq) * one(x_last)) : @inbounds 0 * y[aq.idxR] * one(aq.xq) * one(x_last))
     @inbounds return _constant_kernel(op, y[aq.idxL], y[aq.idxR], aq.h, aq.dL, side_param)
 end
 
@@ -319,21 +319,16 @@ end
 # Interpolant Anchor Dispatch (thin wrappers)
 # ========================================
 
-# No extrapolation: enriched error message with domain bounds
+# No extrapolation: enriched error message with domain bounds, then delegate
+# to the shared cell-local path (right-edge & in-cell carrier handling).
 @inline function _constant_anchor_dispatch(
-        itp::ConstantInterpolant{T},
-        aq::_ConstantAnchoredQuery{T},
-        op::O,
-        ::NoExtrap
-    ) where {T, O <: AbstractEvalOp}
+        itp::ConstantInterpolant, aq::_ConstantAnchoredQuery, op::AbstractEvalOp, ::NoExtrap
+    )
     if aq.state != IN_DOMAIN
         x_min, x_max = first(itp.x), last(itp.x)
         throw(DomainError(aq.xq, "query point outside domain [$x_min, $x_max]"))
     end
-    if aq.xq == last(itp.x)
-        return op isa EvalValue ? (@inbounds itp.y[aq.idxR] * one(aq.xq)) : zero(T) * one(aq.xq)
-    end
-    @inbounds return _constant_kernel(op, itp.y[aq.idxL], itp.y[aq.idxR], aq.h, aq.dL, itp.side)
+    return _constant_eval_at_anchor(itp.y, last(itp.x), aq, op, itp.side, NoExtrap())
 end
 
 # Inside domain or extension mode: delegate to shared

src/constant/constant_oneshot.jl

@@ -30,10 +30,13 @@
         searcher::S
     ) where {Tg, Tv, Tq <: Real, S <: Searcher}
     if _extract_primal(xi) == _extract_primal(last(x))
-        # `last(y)` covers both raw vectors and `_ExclusivePeriodicData` (cyclic
-        # `inner[1]`). `* one(xi)` propagates Tq carrier to match the kernel
-        # path below (without it, Int y + Float xq returns Union{Int,Float}).
-        return op isa EvalValue ? last(y) * one(xi) : 0 * first(y) * one(xi)
+        # `last(y)` for both raw vectors and `_ExclusivePeriodicData` (cyclic
+        # `inner[1]`). `one(Tq) * one(Tg)` threads both query and grid carriers
+        # — must match the kernel branch (which threads Tg via `dL`), else
+        # inference becomes `Union{Tv, Dual}` when grid is Dual and query is Float.
+        return op isa EvalValue ?
+            last(y) * one(Tq) * one(Tg) :
+            0 * last(y) * one(Tq) * one(Tg)
     end
     idx, idx_R, xL, xR = search_interval(searcher, x, xi)
     dL = xi - xL
@@ -106,7 +109,9 @@ end
     # the exact boundary. `last(_ExclusivePeriodicData) = inner[1]` so `:exclusive`
     # cyclic wrap is preserved; raw Vector yields `y[n]`.
     _extract_primal(xi_wrapped) == _extract_primal(last(x)) &&
-        return op isa EvalValue ? last(y) * one(xi_wrapped) : 0 * first(y) * one(xi_wrapped)
+        return op isa EvalValue ?
+        last(y) * one(Tq) * one(Tg) :
+        0 * last(y) * one(Tq) * one(Tg)
     idx, idx_R, xL, xR = search_interval(searcher, x, xi_wrapped)
     dL = xi_wrapped - xL
     # Unwrap data once: `search_interval` already resolved the seam (idx_R = 1
@@ -148,7 +153,7 @@ Constant (step/piecewise constant) interpolation at a single point.
   - `NearestSide()` (default): nearest neighbor (left tie-breaking at midpoint)
   - `LeftSide()`: always use left value
   - `RightSide()`: use right value (except at grid points)
-- `deriv::DerivOp`: Derivative order (`EvalValue()`, `DerivOp(1)`, or `DerivOp(2)`). Derivatives are always 0.
+- `deriv::DerivOp`: Derivative order — `EvalValue()` for value, or any `DerivOp(n)` for nth derivative (constant interpolation has zero derivative at all orders ≥ 1).
 - `search::AbstractSearchPolicy`: Search algorithm for interval finding
   - `BinarySearch()` (default): O(log n) binary search, stateless
   - `LinearBinarySearch(linear_window=0)`: O(1) if hint valid, O(log n) fallback

src/constant/constant_oneshot_series.jl

@@ -285,7 +285,7 @@ function constant_interp(
     K = n_series(s)
     Tv = _series_eltype(s)
     Tv_out = _output_eltype(_constant_kernel_shape, Tg, Tv, Tq)
-    outputs = [Vector{Tv_out}(undef, length(xqs)) for _ in 1:K]
+    outputs = _alloc_series_batch_outputs(Tv_out, K, length(xqs))
     constant_interp!(outputs, x, s, xqs; bc, side, extrap, deriv, search)
     return outputs
 end

src/constant/constant_series_interp.jl

@@ -3,8 +3,7 @@
 # ║         Multiple y-data series sharing the same x-grid                    ║
 # ╚═══════════════════════════════════════════════════════════════════════════╝
 #
-# Phase E.2: Unified matrix storage for optimal performance.
-# Key optimization: Adaptive layout with lazy transpose for scalar queries.
+# Unified matrix storage with adaptive layout — lazy transpose for scalar queries.
 #
 # Include order: ... → constant_anchor.jl → constant_series_interp.jl
 #
@@ -211,18 +210,16 @@ Outside-domain delegates to `_eval_series_at_anchor!` for extrapolation.
     y_point = _ensure_point_layout!(sitp)
     n_pts = n_points(sitp)
 
-    # Special case: at right boundary (use primal for comparison)
+    # Right boundary: yL == yR == y_point[k, n_pts]. `_constant_kernel(op, ...)`
+    # produces `y_at_max * one(dL)` for EvalValue, `0 * y_at_max * one(dL)` for
+    # deriv — cell-local NaN and Tq carrier propagate through both.
     xq_primal = _extract_primal(xq)
     if xq_primal == _extract_primal(last(sitp.x))
-        if op isa EvalValue
-            @inbounds @simd for k in axes(output, 1)
-                output[k] = y_point[k, n_pts]
-            end
-        else
-            z = 0 * first(y_point)
-            @inbounds @simd for k in axes(output, 1)
-                output[k] = z
-            end
+        h = aq.h
+        dL = aq.dL
+        @inbounds @simd for k in axes(output, 1)
+            y_at_max = y_point[k, n_pts]
+            output[k] = _constant_kernel(op, y_at_max, y_at_max, h, dL, sitp.side)
         end
         return output
     end
@@ -273,13 +270,13 @@ end
         n_pts::Int,
         ::Tg,
         ::Tg,
-        ::_ConstantAnchoredQuery{Tg},
+        aq::_ConstantAnchoredQuery{Tg},
         extrap::_ClampOrFill,
         ::AbstractSide,
         op::AbstractEvalOp,
         side::UInt8
     ) where {Tg, Tv}
-    return _fill_constant_extrap_simd!(out, y_point, side, n_pts, op, extrap)
+    return _fill_constant_extrap_simd!(out, y_point, side, n_pts, op, extrap, aq)
 end
 
 # ExtendExtrap - extend using same constant value at boundary interval
@@ -588,13 +585,12 @@ Internal: Evaluate single series at single query point with extrapolation handli
         side_val::AbstractSide,
         op::AbstractEvalOp
     ) where {Tg, Tv}
-    # Special case: at right boundary (MUST be preserved!)
+    # Right boundary: yL == yR == y[n_pts, k]. Kernel handles op-dispatch
+    # (value = y_at_max * one(dL), deriv = 0 * y_at_max * one(dL)) so the
+    # cell-local NaN and Tq carrier both propagate.
     if aq.xq == x_max
-        if op isa EvalValue
-            @inbounds return y[n_pts, k]
-        else
-            return 0 * first(y)  # Derivatives of step function are zero
-        end
+        @inbounds y_at_max = y[n_pts, k]
+        return _constant_kernel(op, y_at_max, y_at_max, aq.h, aq.dL, side_val)
     end
 
     # Inside domain: normal evaluation
@@ -606,14 +602,16 @@ Internal: Evaluate single series at single query point with extrapolation handli
     if extrap isa ExtendExtrap || extrap isa WrapExtrap
         return _eval_constant_series_anchored(y, k, aq, side_val, op)
     elseif extrap isa _ClampOrFill
-        return _constant_extrap_boundary_value(y, aq.state, n_pts, k, op, extrap)
+        return _constant_extrap_boundary_value(y, aq.state, n_pts, k, op, extrap, aq)
     else
         _throw_extrap_domain_error(aq.xq, x_min, x_max)
     end
 end
 
 """
 Internal: Core constant evaluation for series k at anchored query point.
+Value and deriv share `_constant_kernel(op, ...)` — deriv returns
+`0 * y_left * one(dL)` so cell-local NaN and Tq carrier propagate.
 """
 @inline function _eval_constant_series_anchored(
         y::Matrix{Tv},
@@ -622,16 +620,9 @@ Internal: Core constant evaluation for series k at anchored query point.
         side_val::AbstractSide,
         op::AbstractEvalOp
     ) where {Tg, Tv}
-    # Derivatives of constant (step) function are zero
-    if !(op isa EvalValue)
-        return 0 * first(y)
-    end
-
-    idxL = aq.idxL
-    idxR = aq.idxR
     @inbounds begin
-        y_left = y[idxL, k]
-        y_right = y[idxR, k]
+        y_left = y[aq.idxL, k]
+        y_right = y[aq.idxR, k]
     end
-    return _constant_kernel(EvalValue(), y_left, y_right, aq.h, aq.dL, side_val)
+    return _constant_kernel(op, y_left, y_right, aq.h, aq.dL, side_val)
 end

src/constant/nd/constant_nd_eval.jl

@@ -26,14 +26,13 @@ end
 # In-place + allocating batch use the inherited `AbstractInterpolantND`
 # protocol (trait-sized allocator via `_output_eltype`). Scalar routes
 # through `_eval_nd_at_point`; Constant's "any derivative → 0" rule is
-# wired via `_deriv_zero_fill` below.
+# applied inside `_eval_at_cell` (see below) via `_constant_nd_evaluate`'s
+# multi-dispatch — the cell-local kernel result is multiplied by `0` when
+# any deriv operator is non-EvalValue, so NaN/Inf in the queried cell's
+# data propagates through IEEE `NaN * 0 = NaN`.
 
-# Derivative zero-fill trait: constant has zero derivative at all orders
-@inline _deriv_zero_fill(::ConstantInterpolantND, ops::NTuple{N, AbstractEvalOp}, ::Val{N}) where {N} =
-    _has_any_derivative(ops, Val(N))
-
-# Zero-ref for fill-value derivative computation (duck-typed zero via 0 * data_element)
-@inline _zero_ref(itp::ConstantInterpolantND) = @inbounds first(itp.data)
+# Per-method sample of `Tv` for fill-value paths (e.g. `_try_fill_oob`).
+@inline _sample_data(itp::ConstantInterpolantND) = @inbounds first(itp.data)
 
 # ========================================
 # CELL LOCATION (locate once, evaluate many)
@@ -94,13 +93,27 @@ end
         ops::NTuple{N, AbstractEvalOp}
     ) where {Tg, Tv, N}
     data, stencils, hs, sides, q_eval, Ls = cell
-    if _has_any_derivative(ops, Val(N))
-        # `prod(one, q_eval)` folds carrier over every axis.
-        return 0 * first(itp.data) * prod(one, q_eval)
-    end
-    return _constant_nd_kernel(data, stencils, hs, sides, q_eval, Ls)
+    return _constant_nd_evaluate(data, stencils, hs, sides, q_eval, Ls, ops, Val(N))
 end
 
+# Deriv-aware Constant ND evaluation via two-method dispatch (mirrors Linear's
+# `_linear_weight` pattern). `_constant_nd_kernel` only handles the EvalValue
+# path; the deriv fallback multiplies the cell-local kernel result by `0`.
+# The cell-local NaN/Inf carrier survives `* 0` via IEEE (`NaN * 0 = NaN`),
+# so NaN data in the queried cell propagates through value and partials slots.
+# Intentionally NOT a `fill!`/`zero(Tv)` shortcut: NaN propagation, Dual
+# carrier, and duck-typed Tq all require running the kernel — the cost
+# matches the value path.
+@inline _constant_nd_evaluate(
+    data, stencils, hs, sides, q_eval, Ls,
+    ::NTuple{N, EvalValue}, ::Val{N}
+) where {N} = _constant_nd_kernel(data, stencils, hs, sides, q_eval, Ls)
+
+@inline _constant_nd_evaluate(
+    data, stencils, hs, sides, q_eval, Ls,
+    ::NTuple{N, AbstractEvalOp}, ::Val{N}
+) where {N} = _constant_nd_kernel(data, stencils, hs, sides, q_eval, Ls) * 0
+
 # ========================================
 # Derivative Check
 # ========================================