ProjectTorreyPines · mgyoo86 · May 19, 2026 · May 18, 2026 · May 18, 2026 · May 18, 2026
diff --git a/docs/src/boundary-conditions/overview.md b/docs/src/boundary-conditions/overview.md
@@ -3,7 +3,7 @@
 Boundary conditions (BCs) determine how spline interpolants behave at the domain endpoints. They are essential for **quadratic** and **cubic** splines, which require constraints to uniquely determine the spline coefficients.
 
 !!! note "Not Needed for All Methods"
-    - **Constant** and **Linear** interpolation do not require boundary conditions
+    - **Constant** and **Linear** interpolation do not require boundary conditions for construction, but accept `PeriodicBC` as a periodic domain convention
     - **Quadratic** splines require one BC (single endpoint)
     - **Cubic** splines require two BCs (both endpoints)
 
@@ -22,7 +22,7 @@ AbstractBC{T}
 │       ├── QuadraticFit     # = PolyFit{2} (3 points, O(h²))
 │       └── CubicFit        # = PolyFit{3} (4 points, O(h³))
 ├── BCPair{T,L,R}           # Both endpoints (cubic only)
-├── PeriodicBC{T}           # Periodic BC (cubic only)
+├── PeriodicBC{T}           # Periodic axis convention / method-specific BC
 ├── ZeroCurvBC{T}            # Zero curvature at both ends (cubic only)
 ├── ZeroSlopeBC{T}            # Zero slope at both ends (cubic only)
 ├── MinCurvFit{T}           # Minimum curvature (quadratic only)
@@ -118,6 +118,6 @@ PeriodicBC()     # S(x) = S(x + τ) with C² continuity
 ## See Also
 
 - [PointBC Details](pointbc.md) — In-depth explanation of PolyFit with visualizations
-- [PeriodicBC Details](periodicbc.md) — Inclusive vs exclusive endpoints, period inference, comparison with `WrapExtrap()`
+- [PeriodicBC Details](periodicbc.md) — `PeriodicBC` vs `WrapExtrap()`, endpoint policy, per-method behavior (Constant/Linear, PCHIP/Cardinal/Akima, Cubic)
 - [Quadratic Interpolation](../interpolation/quadratic.md) — BC examples in context
 - [Cubic Interpolation](../interpolation/cubic.md) — BCPair and PeriodicBC details
diff --git a/docs/src/boundary-conditions/periodicbc.md b/docs/src/boundary-conditions/periodicbc.md
diff --git a/src/akima/akima_oneshot.jl b/src/akima/akima_oneshot.jl
@@ -27,7 +27,7 @@
     Tdy = _output_eltype(Tv, float(eltype(x_eff)))
     dy = acquire!(pool, Tdy, length(y_ext))
     _akima_slopes!(dy, x_eff, y_ext; bc = bc_eff)
