(feat): duck-typed grid support for linear interpolation (#116)

* (feat): Linear 1D duck-typed grid support — Vector{ForwardDiff.Dual} (#81)

Enable LinearInterpolant to accept duck-typed grid scalars (e.g. ForwardDiff.Dual)
following the same Tv duck-typing pattern. Grid-side Dual propagates derivative
partials through spacing cache, binary search, and kernel evaluation.

Core changes:
- Add _PromotableGrid union (Integer, AbstractFloat, Rational) for fast-path gating
- Relax AbstractInterpolant, AbstractInterpolant1D type params (Tg unconstrained)
- Relax AbstractGridSpacing, ScalarSpacing, VectorSpacing (T unconstrained)
- Duck branch in _promote_itp_inputs: _PromotableGrid → Float path, else pass-through
- Relax scalar callable where-clause in interpolant_protocol.jl

Linear changes:
- LinearInterpolant struct: Tg unconstrained (duck grids get VectorSpacing{Tg})
- All 5 _linear_kernel overloads: drop Tg<:AbstractFloat
- Unified scalar one-shot entry: absorbs _promote_itp_inputs inline, removing the
  redundant Real→Float wrapper that caused infinite recursion on duck grids

Extension:
- _to_grid_type(::Real, ::Type{<:Dual}) override avoids Dual(Float) MethodError
  in _search_binary; returns primal for index comparison

Scope: Phase 1 — scalar eval + scalar one-shot only.
Vector queries (Phase 2), linear series (Phase 3), Range{Dual} (Phase 4),
adjoint, ND, and non-linear methods are out of scope.

Performance: zero regression on Float/Int/Range paths (12 benchmark cases,
byte-identical allocation counts vs master).

* (feat): Range{Dual} support via _CachedRange + DirectSearch O(1)

Extend the duck-typed grid infrastructure from the previous commit (Vector{Dual})
to also support Range{Dual} grids, preserving the O(1) DirectSearch fast path.

Core changes:
- Relax _CachedRange{T}, all constructors, _to_float overloads, and
  _to_float_adding_endpoint from T<:AbstractFloat to unconstrained T.
  x86 TwicePrecision specialization retains FT<:AbstractFloat and Dual
  ranges naturally bypass it via the generic _to_float fallback.
- _search_direct: extract primal before unsafe_trunc for index calculation;
  xL/xR retain full Dual type for kernel partial propagation.
- _create_spacing: relax AbstractRange{T<:AbstractFloat} to AbstractRange{T}.
  Range{Dual} gets ScalarSpacing{Dual} — O(1) cached h/inv_h.
- Simplify _promote_itp_inputs: remove _PromotableGrid branch entirely.
  All grid types (Int, Float, Dual, custom) now flow through a single path
  via _promote_grid_float → _to_float. Duck grids (Dual) are normalized to
  _CachedRange{Dual} just like Float ranges, with no y-value widening when
  Tg is not AbstractFloat.
- Remove unused _PromotableGrid constant and all references.
- Relax _to_float Vector identity/broadcast paths from FT<:AbstractFloat to T.

Pipeline for Range{Dual}:
  StepRangeLen{Dual} → _to_float → _CachedRange{Dual}
    → ScalarSpacing{Dual} → DirectSearch O(1) (primal-based trunc)
    → xL/xR as Dual → kernel partial propagation → Dual result

Tests:
- Range{Dual} construction, all 6 range source types (StepRangeLen, LinRange,
  UnitRange, StepRange, broadcast and scalar-mult variants)
- Range vs Vector path mid-interval parity (primal match, partial within 1e-10)
- Exact-knot behavior: analytic chain-rule verification + one-sided FD cross-check
  (DirectSearch→RIGHT, BinarySearch→LEFT — both match their respective one-sided FD)
- ForwardDiff.derivative MWE through Range grid
- TDG duck type: added float() protocol requirement

* (feat): complete Linear 1D Dual-grid support — vector, anchor, one-shot paths

Open all remaining Linear 1D paths to duck-typed grids (Vector{Dual}, Range{Dual}).

Protocol + shared:
- interpolant_protocol.jl: relax vector/in-place callable Tg constraint; include Tg
  in output eltype via promote_type(Tv, Tg, Tq)
- anchor_common.jl: relax _AnchorLoc{Tg} constraint
- utils.jl: relax 3-arg _promote_itp_inputs (query stays Float when Tg is duck);
  add generic _promote_for_anchor fallback for duck grids;
  relax _promote_value_type guard (skip y-widening when Tg is not AbstractFloat)

Linear:
- linear_interpolant.jl: relax _linear_vector_loop! and generic constructor
- linear_oneshot.jl: unify linear_interp! as single entry (absorbs promotion,
  removes redundant Real/Range wrappers that caused infinite recursion on duck
  grids); unify allocating vector one-shot similarly; relax _linear_interp_loop!
  to accept heterogeneous grid/query types
- linear_anchor.jl: relax _LinearAnchoredQuery and all overloads; add outer
  constructor that infers Tq from alpha's arithmetic type and promotes xq via
  convert — eliminates _promote_for_anchor calls from linear anchor path;
  merge two scalar _anchor_query overloads into one

Net effect: -54 lines (removed redundant wrappers, merged overloads).
All Linear 1D paths now accept duck-typed grids: constructor, scalar/vector
callable, scalar/vector one-shot, in-place, and anchor-based evaluation.

* (feat): Linear Series — complete Dual-grid support

Relax Tg<:AbstractFloat constraints on LinearSeriesInterpolant, all eval
overloads, oneshot series, and AbstractSeriesInterpolant abstract type.

Key changes beyond type constraint relaxation:
- _make_anchor: accept any xq type (was xq::Tg); _anchor_loc + outer constructor
  handle primal extraction and type widening internally
- Remove explicit _extract_primal + Tg(xq_primal) conversion from scalar eval
  caller — this caused MethodError on Dual grids (Dual(Float) constructor absent)
- _eval_linear_series_point!: use aq.xq instead of separate xq argument; aq.xq
  is already widened by the outer constructor to carry both query-side and grid-side
  Dual information
- Series output type: include Tg via promote_type(Tv, Tg) in _series_output_type
  so that Dual-grid results allocate the correct output vector

Both query-side Dual (ForwardDiff.derivative through series) and grid-side Dual
(Vector{Dual}/Range{Dual} grids) are now functional.

* (feat): Linear Adjoint — Dual-grid support

Relax Tg<:AbstractFloat on LinearAdjoint, AbstractAdjoint, AbstractAdjoint1D.

Key change: query normalization in linear_adjoint() avoids converting Float queries
to Dual (which would inject false grid-parameter partials). Queries stay at their
own float precision; the anchor outer constructor widens xq to match alpha's
arithmetic type during scatter.

Dot-product identity ⟨Wf, ȳ⟩ = ⟨f, Wᵀȳ⟩ verified with Dual grids.

* (test): add Series, Adjoint, and full-path coverage for Dual grids

Extend test/ext/test_linear_dual_grid.jl with:
- Linear Series: Vector{Dual} + Range{Dual} grids, primal parity with Float path
- Linear Adjoint: Vector{Dual} + Range{Dual} grids, dot-product identity
- All-paths coverage: scalar/vector callable, in-place, scalar/vector one-shot,
  in-place one-shot — for both Vector{Dual} and Range{Dual}
- Float in-place zero-alloc gate: itp(out, xq_vec) and linear_interp!(out, ...)
  both verified @Allocations == 0 after warmup

* (feat): Linear ND — complete Dual-grid support

Extend duck-typed grid support to all Linear ND paths: interpolant constructor,
scalar/batch eval (via interpolant_protocol.jl), one-shot scalar/batch, and
ND adjoint operator.

Changes:
- Relax AbstractInterpolantND, AbstractAdjointND abstract types
- Relax LinearInterpolantND, LinearAdjointND, _LinearNDAdjointAnchor structs
- Replace `Tg = Tg <: AbstractFloat ? Tg : Float64` promotion pattern with
  `Tg = float(Tg)` throughout (5 sites: interpolant, oneshot×3, adjoint×2)
- Add _value_type duck-Tg fallback: `where {T, Tg} = T` — when grid is Dual,
  y values are not promoted to grid type (no grid-parameter partials in data)
- Relax _linear_nd_oneshot_eval! where-clauses

Tests:
- 2D homogeneous Dual grids: primal parity + FD cross-check
- 2D heterogeneous axis (Dual×Float): per-axis FD cross-check
- 2D Range{Dual} + Vector{Float} mix: CachedRange verification + primal parity

* (fix): oneshot series Real wrapper infinite recursion on Julia 1.10 + stale docstrings

Bug fix:
- linear_oneshot_series.jl: remove 4 Real type promotion wrappers that caused
  infinite recursion on Julia 1.10 LTS (same pattern as the scalar/vector one-shot
  fix from earlier commits). Typed methods now handle promotion internally via
  _promote_grid_float + _to_float.
- Fix output type computation to use promoted Tg (eltype after _to_float), not
  the original raw Tg.
- Fix pool anchor type to use promoted Tq = promote_type(Tq, Tg_actual).

Docstring cleanup (from code review):
- cached_range.jl: docstring header now matches struct `_CachedRange{T}`
- linear_anchor.jl: Tg/Tq descriptions updated for duck-typed grids
- linear_series_interp.jl: _eval_linear_series_point! docstring updated
  (removed stale xq argument references)

* (fix): 4 Codex-identified issues — ND oneshot spacing, hinted search, series output types

TDD RED→GREEN for 4 identified issues:

P1-1: _create_spacing_pooled in nd_utils.jl had T<:AbstractFloat constraint,
  blocking ND one-shot paths for Dual grids. Fixed: T<:AbstractFloat → T.

P1-2: _search_direct! (hinted mutating variant) had T<:AbstractFloat on both
  overloads, blocking Range{Dual} grids with hint-based eval. Fixed: T → T.

P2-3: LinearSeriesInterpolant vector eval allocated Vector{Vector{Tv}} without
  including Tg in output type. Dual grid results couldn't be stored. Fixed:
  _series_output_type(Tv, Tq) → _series_output_type(promote_type(Tv, Tg), Tq).

P2-4: Oneshot series vector allocator used _value_type(Tv, Tg_p) which returns
  Tv unchanged for duck Tg. Fixed: use promote_type(_series_eltype(s), Tg_p)
  to include grid type in output allocation.

All 4 issues verified with dedicated test cases (RED→GREEN cycle).

* Runic formatting

* Runic formatting

* (fix): use _output_eltype instead of raw promote_type for output allocation

promote_type(MyDuck, Float64) returns Any on LTS (no promote_rule defined for
custom duck types), causing Vector{Any} output allocation. _output_eltype already
has the isconcretetype fallback: returns Tv when promote_type gives non-concrete.

Replaced promote_type(Tv, Tg, ...) with _output_eltype(Tv, Tg, ...) in:
- interpolant_protocol.jl: vector callable + ND batch allocator
- linear_oneshot.jl: vector one-shot allocator
- linear_series_interp.jl: scalar + vector series allocator
- linear_anchor.jl: anchor vector allocator
- linear_oneshot_series.jl: scalar + vector oneshot series allocator

Both duck-Tv (MyDuck) and duck-Tg (Dual) paths verified.

* (chore): PR review fixes — include test in runtests, update stale headers/docstrings

- Add test_linear_dual_grid.jl to runtests.jl default suite (Copilot review)
- Update ext test header: remove stale "Phase 1 scalar-only" scope note,
  reflect full coverage (1D/ND/series/adjoint/vector/anchor)
- Update LinearInterpolantND docstring: Tg description matches unconstrained code
- Runic formatting pass

* test: update allocation test to use @allocated with threshold

Author: mgyoo86

ext/FastInterpolationsForwardDiffExt.jl

@@ -14,7 +14,7 @@ module FastInterpolationsForwardDiffExt
     using FastInterpolations
     using ForwardDiff: Dual, value
 
-    import FastInterpolations: _extract_primal, _promote_for_anchor
+    import FastInterpolations: _extract_primal, _promote_for_anchor, _to_grid_type
 
     # ForwardDiff support: extract primal value from Dual for index search
     # - Use _extract_primal(xq) for comparisons and index lookup
@@ -27,5 +27,19 @@ module FastInterpolationsForwardDiffExt
     # This enables AD through anchor-based series evaluation with proper type consistency
     @inline _promote_for_anchor(xq::Dual{T, V, N}, ::Type{Tg}) where {T, V, N, Tg <: AbstractFloat} = xq + zero(Tg)
 
+    # ── Grid-side Dual support ─────────────────────────────────────────────
+    #
+    # When the grid is Vector{Dual}, binary search calls `_to_grid_type(xq, Tg)`
+    # with `Tg <: Dual`. The package-level fallback at `src/core/utils.jl:309`
+    # does `Tg(_extract_primal(xq))` which becomes `Dual(Float64)` — a MethodError
+    # because abstract `Dual` has no constructor from a single float.
+    #
+    # Fix: for a Dual-typed grid, return the primal of xq (a plain float). The
+    # subsequent `x[mid] <= xq` comparison inside `_search_binary` then compares
+    # `Dual <= Float`, which ForwardDiff forwards to the primal of `x[mid]` —
+    # correct and cheap. Arithmetic later (xq - xL, dL * inv_h) uses the original
+    # grid scalars and therefore carries the grid-side Dual partials through.
+    @inline _to_grid_type(xq::Real, ::Type{<:Dual}) = _extract_primal(xq)
+
 end # module
 0

src/core/abstract_types.jl

@@ -17,26 +17,29 @@
 # ========================================
 
 """
-    AbstractInterpolant{Tg<:AbstractFloat, Tv}
+    AbstractInterpolant{Tg, Tv}
 
 Abstract supertype for all interpolant objects.
 
 # Type Parameters
-- `Tg`: Grid/coordinate type (Float32 or Float64) - used for x-coordinates, spacing, search
+- `Tg`: Grid/coordinate type — normally `Float32`/`Float64`, but **unconstrained**
+        at this abstract level to allow duck-typed grid scalars (e.g. `ForwardDiff.Dual`)
+        that satisfy the grid arithmetic/ordering protocol (`-`, `inv`, `*`, `<`).
+        Non-linear concrete subtypes still enforce `Tg<:AbstractFloat` individually.
 - `Tv`: Value type (unconstrained) - used for y-values, coefficients, return values.
         Any type supporting the 5 core operations (+, -, Tg*Tv, Tv*Tg, Int*Tv).
 
 # Design Invariant
-- Grid operations (search, spacing) always use Tg
+- Grid operations (search, spacing) use Tg; duck grids use primal extraction for comparisons
 - Value operations (kernel, coefficients) use Tv
-- Evaluation returns type based on promote_type(Tv, query_type)
+- Evaluation returns type based on promote_type(Tv, Tg, query_type)
 
 # Subtypes
 - `AbstractInterpolant1D{Tg, Tv}`: 1D interpolants (shared callable protocol)
 - `AbstractSeriesInterpolant{Tg, Tv}`: Multi-series interpolants
 - `AbstractInterpolantND{Tg, Tv, N}`: N-dimensional interpolants
 """
-abstract type AbstractInterpolant{Tg <: AbstractFloat, Tv} end
+abstract type AbstractInterpolant{Tg, Tv} end
 
 """
     AbstractInterpolant1D{Tg<:AbstractFloat, Tv} <: AbstractInterpolant{Tg, Tv}
@@ -63,7 +66,7 @@ by implementing the required interface:
 - `CubicInterpolant{Tg, Tv}`: C2 cubic spline
 - `AbstractHermiteInterpolant1D{Tg, Tv}`: Cubic Hermite family (see subtypes below)
 """
-abstract type AbstractInterpolant1D{Tg <: AbstractFloat, Tv} <: AbstractInterpolant{Tg, Tv} end
+abstract type AbstractInterpolant1D{Tg, Tv} <: AbstractInterpolant{Tg, Tv} end
 
 """
     AbstractHermiteInterpolant1D{Tg, Tv} <: AbstractInterpolant1D{Tg, Tv}
@@ -119,7 +122,7 @@ sitp(output, 0.5)           # In-place evaluation
 This is a pure type hierarchy - no methods are defined on `AbstractSeriesInterpolant` itself.
 All functionality is implemented in concrete subtypes.
 """
-abstract type AbstractSeriesInterpolant{Tg <: AbstractFloat, Tv} <: AbstractInterpolant{Tg, Tv} end
+abstract type AbstractSeriesInterpolant{Tg, Tv} <: AbstractInterpolant{Tg, Tv} end
 
 """
     AbstractInterpolantND{Tg<:AbstractFloat, Tv, N}
@@ -150,7 +153,7 @@ itp((0.5, 0.5); deriv=(1, 0))      # ∂f/∂x
 gradient(itp, (0.5, 0.5))          # (∂f/∂x, ∂f/∂y)
 ```
 """
-abstract type AbstractInterpolantND{Tg <: AbstractFloat, Tv, N} <: AbstractInterpolant{Tg, Tv} end
+abstract type AbstractInterpolantND{Tg, Tv, N} <: AbstractInterpolant{Tg, Tv} end
 
 Base.size(itp::AbstractInterpolantND) = map(length, itp.grids)
 
@@ -184,7 +187,7 @@ Subtypes automatically inherit 6 callable overloads (Vector/Real/Tuple × alloc/
 - [`AbstractAdjoint1D`](@ref): 1D adjoint operators
 - [`AbstractAdjointND`](@ref): N-dimensional adjoint operators
 """
-abstract type AbstractAdjoint{Tg <: AbstractFloat} end
+abstract type AbstractAdjoint{Tg} end
 
 """
     AbstractAdjoint1D{Tg<:AbstractFloat} <: AbstractAdjoint{Tg}
@@ -206,7 +209,7 @@ Subtypes automatically inherit 6 callable overloads (Vector/Real/Tuple × alloc/
 - `LinearAdjoint{Tg, EP}`: Adjoint of linear interpolation (1D, pure scatter)
 - `CubicAdjoint{Tg, C, BC}`: Adjoint of cubic spline interpolation (1D)
 """
-abstract type AbstractAdjoint1D{Tg <: AbstractFloat} <: AbstractAdjoint{Tg} end
+abstract type AbstractAdjoint1D{Tg} <: AbstractAdjoint{Tg} end
 
 """
     AbstractAdjointND{Tg<:AbstractFloat, N} <: AbstractAdjoint{Tg}
@@ -225,7 +228,7 @@ Subtypes automatically inherit shared callable dispatch from `nd_adjoint_protoco
 - `CubicAdjointND{Tg, N, ...}`: Adjoint of cubic spline interpolation (ND)
 - `LinearAdjointND{Tg, N, ...}`: Adjoint of linear interpolation (ND)
 """
-abstract type AbstractAdjointND{Tg <: AbstractFloat, N} <: AbstractAdjoint{Tg} end
+abstract type AbstractAdjointND{Tg, N} <: AbstractAdjoint{Tg} end
 
 # ========================================
 # Type Helper Functions

src/core/anchor_common.jl

@@ -38,7 +38,7 @@ with NO geometry (h, inv_h, dL, dR). Geometry is each method's internal concern.
 - `xL::Tg`: left node x[idx]
 - `xR::Tg`: right node x[idx+1]
 """
-struct _AnchorLoc{Tg <: AbstractFloat, Tq <: Real}
+struct _AnchorLoc{Tg, Tq <: Real}
     idx::Int
     xq::Tq
     state::UInt8
@@ -74,7 +74,7 @@ Dual type. The interval search uses `_extract_primal(xq)` for comparisons.
         xq::Tq,
         wrap::Bool,
         policy::P = DEFAULT_SEARCHER
-    ) where {Tg <: AbstractFloat, Tq <: Real, P <: Searcher}
+    ) where {Tg, Tq <: Real, P <: Searcher}
     x_min, x_max = first(x), last(x)
 
     # Use primal value for comparisons (supports ForwardDiff.Dual)

src/core/cached_range.jl

@@ -9,7 +9,7 @@
 # Include order: grid_spacing.jl → cached_range.jl → search.jl → utils.jl
 
 """
-    _CachedRange{T <: AbstractFloat} <: AbstractRange{T}
+    _CachedRange{T} <: AbstractRange{T}
 
 `AbstractRange` subtype that pre-caches `first`, `last`, `step`, and `inv(step)` as
 plain `T` fields, enabling multiply-instead-of-divide search and avoiding
@@ -36,7 +36,7 @@ Because `_CachedRange <: AbstractRange`, all existing `AbstractRange` dispatch
 (DirectSearch routing, `_resolve_search_policy`, `search_interval`, etc.) propagates
 automatically without any changes at call sites.
 """
-struct _CachedRange{T <: AbstractFloat} <: AbstractRange{T}
+struct _CachedRange{T} <: AbstractRange{T}
     lo::T
     hi::T
     h::T
@@ -47,13 +47,13 @@ struct _CachedRange{T <: AbstractFloat} <: AbstractRange{T}
 end
 
 # Exact constructor: domain_lo == lo, domain_hi == hi (default for non-TwicePrecision paths)
-@inline function _CachedRange{T}(lo::T, hi::T, h::T, inv_h::T, len::Int) where {T <: AbstractFloat}
+@inline function _CachedRange{T}(lo::T, hi::T, h::T, inv_h::T, len::Int) where {T}
     return _CachedRange{T}(lo, hi, h, inv_h, len, lo, hi)
 end
 
-# Convenience: construct from any AbstractRange{T} where T <: AbstractFloat, using its own eltype.
+# Convenience: construct from any AbstractRange{T}, using its own eltype.
 # Internal code should prefer _to_float(x, Tg) when the desired target type Tg differs from eltype(x).
-_CachedRange(x::AbstractRange{T}) where {T <: AbstractFloat} = _to_float(x, T)
+_CachedRange(x::AbstractRange{T}) where {T} = _to_float(x, T)
 _CachedRange(x::_CachedRange) = x
 
 Base.length(r::_CachedRange) = r.len
@@ -78,7 +78,7 @@ end
 # Without this method, any code that windows or slices a `_CachedRange` (e.g. the
 # Hermite ND cell-local OnTheFly path in hetero_eval.jl) would silently degrade
 # its grid to a non-range type.
-@inline function Base.getindex(r::_CachedRange{T}, idx::AbstractUnitRange{<:Integer}) where {T <: AbstractFloat}
+@inline function Base.getindex(r::_CachedRange{T}, idx::AbstractUnitRange{<:Integer}) where {T}
     @boundscheck checkbounds(r, idx)
     new_len = length(idx)
     # Empty slice: return a length-0 _CachedRange anchored at r.lo (callers that
@@ -112,9 +112,9 @@ This is the universal normalizer for range grids. After `_to_float`, the grid ty
 space is exactly `{_CachedRange{FT}, Vector{FT}}`, eliminating downstream dispatch
 on `StepRangeLen`, `LinRange`, `OrdinalRange`, etc.
 """
-function _to_float(x::AbstractRange, ::Type{FT}) where {FT <: AbstractFloat}
-    h = FT(step(x))
-    return _CachedRange{FT}(FT(first(x)), FT(last(x)), h, inv(h), length(x))
+function _to_float(x::AbstractRange, ::Type{T}) where {T}
+    h = T(step(x))
+    return _CachedRange{T}(T(first(x)), T(last(x)), h, inv(h), length(x))
 end
 
 # x86_64: TwicePrecision first()/last() ~9ns each on Intel — bypass via plain-T muladd.
@@ -136,15 +136,15 @@ end
 end
 
 # _CachedRange same-type pass-through: already normalized, return as-is.
-_to_float(x::_CachedRange{T}, ::Type{T}) where {T <: AbstractFloat} = x
+_to_float(x::_CachedRange{T}, ::Type{T}) where {T} = x
 
 # _CachedRange type-mismatch (e.g. Float32 → Float64 via _convert_grid):
 # Uses 5-arg constructor (domain = exact). Any x86_64 domain widening from the source
 # is intentionally dropped: the type conversion itself introduces fresh rounding,
 # so re-widening would need to be based on the new FT, not the old T.
-function _to_float(x::_CachedRange, ::Type{FT}) where {FT <: AbstractFloat}
-    h = FT(x.h)
-    return _CachedRange{FT}(FT(x.lo), FT(x.hi), h, inv(h), x.len)
+function _to_float(x::_CachedRange, ::Type{T}) where {T}
+    h = T(x.h)
+    return _CachedRange{T}(T(x.lo), T(x.hi), h, inv(h), x.len)
 end
 
 """
@@ -159,25 +159,25 @@ Dispatch:
 """
 # domain_hi = hi_new (exact): the extension uses cached plain-T fields only
 # (no TwicePrecision involved), so no additional rounding uncertainty.
-@inline function _to_float_adding_endpoint(x::_CachedRange{FT}, ::Type{FT}) where {FT <: AbstractFloat}
+@inline function _to_float_adding_endpoint(x::_CachedRange{T}, ::Type{T}) where {T}
     hi_new = x.hi + x.h
-    return _CachedRange{FT}(
+    return _CachedRange{T}(
         x.lo, hi_new, x.h, x.inv_h, x.len + 1,
         x.domain_lo, hi_new
     )
 end
 
-@inline function _to_float_adding_endpoint(x::AbstractRange, ::Type{FT}) where {FT <: AbstractFloat}
-    r = _to_float(x, FT)
-    return _to_float_adding_endpoint(r, FT)
+@inline function _to_float_adding_endpoint(x::AbstractRange, ::Type{T}) where {T}
+    r = _to_float(x, T)
+    return _to_float_adding_endpoint(r, T)
 end
 
 # ========================================
 # _create_spacing: _CachedRange specialization
 # ========================================
 
 # _CachedRange already has h and inv_h cached — trivial field copy, no recomputation.
-function _create_spacing(x::_CachedRange{T}) where {T <: AbstractFloat}
+function _create_spacing(x::_CachedRange{T}) where {T}
     return ScalarSpacing{T}(x.h, x.inv_h)
 end
 

src/core/eval_ops.jl

@@ -150,6 +150,7 @@ end
 """Standard Julia numeric types that should be auto-promoted in convenience wrappers."""
 const _PromotableValue = Union{Integer, AbstractFloat, Rational, Complex}
 
+
 # ========================================
 #
 # Compile-time type tags for extrapolation mode selection.

src/core/grid_spacing.jl

@@ -16,9 +16,15 @@ Concrete subtypes:
 - `ScalarSpacing{T}`: Stores uniform spacing as scalars (O(1) memory)
 - `VectorSpacing{T}`: Stores non-uniform spacing as vectors (O(N) memory)
 
-The type parameter `T` is the floating-point type (Float32 or Float64).
+The type parameter `T` is normally a floating-point type (Float32 or Float64),
+but it is **unconstrained** at the abstract level to support duck-typed grids
+whose scalar type satisfies the required arithmetic protocol (`-`, `inv`, `*`)
+— e.g. `ForwardDiff.Dual`. Those grids land in `VectorSpacing{Dual}` via the
+unconstrained `_create_spacing(::AbstractVector)` method below, which uses only
+`Vector{T}(undef, ...)`, subtraction, and `inv` — all of which Dual supports.
+The Float path is unchanged: concrete allocations still specialize per-eltype.
 """
-abstract type AbstractGridSpacing{T <: AbstractFloat} end
+abstract type AbstractGridSpacing{T} end
 
 """
     ScalarSpacing{T} <: AbstractGridSpacing{T}
@@ -43,12 +49,14 @@ x = range(0.0, 1.0, 1001)  # 1000 intervals, uniform spacing
 spacing = _create_spacing(x)  # ScalarSpacing{Float64}(0.001, 1000.0)
 ```
 """
-struct ScalarSpacing{T <: AbstractFloat} <: AbstractGridSpacing{T}
+struct ScalarSpacing{T} <: AbstractGridSpacing{T}
     h::T
     inv_h::T
 end
 
-# Outer constructor for type promotion with mixed Real types
+# Outer constructor for type promotion with mixed Real types.
+# Standard numerics coerce through Float64; duck types (Dual, Measurement, etc.)
+# are never routed here — they reach VectorSpacing via _create_spacing directly.
 function ScalarSpacing(h::Real, inv_h::Real)
     T = promote_type(typeof(h), typeof(inv_h))
     T = T <: AbstractFloat ? T : Float64
@@ -66,15 +74,24 @@ where intervals may have different spacings.
 - `inv_h::Vector{T}`: Precomputed reciprocals inv_h[i] = 1/h[i]
 
 # Memory
-O(N) - stores 2*(n-1) floating-point values.
+O(N) - stores 2*(n-1) values.
+
+# Type parameter
+`T` is normally a float (Float32/Float64) but may also be any scalar type supporting
+`-`, `inv`, and `Vector{T}(undef, n)` — e.g. `ForwardDiff.Dual`. Duck-typed grids
+carry their derivative partials through `h`/`inv_h`, so the cached spacing is
+differentiable and reusable across queries.
 
 # Example
 ```julia
-x = [0.0, 0.3, 0.7, 1.0]  # Non-uniform spacing
-spacing = _create_spacing(x)  # VectorSpacing with h=[0.3, 0.4, 0.3]
+x = [0.0, 0.3, 0.7, 1.0]                 # Non-uniform Float spacing
+spacing = _create_spacing(x)              # VectorSpacing{Float64} with h=[0.3, 0.4, 0.3]
+
+x_dual = ForwardDiff.Dual{:tag}.(x, ...)  # Dual-valued grid
+spacing_d = _create_spacing(x_dual)       # VectorSpacing{Dual{:tag,Float64,1}}
 ```
 """
-struct VectorSpacing{T <: AbstractFloat} <: AbstractGridSpacing{T}
+struct VectorSpacing{T} <: AbstractGridSpacing{T}
     h::Vector{T}
     inv_h::Vector{T}
 end
@@ -133,7 +150,7 @@ Extracts the constant step size and precomputes its reciprocal.
 Defensive fallback for non-normalized Range inputs; the primary path
 uses `_create_spacing(::_CachedRange)` in `cached_range.jl`.
 """
-function _create_spacing(x::AbstractRange{T}) where {T <: AbstractFloat}
+function _create_spacing(x::AbstractRange{T}) where {T}
     # step(x) already returns T for AbstractRange{T}, avoid redundant conversion
     h = step(x)
     inv_h = inv(h)
@@ -153,7 +170,7 @@ Create vector spacing for non-uniform grids (Vector inputs).
 
 Computes h[i] = x[i+1] - x[i] and inv_h[i] = 1/h[i] for each interval.
 """
-function _create_spacing(x::AbstractVector{T}) where {T <: AbstractFloat}
+function _create_spacing(x::AbstractVector{T}) where {T}
     n = length(x)
     h = Vector{T}(undef, n - 1)
     inv_h = Vector{T}(undef, n - 1)

src/core/interpolant_protocol.jl

@@ -35,7 +35,7 @@
         deriv::DerivOp = EvalValue(),
         search::AbstractSearchPolicy = _itp_search(itp),
         hint::Union{Nothing, Base.RefValue{Int}} = nothing
-    ) where {Tg <: AbstractFloat, Tv}
+    ) where {Tg, Tv}
     grid = _itp_grid(itp)
     extrap = _itp_extrap(itp)
     @boundscheck _check_domain(grid, xq, extrap)
@@ -53,10 +53,12 @@ function (itp::AbstractInterpolant1D{Tg, Tv})(
         deriv::DerivOp = EvalValue(),
         search::AbstractSearchPolicy = _itp_search(itp),
         hint::Union{Nothing, Base.RefValue{Int}} = nothing
-    ) where {Tg <: AbstractFloat, Tv, Tq <: Real}
+    ) where {Tg, Tv, Tq <: Real}
     grid = _itp_grid(itp)
     extrap = _itp_extrap(itp)
-    T_out = promote_type(Tv, Tq)
+    # Include Tg in promotion: when the grid is Dual, interpolation weights
+    # carry grid partials, so the output type must reflect Tg arithmetic.
+    T_out = _output_eltype(Tv, Tg, Tq)
     output = Vector{T_out}(undef, length(xq))
     searcher = _resolve_search(grid, xq, search, hint)
     _itp_vector_loop!(output, itp, xq, extrap, deriv, searcher)
@@ -73,7 +75,7 @@ function (itp::AbstractInterpolant1D{Tg, Tv})(
         deriv::DerivOp = EvalValue(),
         search::AbstractSearchPolicy = _itp_search(itp),
         hint::Union{Nothing, Base.RefValue{Int}} = nothing
-    ) where {Tg <: AbstractFloat, Tv, Tq <: Real}
+    ) where {Tg, Tv, Tq <: Real}
     @assert length(output) == length(xq) "output length must match xq length"
     grid = _itp_grid(itp)
     extrap = _itp_extrap(itp)