-    searcher = _resolve_search(x_eff, xq, search, hint, 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
 
@@ -50,7 +50,7 @@ end
     Tdy = _output_eltype(Tv, float(eltype(x_eff)))
     dy = acquire!(pool, Tdy, length(y_ext))
     _akima_slopes!(dy, x_eff, y_ext; bc = bc_eff)
-    searcher = _resolve_search(x_eff, x_query, search, hint, 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
 
@@ -74,7 +74,7 @@ end
     x_eff = _resolve_axis(x, bc)
     y_eff = _resolve_data(y, bc)
 
-    searcher = _resolve_search(x_eff, xq, search, hint, NoBC())
+    searcher = _resolve_search(x_eff, xq, search, hint)
     return _hermite_eval_at_point(x_eff, y_eff, AkimaSlopes(bc), xq, extrap, deriv, searcher)
 end
 
@@ -97,7 +97,7 @@ end
     y_eff = _resolve_data(y, bc)
 
 
-    searcher = _resolve_search(x_eff, x_query, search, hint, NoBC())
+    searcher = _resolve_search(x_eff, x_query, search, hint)
     return _hermite_vector_loop!(output, x_eff, y_eff, AkimaSlopes(bc), x_query, extrap, deriv, searcher)
 end
 

diff --git a/src/cardinal/cardinal_oneshot.jl b/src/cardinal/cardinal_oneshot.jl
@@ -33,7 +33,7 @@
     Tdy = _output_eltype(Tv, float(eltype(x_eff)))
     dy = acquire!(pool, Tdy, length(y_ext))
     _cardinal_slopes!(dy, x_eff, y_ext, tension; bc = bc_eff)
-    searcher = _resolve_search(x_eff, xq, search, hint, 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
 
@@ -57,7 +57,7 @@ end
     Tdy = _output_eltype(Tv, float(eltype(x_eff)))
     dy = acquire!(pool, Tdy, length(y_ext))
     _cardinal_slopes!(dy, x_eff, y_ext, tension; bc = bc_eff)
-    searcher = _resolve_search(x_eff, x_query, search, hint, 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
 
@@ -92,7 +92,7 @@ end
     # `x[n+1]`, so `Base.getindex` raw passthrough on the wrapper is safe.
     x_eff = _resolve_axis(x, bc)
     y_eff = _resolve_data(y, bc)
-    searcher = _resolve_search(x_eff, xq, search, hint, NoBC())
+    searcher = _resolve_search(x_eff, xq, search, hint)
     return _hermite_eval_at_point(x_eff, y_eff, CardinalSlopes(tension, bc), xq, extrap, deriv, searcher)
 end
 
@@ -114,7 +114,7 @@ end
     @boundscheck length(output) == length(x_query) || _throw_length_mismatch(length(x_query), length(output), "x_query", "output")
     x_eff = _resolve_axis(x, bc)
     y_eff = _resolve_data(y, bc)
-    searcher = _resolve_search(x_eff, x_query, search, hint, NoBC())
+    searcher = _resolve_search(x_eff, x_query, search, hint)
     return _hermite_vector_loop!(output, x_eff, y_eff, CardinalSlopes(tension, bc), x_query, extrap, deriv, searcher)
 end
 

diff --git a/src/constant/constant_oneshot.jl b/src/constant/constant_oneshot.jl
@@ -197,7 +197,7 @@ vals = constant_interp(x, y, sorted_queries; search=LinearBinarySearch(linear_wi
     x_eff = _resolve_axis(x, bc)
     y_eff = _resolve_data(y, bc)
     extrap_eff = _resolve_extrap(extrap, bc, x_eff, y_eff)
-    searcher = _resolve_search(x_eff, xi, search, hint, NoBC())
+    searcher = _resolve_search(x_eff, xi, search, hint)
     return _constant_eval_at_point(x_eff, y_eff, xi, extrap_eff, side, deriv, searcher)
 end
 
@@ -250,7 +250,7 @@ function constant_interp!(
     x_eff = _resolve_axis(x, bc)
     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, NoBC())
+    searcher = _resolve_search(x_eff, x_targets, search, nothing)
     _constant_vector_loop!(output, x_eff, y_eff, x_targets, extrap_eff, side, deriv, searcher)
     return output
 end

diff --git a/src/constant/constant_oneshot_series.jl b/src/constant/constant_oneshot_series.jl
@@ -83,7 +83,9 @@ end
     Tv_out = Tq <: _PromotableValue ? Tv : promote_type(Tv, Tq)
     output = Vector{Tv_out}(undef, K)
     if _is_periodic_bc(bc)
-        searcher = _resolve_search(x, xq, search, hint, bc)
+        # Helper wraps `x` via `_resolve_axis(x, bc)` and searches against the
+        # wrapped axis — axis dispatch handles seam, no `bc` thread needed.
+        searcher = _resolve_search(x, xq, search, hint)
         return _constant_oneshot_series_periodic!(output, x, s, xq, bc, deriv, side, searcher)
     end
     _check_domain(x, xq, extrap)
@@ -115,7 +117,7 @@ end
     length(output) == n_series(s) || _throw_series_dim_mismatch(length(output), n_series(s))
     x = _to_float(x, Tg)
     if _is_periodic_bc(bc)
-        searcher = _resolve_search(x, xq, search, hint, bc)
+        searcher = _resolve_search(x, xq, search, hint)
         return _constant_oneshot_series_periodic!(output, x, s, xq, bc, deriv, side, searcher)
     end
     _check_domain(x, xq, extrap)
@@ -154,7 +156,7 @@ end
     if _is_periodic_bc(bc)
         x_eff = _resolve_axis(x, bc)
         extrap_p = _resolve_extrap(NoExtrap(), bc, x_eff)
-        searcher = _resolve_search(x_eff, xqs, search, nothing, NoBC())
+        searcher = _resolve_search(x_eff, xqs, search, nothing)
         x_last = @inbounds Tg_actual(last(x_eff))
         @inbounds for j in 1:NQ
             xq_wrapped = _wrap_to_domain(xqs[j], x_eff)
@@ -205,7 +207,7 @@ end
     if _is_periodic_bc(bc)
         x_eff = _resolve_axis(x, bc)
         extrap_p = _resolve_extrap(NoExtrap(), bc, x_eff)
-        searcher = _resolve_search(x_eff, xqs, search, nothing, NoBC())
+        searcher = _resolve_search(x_eff, xqs, search, nothing)
         x_last = @inbounds Tg_actual(last(x_eff))
         aq_vec = acquire!(pool, _ConstantAnchoredQuery{Tg_actual, Tqp}, NQ)
         @inbounds for j in 1:NQ

diff --git a/src/constant/nd/constant_nd_oneshot.jl b/src/constant/nd/constant_nd_oneshot.jl
@@ -33,15 +33,14 @@ function _constant_interp_nd_oneshot(
     # passthrough or cached float form. Wrapper search returns post-fold
     # `idx_R = 1` at seam so `data[..., idx_R, ...]` indexes raw data.
     grids_eff = map(_resolve_axis, grids, bcs)
-    nobcs = ntuple(_ -> NoBC(), Val(N))
     # NoExtrap domain check must precede FillExtrap short-circuit
     _validate_nd_domain(grids_eff, query, extraps_val)
     oob_result = _try_fill_oob(query, grids_eff, extraps_val, EvalValue(), @inbounds first(data))
     oob_result !== nothing && return oob_result
 
     extraps_eff = _resolve_extrap(extraps_val, bcs, grids_eff, data, Val(N))
     q_eval = _handle_all_extraps(query, grids_eff, extraps_eff)
-    stencils, Ls, Rs = _search_all_intervals_stencil(q_eval, grids_eff, searches, hints, nobcs)
+    stencils, Ls, Rs = _search_all_intervals_stencil(q_eval, grids_eff, searches, hints)
     # 4-arg `_get_h(g, idx, xL, xR)` — cached path for `_CachedVector` (idx)
     # / `_CachedRange` (scalar field); raw `Vector` falls back to `xR - xL`.
     idxLs = map(first, stencils)
@@ -71,7 +70,6 @@ function _constant_interp_nd_oneshot_batch!(
     length(output) == nq || _throw_query_output_mismatch(nq, length(output))
     _query_validate(queries)
     grids_eff = map(_resolve_axis, grids, bcs)
-    nobcs = ntuple(_ -> NoBC(), Val(N))
     extraps_eff = _resolve_extrap(extraps_val, bcs, grids_eff, data, Val(N))
     # Batch-level InBounds promotion: see cubic_nd_oneshot.jl / linear_nd_oneshot.jl
     # for the same pattern. Replaces `_validate_nd_domain` (throw via 1D
@@ -84,7 +82,7 @@ function _constant_interp_nd_oneshot_batch!(
             output[k] = oob_val; continue
         end
         q_eval = _handle_all_extraps(query_k, grids_eff, extraps_eff)
-        stencils, Ls, Rs = _search_all_intervals_stencil(q_eval, grids_eff, policies, hints, nobcs)
+        stencils, Ls, Rs = _search_all_intervals_stencil(q_eval, grids_eff, policies, hints)
         idxLs = map(first, stencils)
         hs = map(_get_h, grids_eff, idxLs, Ls, Rs)
         output[k] = _constant_nd_kernel(data, stencils, hs, side_vals, q_eval, Ls)

diff --git a/src/core/nd_utils.jl b/src/core/nd_utils.jl
@@ -149,21 +149,21 @@ end
     error("unreachable: _first_fill_value called without FillExtrap")
 end
 
-# ── _resolve_extrap: ND variants (expand + promote [+ materialize]) ──
+# ── _resolve_extrap: ND variants (expand + promote [+ per-axis resolve]) ──
 #
 # Continues the `_resolve_extrap` family from `src/core/periodic.jl` (primitive
 # per-axis + 1D bundled + ND bundled-with-data). These ND methods handle the
 # scalar→NTuple expansion, periodic-BC override, and FillExtrap value-type
 # promotion. Two shapes by arity:
 #
-# - 4-arg (extrap, bcs, Val(N), Tv): expand + promote. Returns NTuple with
-#   possibly-unmaterialized `WrapExtrap{Nothing}` on periodic axes. Used by
-#   callers that materialize separately (post-extension persistent paths
-#   where bc-aware materialize would trip the pre-extension `<` check).
+# - 4-arg (extrap, bcs, Val(N), Tv): expand + promote. Returns NTuple per axis,
+#   forcing `WrapExtrap` where `bcs[d]` is periodic. Used by callers that
+#   resolve against the grid separately (post-extension persistent paths
+#   where the BC-aware check would trip the pre-extension `<` invariant).
 #
-# - 5-arg (extrap, bcs, grids, Val(N), Tv): above + per-axis materialize.
-#   `bcs::NTuple` → 3-arg primitive (bc-aware); `bcs::Nothing` → 2-arg primitive
-#   (grid-span only, no periodic override needed).
+# - 5-arg (extrap, bcs, grids, Val(N), Tv): above + per-axis passthrough
+#   against `grids`. `bcs::NTuple` → 3-arg primitive (bc-aware);
+#   `bcs::Nothing` → 2-arg primitive (no periodic override).
 
 # ── 4-arg: expand + promote (no materialize) ──
 
@@ -824,35 +824,31 @@ end
 # by Julia escape analysis → 0 heap bytes/call.
 
 # Per-axis inline: build Searcher + run search_interval in one body.
-# 3-arg `search_interval` (no spacing) — Range uses _search_direct's own step,
-# Vector uses _search_binary; the
-# stencil-using callers compute `h` from the search-returned `(xL, xR)` via
-# 3-arg `_get_h(grid, xL, xR)` dispatch.
-@inline _search_axis_stencil(grid, q, search, hint, bc) =
-    @inbounds search_interval(_resolve_search(grid, q, search, hint, bc), grid, q)
+# Callers pre-wrap grids via `_resolve_axis` (or pass already-wrapped axes),
+# so seam handling is via axis-level dispatch in `periodic_axis.jl`.
+@inline _search_axis_stencil(grid, q, search, hint) =
+    @inbounds search_interval(_resolve_search(grid, q, search, hint), grid, q)
 
 @inline function _search_all_intervals_stencil(
         q_evals::Tuple{Vararg{Real, N}},
         grids::Tuple{Vararg{AbstractVector, N}},
         searches::Tuple{Vararg{AbstractSearchPolicy, N}},
         hints::Tuple{Vararg{Base.RefValue{Int}, N}},
-        bcs::Tuple{Vararg{AbstractBC, N}},
     ) where {N}
-    results = map(_search_axis_stencil, grids, q_evals, searches, hints, bcs)
+    results = map(_search_axis_stencil, grids, q_evals, searches, hints)
     return _project_search_results(results, _getstencil)
 end
 
 # Nothing-hint overload — scalar oneshot entries only. Batch must use the
-# 5-arg NTuple form (hint allocation hoisted via `_resolve_oneshot_search_nd`).
+# 4-arg NTuple form (hint allocation hoisted via `_resolve_oneshot_search_nd`).
 @inline function _search_all_intervals_stencil(
         q_evals::Tuple{Vararg{Real, N}},
         grids::Tuple{Vararg{AbstractVector, N}},
         searches::Tuple{Vararg{AbstractSearchPolicy, N}},
         ::Nothing,
-        bcs::Tuple{Vararg{AbstractBC, N}},
     ) where {N}
     hints = _ensure_hint_nd(nothing, Val(N))
-    return _search_all_intervals_stencil(q_evals, grids, searches, hints, bcs)
+    return _search_all_intervals_stencil(q_evals, grids, searches, hints)
 end
 
 

diff --git a/src/core/periodic.jl b/src/core/periodic.jl
@@ -12,49 +12,18 @@
 # ========================================
 
 """
-    _wrap_to_domain(xi::FT, x_min::FT, x_max::FT) where {FT<:AbstractFloat}
+    _wrap_to_domain(xi, x_min, x_max)
 
-Wrap a query point `xi` to the domain [x_min, x_max].
-Used for periodic boundary conditions and extrap=WrapExtrap().
-
-Closed-domain convention: `xi == x_max` is an in-domain boundary query
-(returns `xi` unchanged); only strictly-OOB queries (`xi < x_min` or
-`xi > x_max`) take the cold `mod()` path. `PeriodicBC{:inclusive}` is
-forward-**value**-invariant because `y[1] ≈ y[end]` by construction; the
-adjoint sensitivity at the seam now scatters to slot `n` instead of slot `1`
-(delta-equivalent under the cycle constraint). `:exclusive` is fully invariant
-(forward + adjoint) because `_ExclusivePeriodicAxis.search_interval` already
-returns `idx_R = 1` for `xq >= inner[n]` at the seam.
-
-Optimized: skips expensive `mod()` when xi is already in domain.
+Wrap `xi` into `[x_min, x_max]` for `WrapExtrap` / periodic BC. Closed-domain:
+`xi == x_max` returns unchanged; only strictly-OOB queries take the `mod()`
+path. `_extract_primal` is identity on plain numerics, so the in-domain branch
+returns the original `xi` and preserves AD `Dual` carriers.
 """
-@inline function _wrap_to_domain(xi::Tg, x_min::Tg, x_max::Tg) where {Tg}
-    # Hot path: already in domain — return as-is (no arith). Cold `mod()`
-    # work goes through `@noinline` `_wrap_to_domain_slow` so it doesn't
-    # bloat the caller (every WrapExtrap eval kernel) with mod-related
-    # asm. On constant rng+perEx persistent (3-4 ns baseline, 138 lines
-    # before split), this collapses the eval kernel to ~75 lines.
-    if (xi >= x_min) && (xi <= x_max)
-        return xi
-    end
-    return _wrap_to_domain_slow(xi, x_min, x_max)
-end
-
-# Generic wrapper: handles Dual, Int, Float32 on Float64 grid, etc.
-# IMPORTANT: Preserves AD Dual type through the entire operation.
-# Same hot/cold split as the AbstractFloat overload above.
 @inline function _wrap_to_domain(xi::Real, x_min::Tg, x_max::Tg) where {Tg}
     xi_primal = _extract_primal(xi)
-    # Fast path: already in domain, return original xi (preserves Dual type for AD)
     if (xi_primal >= x_min) && (xi_primal <= x_max)
         return xi
     end
-    return _wrap_to_domain_slow(xi, x_min, x_max)
-end
-
-# Cold path — `mod()` work hoisted out of the inlined hot path.
-# `mod()` works correctly with `ForwardDiff.Dual`: d/dx[mod(x,p)] = 1.
-@noinline function _wrap_to_domain_slow(xi, x_min, x_max)
     period = x_max - x_min
     return x_min + mod(xi - x_min, period)
 end
@@ -456,43 +425,31 @@ end
 # Value extension: append first element
 _extend_values(y::AbstractVector) = vcat(y, first(y))
 
-# ========================================
-# WrapExtrap is a tag struct (eval_ops.jl)
-# ========================================
-#
-# After the surface-API axis resolution (`_resolve_axis` / `_cache_axis`),
-# every supported axis exposes `first/last` matching the canonical wrap
-# domain — including `_ExclusivePeriodicAxis`, whose `last` reports the
-# precomputed virtual endpoint. Eval kernels read those bounds directly
-# from the axis via `_wrap_to_domain(xq, x, ::WrapExtrap)`. No BC-aware
-# `WrapExtrap` constructors, no `WrapExtrap{Nothing}` materialization, no
-# duplicated `_x_min/_x_max` fields.
-
 # ========================================
 # Extrap Resolution (_resolve_extrap)
 # ========================================
 #
-# Single name for every extrap materialization / validation step. `extrap` is
-# always the 1st arg — consistent with `_resolve_search`, `_resolve_coeffs`,
-# `_resolve_grididx` naming family. Layers (dispatched by arg shape):
-#
-# 1D primitives (per-axis):
-#   (extrap, x)                          — grid-only; upgrade {Nothing} or passthrough
-#   (extrap, bc, x)                      — BC-aware; PeriodicBC forces WrapExtrap
+# `WrapExtrap` is a tag struct — the axis is the source of truth for the
+# wrap domain. After surface-API axis resolution (`_resolve_axis` /
+# `_cache_axis`), every supported axis exposes `first/last` matching the
+# canonical wrap domain — including `_ExclusivePeriodicAxis`, whose `last`
+# reports the precomputed virtual endpoint. Eval kernels read those bounds
+# directly via `_wrap_to_domain(xq, x)`.
 #
-# `WrapExtrap` is a tag struct — no `{Nothing}` placeholder, no materialization.
-# The axis carries the canonical wrap domain via `first/last` after the
-# surface-API axis resolution.
+# `_resolve_extrap` is the single name for every extrap normalization step
+# (BC override + FillExtrap value-type promotion). `extrap` is always the
+# 1st arg — consistent with `_resolve_search`, `_resolve_coeffs`,
+# `_resolve_grididx`. Layers (dispatched by arg shape):
 #
-# 1D entries (per-axis):
-#   (extrap, x)                          — passthrough + FillExtrap promote (no Tv → identity)
-#   (extrap, bc, x)                      — BC-aware: PeriodicBC forces WrapExtrap, otherwise passthrough
+# 1D primitives (per-axis):
+#   (extrap, x)                          — passthrough (tag-struct identity)
+#   (extrap, bc, x)                      — BC-aware: PeriodicBC forces WrapExtrap
 #   (extrap, x, Tv)                      — FillExtrap promote (eltype → Tv)
 #
 # 1D bundled: validate + dispatch
 #   (extrap, bc, x, y)                   — `:inclusive` endpoint check + primitive
 #
-# ND bundled (oneshot): slice validation + per-axis materialize
+# ND bundled (oneshot): slice validation + per-axis
 #   (extraps, bcs, grids, data, Val(N))  — zero-copy ND oneshot entry
 
 # ── Primitive: 2-arg (no BC info; non-periodic persistent / Hermite family) ──