@@ -191,6 +193,6 @@ function (itp::AbstractInterpolantND{Tg, Tv, N})(
         hint::Union{Nothing, NTuple{N, Base.RefValue{Int}}} = nothing
     ) where {Tg, Tv, N}
     Tq = _query_eltype(queries)
-    output = Vector{promote_type(Tv, Tg, Tq)}(undef, _query_length(queries))
+    output = Vector{_output_eltype(Tv, Tg, Tq)}(undef, _query_length(queries))
     return itp(output, queries; deriv = deriv, search = search, hint = hint)
 end

src/core/nd_utils.jl

@@ -864,7 +864,7 @@ end
 Pool-aware spacing for Range grids. Delegates to `_create_spacing` since
 ScalarSpacing is already zero-allocation (two scalar values).
 """
-@inline _create_spacing_pooled(pool::AbstractArrayPool, x::AbstractRange{T}) where {T <: AbstractFloat} = _create_spacing(x)
+@inline _create_spacing_pooled(pool::AbstractArrayPool, x::AbstractRange{T}) where {T} = _create_spacing(x)
 
 """
     _create_spacing_pooled(pool, x::AbstractVector{T}) -> VectorSpacing{T}
@@ -873,7 +873,7 @@ Pool-aware spacing for Vector grids. Acquires `h` and `inv_h` arrays
 from the pool instead of heap-allocating. The pool buffers are released
 automatically when the enclosing `@with_pool` scope exits.
 """
-@inline function _create_spacing_pooled(pool::AbstractArrayPool, x::AbstractVector{T}) where {T <: AbstractFloat}
+@inline function _create_spacing_pooled(pool::AbstractArrayPool, x::AbstractVector{T}) where {T}
     n = length(x)
     h = acquire!(pool, T, n - 1)
     inv_h = acquire!(pool, T, n - 1)