From abcdd080e9c8c5f253c831c598c28794c5116fa5 Mon Sep 17 00:00:00 2001
From: Tharindu Jayasinghe <tharindujay96@gmail.com>
Date: Wed, 1 Jul 2026 14:03:27 -0700
Subject: [PATCH 1/9] perf: fast scatter shims + PDM fine-bin regrouping
 (CPU/torch/CUDA)

- scatter_add: numpy via bincount (buffered; ~10-30x over unbuffered
  add.at), cupy via cupyx.scatter_add
- new row-wise scatter_add_rows / scatter_counts_rows: torch scatters a
  stride-0 expanded view (no (P,N) broadcast copy, no flat index); numpy
  gets one fused bincount pass
- BLS/CE/TLS batch kernels use the row-wise shims; TLS also pads once to
  the widest template and accumulates the correlation in place; BLS
  hoists the per-duration column arange
- PDM bins each point once into n_bins*n_covers fine bins and regroups
  covers exactly (roll + group-sum): 3 scatters instead of 3*n_covers on
  the vectorized path, 3 shared-memory atomics per point instead of
  3*n_covers in the CUDA kernel. PDM numpy: 3652 -> 1289 ms on a
  3k-point/20k-period sweep, bit-identical theta.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
---
 src/cuperiod/core/_arrayapi.py              |  71 ++++++++++++-
 src/cuperiod/methods/_bls_core.py           |  19 ++--
 src/cuperiod/methods/conditional_entropy.py |   8 +-
 src/cuperiod/methods/pdm.py                 | 111 ++++++++++++--------
 src/cuperiod/methods/tls.py                 |  32 +++---
 5 files changed, 158 insertions(+), 83 deletions(-)

diff --git a/src/cuperiod/core/_arrayapi.py b/src/cuperiod/core/_arrayapi.py
index 48e6e0c..24593b2 100644
--- a/src/cuperiod/core/_arrayapi.py
+++ b/src/cuperiod/core/_arrayapi.py
@@ -93,7 +93,10 @@ def scatter_add(target: Any, index: Any, values: Any) -> None:
 
     ``target`` is 1-D and ``index``/``values`` are 1-D and aligned. ``numpy.add.at`` is
     not part of the array-API standard and PyTorch has no equivalent *function*, so this
-    dispatches per backend: torch uses ``Tensor.index_add_``; numpy/cupy use ``add.at``.
+    dispatches per backend: torch uses ``Tensor.index_add_``; cupy uses
+    ``cupyx.scatter_add`` (its fast documented scatter); numpy uses ``bincount``, which
+    computes the same float64 sums in one buffered pass — ``np.add.at`` is unbuffered
+    and an order of magnitude slower on large inputs.
 
     ``values.dtype`` must equal ``target.dtype``: torch ``index_add_`` rejects a
     mismatch (numpy/cupy would silently cast), so callers build both at the same
@@ -102,12 +105,72 @@ def scatter_add(target: Any, index: Any, values: Any) -> None:
     if is_torch_array(target):
         target.index_add_(0, index.long(), values)
         return
+    if is_cupy_array(target):
+        import cupyx
+
+        cupyx.scatter_add(target, index, values)
+        return
+    binned = np.bincount(
+        np.asarray(index), weights=np.asarray(values), minlength=target.size
+    )
+    target += binned.astype(target.dtype, copy=False)
+
+
+def scatter_add_rows(target: Any, index: Any, values: Any) -> None:
+    """In-place ``target[p, index[p, j]] += values[j]`` with repeats accumulated.
+
+    The batch-kernel binning primitive: ``target`` is 2-D ``(P, W)``, ``index`` is
+    ``(P, N)`` int64, and ``values`` is one shared row ``(N,)`` (the usual case — the
+    same light curve binned at ``P`` trial periods) or a full ``(P, N)``. Compared to
+    flattening and calling :func:`scatter_add`, this avoids materializing the broadcast
+    ``(P, N)`` values *and* the flat ``(P, N)`` int64 index on torch — the scatter
+    reads a stride-0 expanded view — and gives numpy one fused ``bincount`` pass.
+    """
+    if is_torch_array(target):
+        src = values if values.ndim == 2 else values.unsqueeze(0).expand_as(index)
+        target.scatter_add_(1, index, src)
+        return
+    if is_cupy_array(target):
+        import cupy
+        import cupyx
+
+        rows = cupy.arange(target.shape[0])[:, None]
+        cupyx.scatter_add(target, (rows, index), values)
+        return
+    n_rows, n_cols = target.shape
+    flat = index + (np.arange(n_rows, dtype=np.int64) * n_cols)[:, None]
+    weights = np.broadcast_to(values, index.shape)
+    binned = np.bincount(
+        flat.reshape(-1), weights=weights.reshape(-1), minlength=target.size
+    )
+    target += binned.reshape(target.shape).astype(target.dtype, copy=False)
+
+
+def scatter_counts_rows(target: Any, index: Any) -> None:
+    """In-place ``target[p, index[p, j]] += 1`` — a per-row histogram count.
+
+    Same layout contract as :func:`scatter_add_rows` but without a values array:
+    numpy uses weightless ``bincount`` and torch scatters a stride-0 view of a single
+    one, so no ``(P, N)`` ones array is ever built. ``target`` is float (counts are
+    accumulated in the working dtype for the entropy/variance math downstream).
+    """
+    if is_torch_array(target):
+        import torch
+
+        one = torch.ones((1, 1), dtype=target.dtype, device=target.device)
+        target.scatter_add_(1, index, one.expand_as(index))
+        return
     if is_cupy_array(target):
         import cupy
+        import cupyx
 
-        cupy.add.at(target, index, values)
+        rows = cupy.arange(target.shape[0])[:, None]
+        cupyx.scatter_add(target, (rows, index), target.dtype.type(1))
         return
-    np.add.at(target, index, values)
+    n_rows, n_cols = target.shape
+    flat = index + (np.arange(n_rows, dtype=np.int64) * n_cols)[:, None]
+    binned = np.bincount(flat.reshape(-1), minlength=target.size)
+    target += binned.reshape(target.shape).astype(target.dtype, copy=False)
 
 
 def to_host(a: Any) -> np.ndarray:
@@ -201,6 +264,8 @@ def to_device_array(host: np.ndarray, *, device: str, dtype: Any) -> Any:
     "resolve_precision",
     "resolve_torch_device",
     "scatter_add",
+    "scatter_add_rows",
+    "scatter_counts_rows",
     "to_device_array",
     "to_host",
 ]
diff --git a/src/cuperiod/methods/_bls_core.py b/src/cuperiod/methods/_bls_core.py
index cf4cebb..5594129 100644
--- a/src/cuperiod/methods/_bls_core.py
+++ b/src/cuperiod/methods/_bls_core.py
@@ -28,7 +28,7 @@
     array_namespace,
     device_ref,
     resolve_precision,
-    scatter_add,
+    scatter_add_rows,
     to_device_array,
     to_host,
 )
@@ -133,7 +133,6 @@ def _bls_search(
     fdtype = periods.dtype
     idtype = xp.int64
     dev = device_ref(periods)
-    n_points = int(t.shape[0])
     yw = y * ivar
     sum_y = float(xp.sum(yw))
     sum_ivar = float(xp.sum(ivar))
@@ -152,6 +151,7 @@ def _bls_search(
     }
     if n_periods == 0 or not dur_bins:
         return out
+    cols_full = xp.arange(width, dtype=idtype, device=dev)
 
     for start in range(0, n_periods, batch):
         stop = min(start + batch, n_periods)
@@ -163,15 +163,10 @@ def _bls_search(
         phase_t = xp.remainder(tau[None, :], pb[:, None])
         ind = xp.astype(phase_t / bin_duration, idtype) + 1
         ind = xp.clip(ind, 0, width - 1)
-        flat = xp.reshape(rows[:, None] * width + ind, (-1,))
-        mean_y = xp.zeros(n_p * width, dtype=fdtype, device=dev)
-        mean_ivar = xp.zeros(n_p * width, dtype=fdtype, device=dev)
-        yw_b = xp.reshape(xp.broadcast_to(yw, (n_p, n_points)), (-1,))
-        ivar_b = xp.reshape(xp.broadcast_to(ivar, (n_p, n_points)), (-1,))
-        scatter_add(mean_y, flat, yw_b)
-        scatter_add(mean_ivar, flat, ivar_b)
-        mean_y = xp.reshape(mean_y, (n_p, width))
-        mean_ivar = xp.reshape(mean_ivar, (n_p, width))
+        mean_y = xp.zeros((n_p, width), dtype=fdtype, device=dev)
+        mean_ivar = xp.zeros((n_p, width), dtype=fdtype, device=dev)
+        scatter_add_rows(mean_y, ind, yw)
+        scatter_add_rows(mean_ivar, ind, ivar)
 
         for j in range(oversample):
             dst = xp.clip(n_bins - oversample + j, 0, width - 1)
@@ -191,7 +186,7 @@ def _bls_search(
             y_in = cy[:, kd:] - cy[:, :-kd]
             ivar_in = cw[:, kd:] - cw[:, :-kd]
             ivar_out = sum_ivar - ivar_in
-            cols = xp.arange(y_in.shape[1], dtype=idtype, device=dev)
+            cols = cols_full[: width - kd]
             valid = (
                 (cols[None, :] <= (n_bins[:, None] - kd))
                 & (ivar_in >= _IVAR_EPS)
diff --git a/src/cuperiod/methods/conditional_entropy.py b/src/cuperiod/methods/conditional_entropy.py
index e3ec696..4408751 100644
--- a/src/cuperiod/methods/conditional_entropy.py
+++ b/src/cuperiod/methods/conditional_entropy.py
@@ -22,7 +22,7 @@
     device_ref,
     resolve_precision,
     resolve_torch_device,
-    scatter_add,
+    scatter_counts_rows,
     to_device_array,
     to_host,
 )
@@ -72,13 +72,11 @@ def _entropy_batch(
         stop = min(start + batch, n_periods)
         pb = periods[start:stop]
         n_p = int(pb.shape[0])
-        rows = xp.arange(n_p, dtype=idtype, device=dev)
         phase = xp.remainder(tau[None, :] / pb[:, None], 1.0)
         phase_bin = xp.clip(xp.astype(phase * n_phase, idtype), 0, n_phase - 1)
         cell = phase_bin * n_mag + mag_bin[None, :]  # (P, N) in [0, n_cells)
-        flat = xp.reshape(rows[:, None] * n_cells + cell, (-1,))
-        count = xp.zeros(n_p * n_cells, dtype=fdtype, device=dev)
-        scatter_add(count, flat, xp.ones(n_p * n_points, dtype=fdtype, device=dev))
+        count = xp.zeros((n_p, n_cells), dtype=fdtype, device=dev)
+        scatter_counts_rows(count, cell)
         count = xp.reshape(count, (n_p, n_phase, n_mag))
 
         phase_total = xp.sum(count, axis=2, keepdims=True)  # (P, n_phase, 1)
diff --git a/src/cuperiod/methods/pdm.py b/src/cuperiod/methods/pdm.py
index 6237d22..ebbc1e8 100644
--- a/src/cuperiod/methods/pdm.py
+++ b/src/cuperiod/methods/pdm.py
@@ -27,7 +27,8 @@
     device_ref,
     resolve_precision,
     resolve_torch_device,
-    scatter_add,
+    scatter_add_rows,
+    scatter_counts_rows,
     to_device_array,
     to_host,
 )
@@ -50,6 +51,24 @@
 DEFAULT_BATCH: Final = 2048
 
 
+def _cover_hists(xp: ModuleType, fine: Any, n_bins: int, n_covers: int) -> Any:
+    """All covers' bin histograms from one fine histogram, laid side by side.
+
+    ``fine`` is ``(P, n_bins*n_covers)``; cover ``c``'s ``n_bins`` histogram is an
+    exact integer regroup of the fine bins — roll right by ``c`` and sum groups of
+    ``n_covers`` — so the folded points are binned **once** rather than per cover.
+    Returns ``(P, n_bins*n_covers)`` with cover ``c`` at columns ``[c*n_bins, ...)``.
+    """
+    n_p = int(fine.shape[0])
+    parts = [
+        xp.sum(
+            xp.reshape(xp.roll(fine, cover, axis=1), (n_p, n_bins, n_covers)), axis=2
+        )
+        for cover in range(n_covers)
+    ]
+    return xp.concat(parts, axis=1)
+
+
 def _theta_batch(
     xp: ModuleType,
     tau: Any,
@@ -65,7 +84,10 @@ def _theta_batch(
     """Stellingwerf Theta for each trial period, vectorized over ``periods``.
 
     ``s^2 = sum_j SSD_j / (N*n_covers - n_nonempty)`` with ``SSD_j`` the within-bin sum
-    of squared deviations across all covers, and ``Theta = s^2 / sigma^2``.
+    of squared deviations across all covers, and ``Theta = s^2 / sigma^2``. The points
+    are binned once into ``n_bins*n_covers`` fine bins; every cover's bin statistics
+    are exact regroups of that histogram (see :func:`_cover_hists`), which cuts the
+    scatter work by ``n_covers``.
     """
     fdtype = periods.dtype
     idtype = xp.int64
@@ -74,34 +96,26 @@ def _theta_batch(
     n_points = int(tau.shape[0])
     n_periods = int(periods.shape[0])
     n_global = n_bins * n_covers
-    cover_step = 1.0 / (n_bins * n_covers)
     theta = xp.empty(n_periods, dtype=fdtype, device=dev)
 
     for start in range(0, n_periods, batch):
         stop = min(start + batch, n_periods)
         pb = periods[start:stop]
         n_p = int(pb.shape[0])
-        rows = xp.arange(n_p, dtype=idtype, device=dev)
         phase = xp.remainder(tau[None, :] / pb[:, None], 1.0)  # (P, N) in [0, 1)
-
-        count = xp.zeros(n_p * n_global, dtype=fdtype, device=dev)
-        ysum = xp.zeros(n_p * n_global, dtype=fdtype, device=dev)
-        ysq = xp.zeros(n_p * n_global, dtype=fdtype, device=dev)
-        ones = xp.ones(n_p * n_points, dtype=fdtype, device=dev)
-        y_b = xp.reshape(xp.broadcast_to(y, (n_p, n_points)), (-1,))
-        y2_b = xp.reshape(xp.broadcast_to(y2, (n_p, n_points)), (-1,))
-        for cover in range(n_covers):
-            offset = cover * cover_step
-            b = xp.astype(xp.remainder(phase + offset, 1.0) * n_bins, idtype)
-            b = xp.clip(b, 0, n_bins - 1)
-            flat = xp.reshape(rows[:, None] * n_global + (b + cover * n_bins), (-1,))
-            scatter_add(count, flat, ones)
-            scatter_add(ysum, flat, y_b)
-            scatter_add(ysq, flat, y2_b)
-
-        count = xp.reshape(count, (n_p, n_global))
-        ysum = xp.reshape(ysum, (n_p, n_global))
-        ysq = xp.reshape(ysq, (n_p, n_global))
+        fine = xp.clip(
+            xp.astype(phase * float(n_global), idtype), 0, n_global - 1
+        )
+        f_count = xp.zeros((n_p, n_global), dtype=fdtype, device=dev)
+        f_ysum = xp.zeros((n_p, n_global), dtype=fdtype, device=dev)
+        f_ysq = xp.zeros((n_p, n_global), dtype=fdtype, device=dev)
+        scatter_counts_rows(f_count, fine)
+        scatter_add_rows(f_ysum, fine, y)
+        scatter_add_rows(f_ysq, fine, y2)
+
+        count = _cover_hists(xp, f_count, n_bins, n_covers)
+        ysum = _cover_hists(xp, f_ysum, n_bins, n_covers)
+        ysq = _cover_hists(xp, f_ysq, n_bins, n_covers)
         mask = count > 0.0
         safe = xp.where(mask, count, 1.0)
         ssd = xp.where(mask, ysq - ysum * ysum / safe, 0.0)
@@ -116,8 +130,11 @@ def _theta_batch(
 #: CUDA threads per block (one block per trial period).
 CUDA_BLOCK: Final = 128
 
-#: One-block-per-period PDM kernel: a block bins its folded points in shared memory
-#: (sum, sum of squares, count per cover) and reduces them to Theta on one thread.
+#: One-block-per-period PDM kernel: a block bins its folded points **once** into
+#: ``n_bins*n_covers`` fine shared-memory bins (sum, sum of squares, count); every
+#: cover's bin statistics are exact regroups of ``n_covers`` consecutive fine bins
+#: (rolled by the cover index), assembled in the one-thread reduction. Cuts the
+#: dominant atomic work per point from ``3*n_covers`` to 3.
 _PDM_CUDA_SRC: Final = r"""
 extern "C" __global__ void pdm_block(
     const double* __restrict__ tau, const double* __restrict__ y,
@@ -133,41 +150,45 @@ def _theta_batch(
     const int M = n_bins * n_covers;
 
     extern __shared__ double sh[];
-    double* s_sum = sh;          // (M) sum of y per bin
-    double* s_sq = sh + M;       // (M) sum of y^2 per bin
-    double* s_cnt = sh + 2 * M;  // (M) point count per bin
+    double* s_sum = sh;          // (M) sum of y per fine bin
+    double* s_sq = sh + M;       // (M) sum of y^2 per fine bin
+    double* s_cnt = sh + 2 * M;  // (M) point count per fine bin
     for (int i = tid; i < M; i += nth) {
         s_sum[i] = 0.0; s_sq[i] = 0.0; s_cnt[i] = 0.0;
     }
     __syncthreads();
 
-    const double cover_step = 1.0 / ((double)n_bins * (double)n_covers);
     for (int j = tid; j < n_points; j += nth) {
         double x = tau[j];
         double ph = (x - period * floor(x / period)) / period;  // mod(tau/period, 1)
         double yj = y[j];
-        for (int c = 0; c < n_covers; ++c) {
-            double pc = ph + (double)c * cover_step;
-            pc -= floor(pc);
-            int b = (int)(pc * n_bins);
-            if (b >= n_bins) b = n_bins - 1;
-            if (b < 0) b = 0;
-            int gid = c * n_bins + b;
-            atomicAdd(&s_sum[gid], yj);
-            atomicAdd(&s_sq[gid], yj * yj);
-            atomicAdd(&s_cnt[gid], 1.0);
-        }
+        int f = (int)(ph * (double)M);
+        if (f >= M) f = M - 1;
+        if (f < 0) f = 0;
+        atomicAdd(&s_sum[f], yj);
+        atomicAdd(&s_sq[f], yj * yj);
+        atomicAdd(&s_cnt[f], 1.0);
     }
     __syncthreads();
 
     if (tid == 0) {
         double ssd = 0.0;
         int nonempty = 0;
-        for (int i = 0; i < M; ++i) {
-            double cnt = s_cnt[i];
-            if (cnt > 0.0) {
-                ssd += s_sq[i] - s_sum[i] * s_sum[i] / cnt;
-                nonempty += 1;
+        for (int c = 0; c < n_covers; ++c) {
+            for (int b = 0; b < n_bins; ++b) {
+                // cover-c bin b = fine bins (b*n_covers - c .. +n_covers-1) mod M
+                double cnt = 0.0, sum = 0.0, sq = 0.0;
+                int f0 = b * n_covers - c;
+                if (f0 < 0) f0 += M;
+                for (int k = 0; k < n_covers; ++k) {
+                    int f = f0 + k;
+                    if (f >= M) f -= M;
+                    cnt += s_cnt[f]; sum += s_sum[f]; sq += s_sq[f];
+                }
+                if (cnt > 0.0) {
+                    ssd += sq - sum * sum / cnt;
+                    nonempty += 1;
+                }
             }
         }
         double den = (double)(n_points * n_covers - nonempty);
diff --git a/src/cuperiod/methods/tls.py b/src/cuperiod/methods/tls.py
index 1cd9f4e..9135baf 100644
--- a/src/cuperiod/methods/tls.py
+++ b/src/cuperiod/methods/tls.py
@@ -29,7 +29,7 @@
     device_ref,
     resolve_precision,
     resolve_torch_device,
-    scatter_add,
+    scatter_add_rows,
     to_device_array,
     to_host,
 )
@@ -137,7 +137,6 @@ def _matched_filter(
         "duration": xp.zeros(n_periods, dtype=fdtype, device=dev),
         "t0": xp.zeros(n_periods, dtype=fdtype, device=dev),
     }
-    n_points = int(tau.shape[0])
     for start in range(0, n_periods, period_batch):
         stop = min(start + period_batch, n_periods)
         pb = periods[start:stop]
@@ -145,15 +144,14 @@ def _matched_filter(
         rows_p = xp.arange(n_p, dtype=idtype, device=dev)
         phase = xp.remainder(tau[None, :] / pb[:, None], 1.0)
         bin_idx = xp.clip(xp.astype(phase * n_bins, idtype), 0, n_bins - 1)
-        flat = xp.reshape(rows_p[:, None] * n_bins + bin_idx, (-1,))
-        a_flat = xp.zeros(n_p * n_bins, dtype=fdtype, device=dev)  # sum w*y' per bin
-        b_flat = xp.zeros(n_p * n_bins, dtype=fdtype, device=dev)  # sum w per bin
-        yw_b = xp.reshape(xp.broadcast_to(yw, (n_p, n_points)), (-1,))
-        w_b = xp.reshape(xp.broadcast_to(w, (n_p, n_points)), (-1,))
-        scatter_add(a_flat, flat, yw_b)
-        scatter_add(b_flat, flat, w_b)
-        a = xp.reshape(a_flat, (n_p, n_bins))
-        b = xp.reshape(b_flat, (n_p, n_bins))
+        a = xp.zeros((n_p, n_bins), dtype=fdtype, device=dev)  # sum w*y' per bin
+        b = xp.zeros((n_p, n_bins), dtype=fdtype, device=dev)  # sum w per bin
+        scatter_add_rows(a, bin_idx, yw)
+        scatter_add_rows(b, bin_idx, w)
+        # One circular pad sized for the widest template serves every width below.
+        max_pad = max(dur_bins) - 1
+        a_ext = xp.concat([a, a[:, :max_pad]], axis=1) if max_pad else a
+        b_ext = xp.concat([b, b[:, :max_pad]], axis=1) if max_pad else b
 
         best_sr = xp.zeros(n_p, dtype=fdtype, device=dev)
         best_depth = xp.zeros(n_p, dtype=fdtype, device=dev)
@@ -161,14 +159,12 @@ def _matched_filter(
         best_width = xp.zeros(n_p, dtype=idtype, device=dev)
         for width in dur_bins:
             g = templates[width]
-            a_ext = xp.concat([a, a[:, : width - 1]], axis=1)
-            b_ext = xp.concat([b, b[:, : width - 1]], axis=1)
-            num = xp.zeros((n_p, n_bins), dtype=fdtype, device=dev)
-            den = xp.zeros((n_p, n_bins), dtype=fdtype, device=dev)
-            for k in range(width):  # correlate the folded data with the template
+            num = float(g[0]) * a_ext[:, :n_bins]
+            den = float(g[0] * g[0]) * b_ext[:, :n_bins]
+            for k in range(1, width):  # correlate the folded data with the template
                 gk = float(g[k])
-                num = num + gk * a_ext[:, k : k + n_bins]
-                den = den + (gk * gk) * b_ext[:, k : k + n_bins]
+                num += gk * a_ext[:, k : k + n_bins]
+                den += (gk * gk) * b_ext[:, k : k + n_bins]
             safe_den = xp.where(den > _W_EPS, den, 1.0)
             # a dip means the in-transit weighted residual is negative -> num < 0
             sr = xp.where((den > _W_EPS) & (num < 0.0), num * num / safe_den, 0.0)

From a6bb69d4993798f835d1b63c9be76138ad521005 Mon Sep 17 00:00:00 2001
From: Tharindu Jayasinghe <tharindujay96@gmail.com>
Date: Wed, 1 Jul 2026 14:08:39 -0700
Subject: [PATCH 2/9] perf(mhaov): assemble the Gram from harmonic trig sums,
 not a design tensor

Every normal-equation entry is analytically a harmonic trig sum
(product-to-sum identities), so compute C_m/S_m for m<=2H via the
Chebyshev recurrence from ONE cos/sin evaluation and assemble the tiny
(F,d,d) Gram from them: O(F*N*d^2) -> O(F*N*H) work, O(d) less transient
memory (no more multi-GB design tensors at 1e5 points), and gemm-free on
every backend by construction, which retires the Blackwell cuBLAS
workaround branch. Same math, same ridge; checksums unchanged.

20k points x 30k frequencies, H=3: numpy 81.6s -> 30.1s, cupy 2.04s ->
0.71s, torch:cuda 2.12s -> 0.66s.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
---
 src/cuperiod/methods/mhaov.py | 109 +++++++++++++++++++++++++---------
 1 file changed, 82 insertions(+), 27 deletions(-)

diff --git a/src/cuperiod/methods/mhaov.py b/src/cuperiod/methods/mhaov.py
index 740af52..2d3474c 100644
--- a/src/cuperiod/methods/mhaov.py
+++ b/src/cuperiod/methods/mhaov.py
@@ -28,8 +28,6 @@
 from cuperiod.core._arrayapi import (
     array_namespace,
     device_ref,
-    is_cupy_array,
-    is_torch_array,
     resolve_precision,
     resolve_torch_device,
     to_device_array,
@@ -64,17 +62,44 @@
 _RIDGE_EPS: Final = 1.0e3
 
 
-def _design(xp: ModuleType, angle: Any, n_harmonics: int) -> Any:
-    """Trig-polynomial design tensor ``(F, N, 2H+1)`` = [1, cos kθ, sin kθ]."""
-    n_freq, n_points = angle.shape
-    d = 2 * n_harmonics + 1
-    dev = device_ref(angle)
-    design = xp.empty((n_freq, n_points, d), dtype=angle.dtype, device=dev)
-    design[:, :, 0] = 1.0
-    for k in range(1, n_harmonics + 1):
-        design[:, :, 2 * k - 1] = xp.cos(k * angle)
-        design[:, :, 2 * k] = xp.sin(k * angle)
-    return design
+def _harmonic_sums(
+    xp: ModuleType, angle: Any, y: Any, n_harmonics: int
+) -> tuple[list[Any], list[Any], list[Any], list[Any]]:
+    """Per-frequency harmonic trig sums from one ``(F, N)`` angle matrix.
+
+    Returns ``(C, S, A, B)`` with ``C[m] = Σ_n cos(m·θ)``, ``S[m] = Σ_n sin(m·θ)`` for
+    ``m = 1..2H`` (index 0 unused: those sums are the constants ``N`` and ``0``) and
+    ``A[k] = Σ_n y_n cos(k·θ)``, ``B[k] = Σ_n y_n sin(k·θ)`` for ``k = 1..H``. The
+    higher harmonics come from the Chebyshev recurrence ``T_m = 2·cosθ·T_{m-1} −
+    T_{m-2}``, so ``cos``/``sin`` are evaluated **once** regardless of ``H``.
+    """
+    m_max = 2 * n_harmonics
+    c1 = xp.cos(angle)
+    s1 = xp.sin(angle)
+    C: list[Any] = [None] * (m_max + 1)
+    S: list[Any] = [None] * (m_max + 1)
+    A: list[Any] = [None] * (n_harmonics + 1)
+    B: list[Any] = [None] * (n_harmonics + 1)
+    two_c1 = 2.0 * c1
+    cm, sm = c1, s1
+    prev_c: Any = None
+    prev_s: Any = None
+    for m in range(1, m_max + 1):
+        C[m] = xp.sum(cm, axis=1)
+        S[m] = xp.sum(sm, axis=1)
+        if m <= n_harmonics:
+            A[m] = xp.sum(y[None, :] * cm, axis=1)
+            B[m] = xp.sum(y[None, :] * sm, axis=1)
+        if m < m_max:
+            if m == 1:  # T_0 is the constant 1 / 0, kept out of the arrays
+                next_c = two_c1 * cm - 1.0
+                next_s = two_c1 * sm
+            else:
+                next_c = two_c1 * cm - prev_c
+                next_s = two_c1 * sm - prev_s
+            prev_c, prev_s = cm, sm
+            cm, sm = next_c, next_s
+    return C, S, A, B
 
 
 def _model_ss_batch(
@@ -93,33 +118,63 @@ def _model_ss_batch(
 
     This is the projection norm of the data onto the ``2H+1`` harmonic basis; the AOV
     F-statistic (single- or multi-band) is formed from it by the callers.
+
+    Every entry of the normal equations is analytically a harmonic trig sum — by the
+    product-to-sum identities, ``Σ cos(kθ)cos(lθ) = ½(C_{|k−l|} + C_{k+l})`` and so on
+    — so instead of materializing the ``(F, N, 2H+1)`` design tensor and forming its
+    ``O(F·N·d²)`` Gram, the sums ``C_m``/``S_m`` (``m ≤ 2H``) are computed in
+    ``O(F·N·H)`` from one ``cos``/``sin`` evaluation (:func:`_harmonic_sums`) and the
+    tiny ``(F, d, d)`` Gram is assembled from them. Mathematically identical (same
+    normal equations, same ridge), gemm-free on every backend, and ``O(d)`` less
+    transient memory.
     """
     d = 2 * n_harmonics + 1
     n_freq = int(frequencies.shape[0])
     fdtype = frequencies.dtype
     dev = device_ref(frequencies)
     ridge = _RIDGE_EPS * float(xp.finfo(fdtype).eps) * float(n_points)
-    eye = xp.eye(d, dtype=fdtype, device=dev) * ridge
     out = xp.empty(n_freq, dtype=fdtype, device=dev)
     two_pi = 2.0 * float(np.pi)
+    sum_y = float(xp.sum(y))
+    h = n_harmonics
 
     for start in range(0, n_freq, batch):
         stop = min(start + batch, n_freq)
         fb = frequencies[start:stop]
+        n_f = int(fb.shape[0])
         angle = (two_pi * fb)[:, None] * tau[None, :]
-        design = _design(xp, angle, n_harmonics)
-        if is_cupy_array(design) or is_torch_array(design):
-            # GPU batched cuBLAS gemm intermittently raises CUBLAS_STATUS_INVALID_VALUE
-            # on some GPUs (Blackwell / sm_120, in both cupy and torch). These gemm-free
-            # reductions match the einsum; the ``d`` loop avoids a (F,N,d,d) blowup.
-            gram = eye + xp.stack(
-                [xp.sum(design * design[:, :, j : j + 1], axis=1) for j in range(d)],
-                axis=-1,
-            )
-            proj = xp.sum(design * y[None, :, None], axis=1)
-        else:
-            gram = xp.einsum("fni,fnj->fij", design, design) + eye
-            proj = xp.einsum("fni,n->fi", design, y)
+        C, S, A, B = _harmonic_sums(xp, angle, y, h)
+        del angle
+
+        gram = xp.zeros((n_f, d, d), dtype=fdtype, device=dev)
+        proj = xp.empty((n_f, d), dtype=fdtype, device=dev)
+        gram[:, 0, 0] = float(n_points) + ridge
+        proj[:, 0] = sum_y
+        for k in range(1, h + 1):
+            gram[:, 0, 2 * k - 1] = C[k]
+            gram[:, 2 * k - 1, 0] = C[k]
+            gram[:, 0, 2 * k] = S[k]
+            gram[:, 2 * k, 0] = S[k]
+            proj[:, 2 * k - 1] = A[k]
+            proj[:, 2 * k] = B[k]
+            for line in range(1, h + 1):
+                # Σ cos·cos, Σ sin·sin, Σ cos·sin over n from the m-sums; C_0 = N,
+                # S_0 = 0, S_{-m} = -S_m.
+                m_diff = abs(k - line)
+                c_diff = C[m_diff] if m_diff else float(n_points)
+                cc = 0.5 * (c_diff + C[k + line])
+                ss = 0.5 * (c_diff - C[k + line])
+                if k == line:
+                    s_diff = 0.0
+                elif k > line:
+                    s_diff = S[m_diff]
+                else:
+                    s_diff = -S[m_diff]
+                cs = 0.5 * (S[k + line] - s_diff)
+                gram[:, 2 * k - 1, 2 * line - 1] = cc + (ridge if k == line else 0.0)
+                gram[:, 2 * k, 2 * line] = ss + (ridge if k == line else 0.0)
+                gram[:, 2 * k - 1, 2 * line] = cs
+                gram[:, 2 * line, 2 * k - 1] = cs
         # numpy 2.x batched solve treats a 2-D RHS as matrices, so add a trailing
         # singleton to keep it a per-frequency vector solve.
         beta = xp.linalg.solve(gram, proj[..., None])[..., 0]

From b7caafdfdecba3586e9fb92781602d5e188ffdd6 Mon Sep 17 00:00:00 2001
From: Tharindu Jayasinghe <tharindujay96@gmail.com>
Date: Wed, 1 Jul 2026 14:10:51 -0700
Subject: [PATCH 3/9] perf(gls): one trig evaluation for all six direct-path
 sums

The base-grid sums for w and w*y share one angle matrix, and the
doubled-frequency sums follow exactly from the double-angle identities
(cos2t = 2cos^2 t - 1, sin2t = 2 sin t cos t), so the portable direct
path now evaluates cos/sin once per frequency instead of six times.
The frequency chunk is also auto-capped by N so the (chunk, N)
transients stay ~64 MB on any device.

5k points x 100k frequencies: torch:cpu 3.06s -> 1.66s, torch:cuda
264ms -> 150ms; max deviation vs finufft unchanged at 2.2e-11.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
---
 src/cuperiod/methods/gls.py | 64 ++++++++++++++++++++++++-------------
 1 file changed, 41 insertions(+), 23 deletions(-)

diff --git a/src/cuperiod/methods/gls.py b/src/cuperiod/methods/gls.py
index a31a320..26d6e95 100644
--- a/src/cuperiod/methods/gls.py
+++ b/src/cuperiod/methods/gls.py
@@ -209,31 +209,51 @@ def lombscargle_power(
     )
 
 
-def _trig_sums_direct(
-    xp: Any, tau: Any, strengths: Any, f0: float, df: float, nf: int, *, freq_batch: int
-) -> tuple[Any, Any]:
-    """Direct (NUFFT-free) trig sums on the uniform grid ``f_k = f0 + k*df``.
-
-    Returns ``(cos_sum, sin_sum)`` real arrays of length ``nf`` with
-    ``cos_sum[k] = sum_j strengths_j cos(2*pi*f_k*tau_j)`` and ``sin_sum`` the matching
-    ``+sin`` sum (the ``isign=+1`` convention of the NUFFT path). This is ``O(N*nf)``
-    rather than the NUFFT's ``O(nf log nf)``, but is pure array-API and so runs on any
-    backend/device — the portable GLS path for AMD/Intel/Mac. Batched over frequency to
-    bound the transient ``(chunk, N)`` angle matrix.
+#: Elements per transient ``(chunk, N)`` matrix in the direct path (64 MB float64);
+#: the frequency chunk is capped so long light curves cannot blow device memory.
+_DIRECT_CHUNK_ELEMS: Final = 1 << 23
+
+
+def _gls_sums_direct(
+    xp: Any, tau: Any, w: Any, wy: Any, f0: float, df: float, nf: int, *,
+    freq_batch: int,
+) -> tuple[Any, Any, Any, Any, Any, Any]:
+    """All six GLS trig sums on ``f_k = f0 + k*df`` from one cos/sin evaluation.
+
+    Returns ``(c, s, yc, ys, c2, s2)`` — the cos/sin sums of the weights ``w`` and the
+    weighted data ``wy`` on the base grid, and of ``w`` on the doubled grid (the
+    ``isign=+1`` convention of the NUFFT path). The base-grid sums share one angle
+    matrix, and the doubled-frequency sums follow from the double-angle identities
+    ``cos2θ = 2cos²θ − 1`` and ``sin2θ = 2·sinθ·cosθ``, so ``cos``/``sin`` are
+    evaluated **once** per frequency instead of six times (twice per grid in the
+    three-call arrangement). ``O(N*nf)`` rather than the NUFFT's ``O(nf log nf)``, but
+    pure array-API — the portable GLS path for AMD/Intel/Mac. Batched over frequency
+    to bound the transient ``(chunk, N)`` matrices, with the chunk auto-capped so the
+    transients stay ~64 MB regardless of ``N``.
     """
     two_pi = 2.0 * float(np.pi)
     fdtype = tau.dtype
     dev = device_ref(tau)
+    n_points = int(tau.shape[0])
+    chunk = max(1, min(freq_batch, _DIRECT_CHUNK_ELEMS // max(1, n_points)))
     freqs = f0 + df * xp.arange(nf, dtype=fdtype, device=dev)
-    cos_sum = xp.empty(nf, dtype=fdtype, device=dev)
-    sin_sum = xp.empty(nf, dtype=fdtype, device=dev)
-    strength_row = strengths[None, :]
-    for start in range(0, nf, freq_batch):
-        stop = min(start + freq_batch, nf)
+    c, s, yc, ys, c2, s2 = (xp.empty(nf, dtype=fdtype, device=dev) for _ in range(6))
+    sum_w = float(xp.sum(w))
+    w_row = w[None, :]
+    wy_row = wy[None, :]
+    for start in range(0, nf, chunk):
+        stop = min(start + chunk, nf)
         ang = (two_pi * freqs[start:stop])[:, None] * tau[None, :]
-        cos_sum[start:stop] = xp.sum(strength_row * xp.cos(ang), axis=1)
-        sin_sum[start:stop] = xp.sum(strength_row * xp.sin(ang), axis=1)
-    return cos_sum, sin_sum
+        cos_a = xp.cos(ang)
+        sin_a = xp.sin(ang)
+        wc = w_row * cos_a
+        c[start:stop] = xp.sum(wc, axis=1)
+        s[start:stop] = xp.sum(w_row * sin_a, axis=1)
+        yc[start:stop] = xp.sum(wy_row * cos_a, axis=1)
+        ys[start:stop] = xp.sum(wy_row * sin_a, axis=1)
+        c2[start:stop] = 2.0 * xp.sum(wc * cos_a, axis=1) - sum_w
+        s2[start:stop] = 2.0 * xp.sum(wc * sin_a, axis=1)
+    return c, s, yc, ys, c2, s2
 
 
 def lombscargle_power_torch(
@@ -267,10 +287,8 @@ def lombscargle_power_torch(
     w_d = to_device_array(w, device=device, dtype=tdtype)
     wy_d = to_device_array(w * y, device=device, dtype=tdtype)
     xp = array_namespace(tau_d)
-    c, s = _trig_sums_direct(xp, tau_d, w_d, f0, df, nf, freq_batch=freq_batch)
-    yc, ys = _trig_sums_direct(xp, tau_d, wy_d, f0, df, nf, freq_batch=freq_batch)
-    c2, s2 = _trig_sums_direct(
-        xp, tau_d, w_d, 2.0 * f0, 2.0 * df, nf, freq_batch=freq_batch
+    c, s, yc, ys, c2, s2 = _gls_sums_direct(
+        xp, tau_d, w_d, wy_d, f0, df, nf, freq_batch=freq_batch
     )
     power = _assemble_power_parts(xp, c, s, yc, ys, c2, s2, y_mean, yy, fit_mean)
     return to_host(power)

From 16c103674d55903f17850a6d9eab252027a690ac Mon Sep 17 00:00:00 2001
From: Tharindu Jayasinghe <tharindujay96@gmail.com>
Date: Wed, 1 Jul 2026 14:12:51 -0700
Subject: [PATCH 4/9] perf(gls): batch the NUFFT base-grid pair; assemble
 one-shot cufinufft on device

The w and w*y transforms share their nonuniform points, so they now run
as one n_trans=2 NUFFT (finufft and cufinufft, plan cache keyed by
(nf, n_trans)). The one-shot cufinufft path keeps all three sums on the
GPU and assembles the power there, so only the final spectrum crosses
back: 194 -> 290 lc/s on 800-point curves; engine path 869 -> 888 lc/s.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
---
 src/cuperiod/methods/gls.py | 111 ++++++++++++++++++++----------------
 1 file changed, 63 insertions(+), 48 deletions(-)

diff --git a/src/cuperiod/methods/gls.py b/src/cuperiod/methods/gls.py
index 26d6e95..452b62f 100644
--- a/src/cuperiod/methods/gls.py
+++ b/src/cuperiod/methods/gls.py
@@ -57,39 +57,47 @@
 
 
 def _trig_sums(
-    tau: FloatArray,
-    strengths: FloatArray,
+    tau: Any,
+    strengths: Any,
     f0: float,
     df: float,
     nf: int,
     backend: NufftBackend,
     eps: float,
-) -> np.ndarray:
-    """``out[k] = sum_j strengths_j exp(2j*pi*(f0 + k*df)*tau_j)`` via a type-1 NUFFT.
-
-    finufft's default mode ordering puts mode ``m = k - nf//2`` at output index ``k``,
-    so the strengths are modulated by ``exp(2j*pi*f_center*tau)`` with
-    ``f_center = f0 + (nf//2)*df`` and the times scaled to ``x = 2*pi*df*tau`` (wrapped
-    to ``[-pi, pi)``); output index ``k`` then lands on frequency ``f0 + k*df``.
+) -> Any:
+    """``out[i, k] = sum_j strengths[i, j] exp(2j*pi*(f0 + k*df)*tau_j)`` via NUFFT.
+
+    ``strengths`` is a ``(n_trans, N)`` stack: transforms that share the same
+    nonuniform points (the base-grid pair ``w`` and ``w*y``) go through **one**
+    batched NUFFT call, sharing the point sort/spread setup. finufft's default mode
+    ordering puts mode ``m = k - nf//2`` at output index ``k``, so the strengths are
+    modulated by ``exp(2j*pi*f_center*tau)`` with ``f_center = f0 + (nf//2)*df`` and
+    the times scaled to ``x = 2*pi*df*tau`` (wrapped to ``[-pi, pi)``); output index
+    ``k`` then lands on frequency ``f0 + k*df``. For ``backend="cufinufft"`` the
+    inputs must already be cupy arrays and the ``(n_trans, nf)`` output **stays on
+    device** so the caller can assemble the power there.
     """
-    f_center = f0 + (nf // 2) * df
-    x = (2.0 * np.pi * df) * tau
-    x = np.mod(x + np.pi, 2.0 * np.pi) - np.pi
-    c = (strengths * np.exp(2j * np.pi * f_center * tau)).astype(np.complex128)
-
     if backend == "finufft":
         import finufft
 
-        out = finufft.nufft1d1(x, c, nf, eps=eps, isign=1)
-        return np.asarray(out, dtype=np.complex128)
-    if backend == "cufinufft":
+        xp: Any = np
+        nufft1d1 = finufft.nufft1d1
+    elif backend == "cufinufft":
         ensure_cuda_dll_path()
         import cufinufft
         import cupy as cp
 
-        out_g = cufinufft.nufft1d1(cp.asarray(x), cp.asarray(c), nf, eps=eps, isign=1)
-        return np.asarray(cp.asnumpy(out_g), dtype=np.complex128)
-    raise ValueError(f"unknown NUFFT backend {backend!r}")
+        xp = cp
+        nufft1d1 = cufinufft.nufft1d1
+    else:
+        raise ValueError(f"unknown NUFFT backend {backend!r}")
+
+    f_center = f0 + (nf // 2) * df
+    x = (2.0 * np.pi * df) * tau
+    x = xp.mod(x + np.pi, 2.0 * np.pi) - np.pi
+    c = (strengths * xp.exp(2j * np.pi * f_center * tau)[None, :]).astype(np.complex128)
+    out = nufft1d1(x, c, nf, eps=eps, isign=1)
+    return out if out.ndim == 2 else out[None, :]
 
 
 def _assemble_power_parts(
@@ -201,12 +209,19 @@ def lombscargle_power(
     if nf <= 0:
         return np.zeros(0, dtype=np.float64)
     tau, w, y, y_mean, yy = _prep(t, y, dy)
-    sw = _trig_sums(tau, w, f0, df, nf, backend, eps)
-    swy = _trig_sums(tau, w * y, f0, df, nf, backend, eps)
-    sw2 = _trig_sums(tau, w, 2.0 * f0, 2.0 * df, nf, backend, eps)
-    return np.asarray(
-        _assemble_power(sw, swy, sw2, y_mean, yy, fit_mean), dtype=np.float64
-    )
+    base = np.stack([w, w * y])
+    tau_b: Any = tau
+    if backend == "cufinufft":
+        # Device inputs in, device sums out: the power is assembled on the GPU and
+        # only the final spectrum crosses back to the host.
+        ensure_cuda_dll_path()
+        import cupy as cp
+
+        tau_b = cp.asarray(tau)
+        base = cp.asarray(base)
+    pair = _trig_sums(tau_b, base, f0, df, nf, backend, eps)
+    sw2 = _trig_sums(tau_b, base[:1], 2.0 * f0, 2.0 * df, nf, backend, eps)
+    return to_host(_assemble_power(pair[0], pair[1], sw2[0], y_mean, yy, fit_mean))
 
 
 #: Elements per transient ``(chunk, N)`` matrix in the direct path (64 MB float64);
@@ -326,21 +341,17 @@ def __init__(self, eps: float = DEFAULT_EPS, nf_bucket: int = 1 << 16) -> None:
         self._cp = cupy
         self._eps = eps
         self._bucket = max(1, int(nf_bucket))
-        self._plans: dict[int, Any] = {}
+        self._plans: dict[tuple[int, int], Any] = {}
 
-    def _plan(self, nf: int) -> Any:
-        plan = self._plans.get(nf)
+    def _plan(self, nf: int, n_trans: int) -> Any:
+        plan = self._plans.get((nf, n_trans))
         if plan is None:
             plan = self._cufinufft.Plan(
-                1, (nf,), eps=self._eps, isign=1, dtype="complex128"
+                1, (nf,), n_trans=n_trans, eps=self._eps, isign=1, dtype="complex128"
             )
-            self._plans[nf] = plan
+            self._plans[(nf, n_trans)] = plan
         return plan
 
-    @staticmethod
-    def _squeeze(arr: Any) -> Any:
-        return arr[0] if getattr(arr, "ndim", 1) == 2 else arr
-
     def power(
         self,
         t: FloatArray,
@@ -352,31 +363,35 @@ def power(
         *,
         fit_mean: bool = True,
     ) -> FloatArray:
-        """GLS power on ``f0 + df*arange(nf)`` via a reused cufinufft plan."""
+        """GLS power on ``f0 + df*arange(nf)`` via reused cufinufft plans.
+
+        The base-grid pair (``w`` and ``w*y`` share the same points) runs as one
+        ``n_trans=2`` batched transform — one ``setpts``/execute instead of two — and
+        the doubled-grid sum uses an ``n_trans=1`` plan of the same bucketed size.
+        """
         if nf <= 0:
             return np.zeros(0, dtype=np.float64)
         cp = self._cp
         tau, w, y, y_mean, yy = _prep(t, y, dy)
         tau_g = cp.asarray(tau)
-        wg = cp.asarray(w)
-        wyg = cp.asarray(w * y)
+        base = cp.asarray(np.stack([w, w * y]))
         nf_plan = ((nf + self._bucket - 1) // self._bucket) * self._bucket
-        plan = self._plan(nf_plan)
         two_pi = 2.0 * np.pi
 
-        def sums(f0_: float, df_: float, strengths: tuple[Any, ...]) -> list[Any]:
+        def sums(f0_: float, df_: float, strengths: Any) -> Any:
             x = (two_pi * df_) * tau_g
             x = cp.mod(x + np.pi, two_pi) - np.pi
             mod = cp.exp(2j * np.pi * (f0_ + (nf_plan // 2) * df_) * tau_g)
+            plan = self._plan(nf_plan, int(strengths.shape[0]))
             plan.setpts(x)
-            return [
-                self._squeeze(plan.execute((c * mod).astype(cp.complex128)))
-                for c in strengths
-            ]
-
-        sw, swy = sums(f0, df, (wg, wyg))
-        (sw2,) = sums(2.0 * f0, 2.0 * df, (wg,))
-        power = _assemble_power(sw[:nf], swy[:nf], sw2[:nf], y_mean, yy, fit_mean)
+            out = plan.execute((strengths * mod[None, :]).astype(cp.complex128))
+            return out if out.ndim == 2 else out[None, :]
+
+        pair = sums(f0, df, base)
+        sw2 = sums(2.0 * f0, 2.0 * df, base[:1])
+        power = _assemble_power(
+            pair[0, :nf], pair[1, :nf], sw2[0, :nf], y_mean, yy, fit_mean
+        )
         return np.asarray(cp.asnumpy(power), dtype=np.float64)
 
 

From 559fcdf27b072a8fc23193fa6a8665fe23439e95 Mon Sep 17 00:00:00 2001
From: Tharindu Jayasinghe <tharindujay96@gmail.com>
Date: Wed, 1 Jul 2026 14:21:27 -0700
Subject: [PATCH 5/9] perf: extend the multicore numba CPU tier to PDM, CE, SL,
 MHAOV, and TLS

One-iteration-per-trial prange kernels mirroring the CUDA designs (PDM
uses the fine-bin regroup; SL uses a stable mergesort argsort so tied
phases pair identically to the array-API backends; MHAOV runs the
Chebyshev harmonic recurrence in scalar registers with a per-frequency
d x d solve). cpu/auto now resolve to numba when the [fast] extra is
installed via a shared fast_cpu_backend attribute, which also replaces
the bespoke BLS resolve_backend override.

3k points, warm kernels, vs the vectorized numpy paths: PDM 306x, CE
135x, MHAOV 57x, TLS 53x, SL 23x; parity <= 2.5e-11 (TLS bit-identical).

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
---
 src/cuperiod/core/config.py                 |  32 +++---
 src/cuperiod/methods/base.py                |  16 ++-
 src/cuperiod/methods/bls.py                 |  25 +----
 src/cuperiod/methods/conditional_entropy.py |  64 ++++++++++-
 src/cuperiod/methods/mhaov.py               | 111 ++++++++++++++++++-
 src/cuperiod/methods/pdm.py                 |  85 ++++++++++++++-
 src/cuperiod/methods/string_length.py       |  60 ++++++++++-
 src/cuperiod/methods/tls.py                 | 113 +++++++++++++++++++-
 8 files changed, 451 insertions(+), 55 deletions(-)

diff --git a/src/cuperiod/core/config.py b/src/cuperiod/core/config.py
index 0cd78d9..135d7ea 100644
--- a/src/cuperiod/core/config.py
+++ b/src/cuperiod/core/config.py
@@ -165,7 +165,7 @@ def _check_bounds(self) -> Self:
         default=20, ge=3, description="Skip if fewer finite points."
     )
     backend: Literal[
-        "auto", "cpu", "gpu", "numpy", "astropy", "cupy", "torch"
+        "auto", "cpu", "gpu", "numba", "numpy", "astropy", "cupy", "torch"
     ] = Field(default="auto", description="Compute backend.")
     batch_periods: int = Field(
         default=2048, ge=1, description="Trial periods per vectorized batch."
@@ -213,9 +213,9 @@ def _check_bounds(self) -> Self:
     min_detections: int = Field(
         default=20, ge=3, description="Skip if fewer finite points."
     )
-    backend: Literal["auto", "cpu", "gpu", "numpy", "cupy", "torch"] = Field(
-        default="auto", description="Compute backend."
-    )
+    backend: Literal[
+        "auto", "cpu", "gpu", "numba", "numpy", "cupy", "torch"
+    ] = Field(default="auto", description="Compute backend.")
     batch_periods: int = Field(
         default=2048, ge=1, description="Trial periods per vectorized batch."
     )
@@ -261,9 +261,9 @@ def _check_bounds(self) -> Self:
     min_detections: int = Field(
         default=20, ge=5, description="Skip if fewer finite points (need > 2H+1)."
     )
-    backend: Literal["auto", "cpu", "gpu", "numpy", "cupy", "torch"] = Field(
-        default="auto", description="Compute backend."
-    )
+    backend: Literal[
+        "auto", "cpu", "gpu", "numba", "numpy", "cupy", "torch"
+    ] = Field(default="auto", description="Compute backend.")
     batch_periods: int = Field(
         default=512, ge=1, description="Trial frequencies per vectorized batch."
     )
@@ -308,9 +308,9 @@ def _check_bounds(self) -> Self:
     min_detections: int = Field(
         default=20, ge=3, description="Skip if fewer finite points."
     )
-    backend: Literal["auto", "cpu", "gpu", "numpy", "cupy", "torch"] = Field(
-        default="auto", description="Compute backend."
-    )
+    backend: Literal[
+        "auto", "cpu", "gpu", "numba", "numpy", "cupy", "torch"
+    ] = Field(default="auto", description="Compute backend.")
     batch_periods: int = Field(
         default=1024, ge=1, description="Trial periods per vectorized batch."
     )
@@ -353,9 +353,9 @@ def _check_bounds(self) -> Self:
     min_detections: int = Field(
         default=20, ge=3, description="Skip if fewer finite points."
     )
-    backend: Literal["auto", "cpu", "gpu", "numpy", "cupy", "torch"] = Field(
-        default="auto", description="Compute backend."
-    )
+    backend: Literal[
+        "auto", "cpu", "gpu", "numba", "numpy", "cupy", "torch"
+    ] = Field(default="auto", description="Compute backend.")
     batch_periods: int = Field(
         default=1024, ge=1, description="Trial periods per vectorized batch."
     )
@@ -425,9 +425,9 @@ def _check_bounds(self) -> Self:
     min_detections: int = Field(
         default=20, ge=3, description="Skip if fewer finite points."
     )
-    backend: Literal["auto", "cpu", "gpu", "numpy", "cupy", "torch"] = Field(
-        default="auto", description="Compute backend."
-    )
+    backend: Literal[
+        "auto", "cpu", "gpu", "numba", "numpy", "cupy", "torch"
+    ] = Field(default="auto", description="Compute backend.")
     period_batch: int = Field(
         default=256, ge=1, description="Trial periods per vectorized batch."
     )
diff --git a/src/cuperiod/methods/base.py b/src/cuperiod/methods/base.py
index 653190d..65b6b4e 100644
--- a/src/cuperiod/methods/base.py
+++ b/src/cuperiod/methods/base.py
@@ -52,8 +52,12 @@ class PeriodogramMethod(ABC):
     natural_domain: ClassVar[Domain | None] = None
     #: The settings model class for this method.
     settings_cls: ClassVar[type[BaseSettings]]
-    #: Best CPU backend name.
+    #: Best always-available CPU backend name.
     cpu_backend: ClassVar[str]
+    #: Preferred CPU backend when its optional dependency is installed (the multicore
+    #: ``numba`` kernels of the ``[fast]`` extra), or ``None`` to always use
+    #: ``cpu_backend``.
+    fast_cpu_backend: ClassVar[str | None] = None
     #: NVIDIA fast-path GPU backend name (cufinufft / cupy), or ``None`` if the method
     #: has no CUDA path.
     gpu_backend: ClassVar[str | None] = None
@@ -108,9 +112,9 @@ def resolve_backend(self, requested: str) -> str:
                 return self.gpu_backend
             if self.portable_gpu_backend is not None and torch_gpu_available():
                 return self.portable_gpu_backend
-            return self.cpu_backend
+            return self._best_cpu_backend(available)
         if requested == "cpu":
-            return self.cpu_backend
+            return self._best_cpu_backend(available)
         if requested == "gpu":
             if self.gpu_backend is not None and cuda_available():
                 return self.gpu_backend
@@ -139,6 +143,12 @@ def resolve_backend(self, requested: str) -> str:
             )
         return requested
 
+    def _best_cpu_backend(self, available: set[str]) -> str:
+        """``fast_cpu_backend`` when its dependency is installed, else ``cpu_backend``."""
+        if self.fast_cpu_backend is not None and self.fast_cpu_backend in available:
+            return self.fast_cpu_backend
+        return self.cpu_backend
+
     def _resolve_torch(self, requested: str) -> str:
         """Validate a concrete ``torch``/``torch:<device>`` request for this method."""
         if self.portable_gpu_backend != "torch":
diff --git a/src/cuperiod/methods/bls.py b/src/cuperiod/methods/bls.py
index abb308c..54ab348 100644
--- a/src/cuperiod/methods/bls.py
+++ b/src/cuperiod/methods/bls.py
@@ -154,36 +154,13 @@ class BLSMethod(PeriodogramMethod):
     natural_domain: ClassVar[Domain] = Domain.FLUX
     settings_cls: ClassVar[type] = BLSSettings
     cpu_backend: ClassVar[str] = "astropy"
+    fast_cpu_backend: ClassVar[str | None] = "numba"
     gpu_backend: ClassVar[str | None] = "cupy"
     portable_gpu_backend: ClassVar[str | None] = "torch"
     all_backends: ClassVar[tuple[str, ...]] = (
         "numba", "numpy", "astropy", "cupy", "torch",
     )
 
-    def resolve_backend(self, requested: str) -> str:
-        """Prefer cupy on NVIDIA, then torch on other GPUs, else the numba CPU search.
-
-        Keeps BLS's CPU preference (multicore numba when installed, else astropy) for
-        ``cpu``/``auto``, while ``auto`` still reaches a GPU: the cupy kernel on CUDA,
-        then the portable torch path on AMD/Intel/Mac. Concrete ``torch``/``torch:*``
-        requests are validated by the base method.
-        """
-        from cuperiod.core.backend import (
-            available_backends,
-            cuda_available,
-            torch_gpu_available,
-        )
-
-        if requested == "auto":
-            if self.gpu_backend is not None and cuda_available():
-                return self.gpu_backend
-            if self.portable_gpu_backend is not None and torch_gpu_available():
-                return self.portable_gpu_backend
-            return "numba" if "numba" in available_backends() else "astropy"
-        if requested == "cpu":
-            return "numba" if "numba" in available_backends() else "astropy"
-        return super().resolve_backend(requested)
-
     def default_grid(self, lc: LightCurve, settings: BLSSettings) -> GridSpec:  # type: ignore[override]
         finite = lc.finite()
         grids = _segment_grids(finite.baseline, settings)
diff --git a/src/cuperiod/methods/conditional_entropy.py b/src/cuperiod/methods/conditional_entropy.py
index 4408751..573612d 100644
--- a/src/cuperiod/methods/conditional_entropy.py
+++ b/src/cuperiod/methods/conditional_entropy.py
@@ -193,6 +193,59 @@ def _ce_cuda(
     return np.asarray(cp.asnumpy(out), dtype=np.float64)
 
 
+# --- CPU fast path: numba-parallel, one loop-iteration per trial period -------
+
+_NUMBA_CE_KERNEL: Any = None
+
+
+def _numba_ce_kernel() -> Any:
+    """Lazily compile (once) and cache the numba CE kernel.
+
+    The same fold-bin-reduce as the CUDA kernel, JIT-run across all cores with
+    ``prange``. The phase is ``q - floor(q)``, the formula of the vectorized numpy
+    path, so the two CPU backends bin identically.
+    """
+    global _NUMBA_CE_KERNEL
+    if _NUMBA_CE_KERNEL is not None:
+        return _NUMBA_CE_KERNEL
+    from numba import njit, prange
+
+    @njit(parallel=True, cache=True, fastmath=False)  # pragma: no cover - njit
+    def _kernel(tau, mag_bin, periods, n_phase, n_mag):  # type: ignore[no-untyped-def]
+        n_periods = periods.shape[0]
+        n_points = tau.shape[0]
+        n_cells = n_phase * n_mag
+        out = np.empty(n_periods)
+        for pidx in prange(n_periods):
+            period = periods[pidx]
+            cnt = np.zeros(n_cells)
+            for j in range(n_points):
+                q = tau[j] / period
+                ph = q - np.floor(q)
+                pb = int(ph * n_phase)
+                if pb > n_phase - 1:
+                    pb = n_phase - 1
+                elif pb < 0:
+                    pb = 0
+                cnt[pb * n_mag + mag_bin[j]] += 1.0
+            h = 0.0
+            for p in range(n_phase):
+                ci = 0.0
+                for m in range(n_mag):
+                    ci += cnt[p * n_mag + m]
+                if ci > 0.0:
+                    lci = np.log(ci)
+                    for m in range(n_mag):
+                        c = cnt[p * n_mag + m]
+                        if c > 0.0:
+                            h += c * (lci - np.log(c))
+            out[pidx] = h / n_points
+        return out
+
+    _NUMBA_CE_KERNEL = _kernel
+    return _kernel
+
+
 def conditional_entropy(
     t: FloatArray,
     y: FloatArray,
@@ -242,6 +295,12 @@ def conditional_entropy(
         return _ce_cuda(
             tau, mag_bin, periods_host, n_phase=n_phase_bins, n_mag=n_mag_bins
         )
+    if backend == "numba":
+        kernel = _numba_ce_kernel()
+        return np.asarray(
+            kernel(tau, mag_bin, periods_host, n_phase_bins, n_mag_bins),
+            dtype=np.float64,
+        )
     if backend == "torch" or backend.startswith("torch:"):
         import torch
 
@@ -264,16 +323,17 @@ def conditional_entropy(
 
 
 class ConditionalEntropyMethod(PeriodogramMethod):
-    """Conditional-entropy period search (numpy CPU, cupy GPU)."""
+    """Conditional-entropy period search (numba/numpy CPU, cupy GPU)."""
 
     name: ClassVar[str] = "CE"
     objective_sense: ClassVar[Literal["max", "min"]] = "min"
     supports_multiband: ClassVar[bool] = False
     settings_cls: ClassVar[type] = CESettings
     cpu_backend: ClassVar[str] = "numpy"
+    fast_cpu_backend: ClassVar[str | None] = "numba"
     gpu_backend: ClassVar[str | None] = "cupy"
     portable_gpu_backend: ClassVar[str | None] = "torch"
-    all_backends: ClassVar[tuple[str, ...]] = ("numpy", "cupy", "torch")
+    all_backends: ClassVar[tuple[str, ...]] = ("numba", "numpy", "cupy", "torch")
 
     def default_grid(self, lc: LightCurve, settings: CESettings) -> GridSpec:  # type: ignore[override]
         finite = lc.finite()
diff --git a/src/cuperiod/methods/mhaov.py b/src/cuperiod/methods/mhaov.py
index 2d3474c..ff7849f 100644
--- a/src/cuperiod/methods/mhaov.py
+++ b/src/cuperiod/methods/mhaov.py
@@ -183,6 +183,105 @@ def _model_ss_batch(
     return out
 
 
+# --- CPU fast path: numba-parallel, one loop-iteration per trial frequency ----
+
+_NUMBA_MHAOV_KERNEL: Any = None
+
+
+def _numba_mhaov_kernel() -> Any:
+    """Lazily compile (once) and cache the numba MHAOV kernel.
+
+    The same harmonic-trig-sum normal equations as :func:`_model_ss_batch`, with the
+    Chebyshev recurrence run per point in scalar registers — no ``(F, N)`` transients
+    at all — and the tiny ``d×d`` solve done per frequency. ``prange`` over the
+    frequency grid.
+    """
+    global _NUMBA_MHAOV_KERNEL
+    if _NUMBA_MHAOV_KERNEL is not None:
+        return _NUMBA_MHAOV_KERNEL
+    from numba import njit, prange
+
+    @njit(parallel=True, cache=True, fastmath=False)  # pragma: no cover - njit
+    def _kernel(tau, y, freqs, n_harmonics, ridge, sum_y, y_mean, total_ss):  # type: ignore[no-untyped-def]
+        n_freq = freqs.shape[0]
+        n_points = tau.shape[0]
+        h = n_harmonics
+        m_max = 2 * h
+        d = 2 * h + 1
+        two_pi = 2.0 * np.pi
+        out = np.empty(n_freq)
+        for fi in prange(n_freq):
+            f = freqs[fi]
+            c_sums = np.zeros(m_max + 1)
+            s_sums = np.zeros(m_max + 1)
+            a_sums = np.zeros(h + 1)
+            b_sums = np.zeros(h + 1)
+            for j in range(n_points):
+                theta = two_pi * f * tau[j]
+                c1 = np.cos(theta)
+                s1 = np.sin(theta)
+                two_c1 = 2.0 * c1
+                cm = c1
+                sm = s1
+                c_prev = 1.0
+                s_prev = 0.0
+                yj = y[j]
+                for m in range(1, m_max + 1):
+                    c_sums[m] += cm
+                    s_sums[m] += sm
+                    if m <= h:
+                        a_sums[m] += yj * cm
+                        b_sums[m] += yj * sm
+                    c_next = two_c1 * cm - c_prev
+                    s_next = two_c1 * sm - s_prev
+                    c_prev = cm
+                    s_prev = sm
+                    cm = c_next
+                    sm = s_next
+            gram = np.zeros((d, d))
+            proj = np.zeros(d)
+            gram[0, 0] = n_points + ridge
+            proj[0] = sum_y
+            for k in range(1, h + 1):
+                gram[0, 2 * k - 1] = c_sums[k]
+                gram[2 * k - 1, 0] = c_sums[k]
+                gram[0, 2 * k] = s_sums[k]
+                gram[2 * k, 0] = s_sums[k]
+                proj[2 * k - 1] = a_sums[k]
+                proj[2 * k] = b_sums[k]
+                for line in range(1, h + 1):
+                    m_diff = k - line if k >= line else line - k
+                    c_diff = c_sums[m_diff] if m_diff else float(n_points)
+                    cc = 0.5 * (c_diff + c_sums[k + line])
+                    ss = 0.5 * (c_diff - c_sums[k + line])
+                    if k == line:
+                        s_diff = 0.0
+                    elif k > line:
+                        s_diff = s_sums[m_diff]
+                    else:
+                        s_diff = -s_sums[m_diff]
+                    cs = 0.5 * (s_sums[k + line] - s_diff)
+                    extra = ridge if k == line else 0.0
+                    gram[2 * k - 1, 2 * line - 1] = cc + extra
+                    gram[2 * k, 2 * line] = ss + extra
+                    gram[2 * k - 1, 2 * line] = cs
+                    gram[2 * line, 2 * k - 1] = cs
+            beta = np.linalg.solve(gram, proj)
+            model_ss = 0.0
+            for i in range(d):
+                model_ss += beta[i] * proj[i]
+            model_ss -= n_points * y_mean * y_mean
+            if model_ss < 0.0:
+                model_ss = 0.0
+            elif model_ss > total_ss:
+                model_ss = total_ss
+            out[fi] = model_ss
+        return out
+
+    _NUMBA_MHAOV_KERNEL = _kernel
+    return _kernel
+
+
 def _compute_model_ss(
     t: FloatArray,
     y: FloatArray,
@@ -210,6 +309,13 @@ def _compute_model_ss(
     if total_ss <= 0.0:
         return np.zeros(freqs.size, dtype=np.float64), 0.0, n
 
+    if backend == "numba":
+        kernel = _numba_mhaov_kernel()
+        ridge = _RIDGE_EPS * float(np.finfo(np.float64).eps) * float(n)
+        out = kernel(
+            tau, y, freqs, n_harmonics, ridge, float(y.sum()), y_mean, total_ss
+        )
+        return np.asarray(out, dtype=np.float64), total_ss, n
     if backend == "cupy":
         ensure_cuda_dll_path()
         import cupy as cp
@@ -349,16 +455,17 @@ def aov_multiband_power(
 
 
 class MHAOVMethod(PeriodogramMethod):
-    """Multiharmonic Analysis of Variance (numpy CPU, cupy GPU)."""
+    """Multiharmonic Analysis of Variance (numba/numpy CPU, cupy GPU)."""
 
     name: ClassVar[str] = "MHAOV"
     objective_sense: ClassVar[Literal["max", "min"]] = "max"
     supports_multiband: ClassVar[bool] = True
     settings_cls: ClassVar[type] = MHAOVSettings
     cpu_backend: ClassVar[str] = "numpy"
+    fast_cpu_backend: ClassVar[str | None] = "numba"
     gpu_backend: ClassVar[str | None] = "cupy"
     portable_gpu_backend: ClassVar[str | None] = "torch"
-    all_backends: ClassVar[tuple[str, ...]] = ("numpy", "cupy", "torch")
+    all_backends: ClassVar[tuple[str, ...]] = ("numba", "numpy", "cupy", "torch")
 
     def default_grid(self, lc: LightCurve, settings: MHAOVSettings) -> GridSpec:  # type: ignore[override]
         finite = lc.finite()
diff --git a/src/cuperiod/methods/pdm.py b/src/cuperiod/methods/pdm.py
index ebbc1e8..1cad0ed 100644
--- a/src/cuperiod/methods/pdm.py
+++ b/src/cuperiod/methods/pdm.py
@@ -249,6 +249,76 @@ def _pdm_cuda(
     return np.asarray(cp.asnumpy(out), dtype=np.float64)
 
 
+# --- CPU fast path: numba-parallel, one loop-iteration per trial period -------
+
+_NUMBA_PDM_KERNEL: Any = None
+
+
+def _numba_pdm_kernel() -> Any:
+    """Lazily compile (once) and cache the numba PDM kernel.
+
+    The same fine-bin search as the CUDA kernel — bin each folded point once into
+    ``n_bins*n_covers`` fine bins, regroup covers exactly in the reduction — JIT-run
+    across all cores with ``prange`` (each trial period is independent). The phase is
+    ``q - floor(q)``, the formula of the vectorized numpy path, so the two CPU
+    backends bin identically.
+    """
+    global _NUMBA_PDM_KERNEL
+    if _NUMBA_PDM_KERNEL is not None:
+        return _NUMBA_PDM_KERNEL
+    from numba import njit, prange
+
+    @njit(parallel=True, cache=True, fastmath=False)  # pragma: no cover - njit
+    def _kernel(tau, y, periods, n_bins, n_covers, sigma2):  # type: ignore[no-untyped-def]
+        n_periods = periods.shape[0]
+        n_points = tau.shape[0]
+        m_fine = n_bins * n_covers
+        out = np.empty(n_periods)
+        for pidx in prange(n_periods):
+            period = periods[pidx]
+            s_sum = np.zeros(m_fine)
+            s_sq = np.zeros(m_fine)
+            s_cnt = np.zeros(m_fine)
+            for j in range(n_points):
+                q = tau[j] / period
+                ph = q - np.floor(q)
+                f = int(ph * m_fine)
+                if f > m_fine - 1:
+                    f = m_fine - 1
+                elif f < 0:
+                    f = 0
+                yj = y[j]
+                s_sum[f] += yj
+                s_sq[f] += yj * yj
+                s_cnt[f] += 1.0
+            ssd = 0.0
+            nonempty = 0
+            for c in range(n_covers):
+                for b in range(n_bins):
+                    cnt = 0.0
+                    su = 0.0
+                    sq = 0.0
+                    f0 = b * n_covers - c
+                    if f0 < 0:
+                        f0 += m_fine
+                    for k in range(n_covers):
+                        f = f0 + k
+                        if f >= m_fine:
+                            f -= m_fine
+                        cnt += s_cnt[f]
+                        su += s_sum[f]
+                        sq += s_sq[f]
+                    if cnt > 0.0:
+                        ssd += sq - su * su / cnt
+                        nonempty += 1
+            den = float(n_points * n_covers - nonempty)
+            out[pidx] = (ssd / den) / sigma2 if den > 0.0 else np.inf
+        return out
+
+    _NUMBA_PDM_KERNEL = _kernel
+    return _kernel
+
+
 def pdm_theta(
     t: FloatArray,
     y: FloatArray,
@@ -272,8 +342,9 @@ def pdm_theta(
         Phase bins per cover.
     n_covers : int, default 3
         Overlapping bin sets, offset by ``1/(n_bins*n_covers)`` in phase.
-    backend : {"numpy", "cupy"}, default "numpy"
-        CPU or GPU.
+    backend : str, default "numpy"
+        ``"numpy"`` (vectorized CPU), ``"numba"`` (multicore CPU), ``"cupy"``
+        (NVIDIA RawKernel), or ``"torch"`` / ``"torch:<device>"`` (portable).
     batch : int, default 2048
         Trial periods per vectorized batch.
 
@@ -300,6 +371,11 @@ def pdm_theta(
         return _pdm_cuda(
             tau, y, periods_host, n_bins=n_bins, n_covers=n_covers, sigma2=sigma2
         )
+    if backend == "numba":
+        kernel = _numba_pdm_kernel()
+        return np.asarray(
+            kernel(tau, y, periods_host, n_bins, n_covers, sigma2), dtype=np.float64
+        )
     if backend == "torch" or backend.startswith("torch:"):
         import torch
 
@@ -324,16 +400,17 @@ def pdm_theta(
 
 
 class PDMMethod(PeriodogramMethod):
-    """Phase Dispersion Minimization (numpy CPU, cupy GPU)."""
+    """Phase Dispersion Minimization (numba/numpy CPU, cupy GPU)."""
 
     name: ClassVar[str] = "PDM"
     objective_sense: ClassVar[Literal["max", "min"]] = "min"
     supports_multiband: ClassVar[bool] = False
     settings_cls: ClassVar[type] = PDMSettings
     cpu_backend: ClassVar[str] = "numpy"
+    fast_cpu_backend: ClassVar[str | None] = "numba"
     gpu_backend: ClassVar[str | None] = "cupy"
     portable_gpu_backend: ClassVar[str | None] = "torch"
-    all_backends: ClassVar[tuple[str, ...]] = ("numpy", "cupy", "torch")
+    all_backends: ClassVar[tuple[str, ...]] = ("numba", "numpy", "cupy", "torch")
 
     def default_grid(self, lc: LightCurve, settings: PDMSettings) -> GridSpec:  # type: ignore[override]
         finite = lc.finite()
diff --git a/src/cuperiod/methods/string_length.py b/src/cuperiod/methods/string_length.py
index 078ec02..cec36d5 100644
--- a/src/cuperiod/methods/string_length.py
+++ b/src/cuperiod/methods/string_length.py
@@ -82,6 +82,58 @@ def _length_batch(
     return length
 
 
+# --- CPU fast path: numba-parallel, one loop-iteration per trial period -------
+
+_NUMBA_SL_KERNEL: Any = None
+
+
+def _numba_sl_kernel() -> Any:
+    """Lazily compile (once) and cache the numba string-length kernel.
+
+    Per period: fold (``q - floor(q)``, the vectorized path's formula), **stable**
+    ``mergesort`` argsort (matching the array-API paths' stable sort, so equal phases
+    keep the same backend-independent neighbour pairing), and accumulate the string.
+    ``prange`` over periods; no ``(P, N)`` transients at all.
+    """
+    global _NUMBA_SL_KERNEL
+    if _NUMBA_SL_KERNEL is not None:
+        return _NUMBA_SL_KERNEL
+    from numba import njit, prange
+
+    @njit(parallel=True, cache=True, fastmath=False)  # pragma: no cover - njit
+    def _kernel(tau, m_scaled, periods):  # type: ignore[no-untyped-def]
+        n_periods = periods.shape[0]
+        n_points = tau.shape[0]
+        out = np.empty(n_periods)
+        for pidx in prange(n_periods):
+            period = periods[pidx]
+            phase = np.empty(n_points)
+            for j in range(n_points):
+                q = tau[j] / period
+                phase[j] = q - np.floor(q)
+            order = np.argsort(phase, kind="mergesort")
+            first = order[0]
+            prev_ph = phase[first]
+            prev_m = m_scaled[first]
+            total = 0.0
+            for j in range(1, n_points):
+                idx = order[j]
+                ph = phase[idx]
+                mm = m_scaled[idx]
+                dphi = ph - prev_ph
+                dmag = mm - prev_m
+                total += np.sqrt(dphi * dphi + dmag * dmag)
+                prev_ph = ph
+                prev_m = mm
+            wrap_phi = (phase[first] + 1.0) - prev_ph
+            wrap_mag = m_scaled[first] - prev_m
+            out[pidx] = total + np.sqrt(wrap_phi * wrap_phi + wrap_mag * wrap_mag)
+        return out
+
+    _NUMBA_SL_KERNEL = _kernel
+    return _kernel
+
+
 def string_length(
     t: FloatArray,
     y: FloatArray,
@@ -128,6 +180,9 @@ def string_length(
             batch=batch,
         )
         return np.asarray(cp.asnumpy(length), dtype=np.float64)
+    if backend == "numba":
+        kernel = _numba_sl_kernel()
+        return np.asarray(kernel(tau, m_scaled, periods_host), dtype=np.float64)
     if backend == "torch" or backend.startswith("torch:"):
         import torch
 
@@ -151,16 +206,17 @@ def string_length(
 
 
 class StringLengthMethod(PeriodogramMethod):
-    """String-length period search (numpy CPU, cupy GPU)."""
+    """String-length period search (numba/numpy CPU, cupy GPU)."""
 
     name: ClassVar[str] = "STRINGLENGTH"
     objective_sense: ClassVar[Literal["max", "min"]] = "min"
     supports_multiband: ClassVar[bool] = False
     settings_cls: ClassVar[type] = StringLengthSettings
     cpu_backend: ClassVar[str] = "numpy"
+    fast_cpu_backend: ClassVar[str | None] = "numba"
     gpu_backend: ClassVar[str | None] = "cupy"
     portable_gpu_backend: ClassVar[str | None] = "torch"
-    all_backends: ClassVar[tuple[str, ...]] = ("numpy", "cupy", "torch")
+    all_backends: ClassVar[tuple[str, ...]] = ("numba", "numpy", "cupy", "torch")
 
     def default_grid(self, lc: LightCurve, settings: StringLengthSettings) -> GridSpec:  # type: ignore[override]
         finite = lc.finite()
diff --git a/src/cuperiod/methods/tls.py b/src/cuperiod/methods/tls.py
index 9135baf..7456344 100644
--- a/src/cuperiod/methods/tls.py
+++ b/src/cuperiod/methods/tls.py
@@ -343,6 +343,109 @@ def _tls_cuda(
     return {name: np.asarray(cp.asnumpy(out[name]), dtype=np.float64) for name in names}
 
 
+# --- CPU fast path: numba-parallel, one loop-iteration per trial period -------
+
+_NUMBA_TLS_KERNEL: Any = None
+
+
+def _numba_tls_kernel() -> Any:
+    """Lazily compile (once) and cache the numba TLS kernel.
+
+    The same fold-bin-correlate sweep as the CUDA kernel — bin the weighted residuals
+    into phase bins, slide every limb-darkened template over every start phase, keep
+    the best signal residue — JIT-run across all cores with ``prange``. The duration
+    then start iteration order and strict ``>`` improvement reproduce the vectorized
+    path's tie-breaking exactly.
+    """
+    global _NUMBA_TLS_KERNEL
+    if _NUMBA_TLS_KERNEL is not None:
+        return _NUMBA_TLS_KERNEL
+    from numba import njit, prange
+
+    @njit(parallel=True, cache=True, fastmath=False)  # pragma: no cover - njit
+    def _kernel(tau, yw, wv, periods, dur_bins, tmpl_off, tmpl, n_bins, w_eps):  # type: ignore[no-untyped-def]
+        n_periods = periods.shape[0]
+        n_points = tau.shape[0]
+        n_dur = dur_bins.shape[0]
+        o_sr = np.zeros(n_periods)
+        o_depth = np.zeros(n_periods)
+        o_duration = np.zeros(n_periods)
+        o_t0 = np.zeros(n_periods)
+        for pidx in prange(n_periods):
+            period = periods[pidx]
+            a = np.zeros(n_bins)
+            b = np.zeros(n_bins)
+            for j in range(n_points):
+                q = tau[j] / period
+                ph = q - np.floor(q)
+                bb = int(ph * n_bins)
+                if bb > n_bins - 1:
+                    bb = n_bins - 1
+                elif bb < 0:
+                    bb = 0
+                a[bb] += yw[j]
+                b[bb] += wv[j]
+            loc_sr = 0.0
+            loc_depth = 0.0
+            loc_start = 0
+            loc_width = 0
+            for di in range(n_dur):
+                width = dur_bins[di]
+                off = tmpl_off[di]
+                for start in range(n_bins):
+                    num = 0.0
+                    den = 0.0
+                    for k in range(width):
+                        bb = start + k
+                        if bb >= n_bins:
+                            bb -= n_bins  # circular wrap
+                        gk = tmpl[off + k]
+                        num += gk * a[bb]
+                        den += gk * gk * b[bb]
+                    if den > w_eps and num < 0.0:
+                        sr = num * num / den
+                        if sr > loc_sr:
+                            loc_sr = sr
+                            loc_start = start
+                            loc_width = width
+                            loc_depth = -num / den
+            o_sr[pidx] = loc_sr
+            o_depth[pidx] = loc_depth
+            o_duration[pidx] = loc_width / n_bins * period
+            centre = (loc_start + loc_width / 2.0) / n_bins
+            o_t0[pidx] = (centre - np.floor(centre)) * period
+        return o_sr, o_depth, o_duration, o_t0
+
+    _NUMBA_TLS_KERNEL = _kernel
+    return _kernel
+
+
+def _tls_numba(
+    tau: FloatArray,
+    yw: FloatArray,
+    w: FloatArray,
+    periods: FloatArray,
+    *,
+    n_bins: int,
+    dur_bins: list[int],
+    templates: dict[int, FloatArray],
+) -> dict[str, FloatArray]:
+    """Numba CPU matched filter; same output dict as :func:`_matched_filter`."""
+    kernel = _numba_tls_kernel()
+    tmpl = np.concatenate([templates[wd] for wd in dur_bins]).astype(np.float64)
+    widths = np.asarray([len(templates[wd]) for wd in dur_bins], dtype=np.int64)
+    offsets = np.zeros(len(dur_bins), dtype=np.int64)
+    if len(dur_bins) > 1:
+        offsets[1:] = np.cumsum(widths[:-1])
+    sr, depth, duration, t0 = kernel(
+        np.ascontiguousarray(tau), np.ascontiguousarray(yw), np.ascontiguousarray(w),
+        np.ascontiguousarray(periods),
+        np.asarray(dur_bins, dtype=np.int64), offsets, tmpl,
+        int(n_bins), float(_W_EPS),
+    )
+    return {"sr": sr, "depth": depth, "duration": duration, "t0": t0}
+
+
 def tls_power(
     t: FloatArray,
     y: FloatArray,
@@ -404,6 +507,11 @@ def tls_power(
             tau, yw, w, periods_host,
             n_bins=n_bins, dur_bins=dur_bins, templates=templates,
         )
+    elif backend == "numba":
+        res = _tls_numba(
+            tau, yw, w, periods_host,
+            n_bins=n_bins, dur_bins=dur_bins, templates=templates,
+        )
     elif backend == "torch" or backend.startswith("torch:"):
         import torch
 
@@ -438,7 +546,7 @@ def tls_power(
 
 
 class TLSMethod(PeriodogramMethod):
-    """Transit Least Squares — limb-darkened matched filter (numpy CPU, cupy GPU)."""
+    """Transit Least Squares — limb-darkened matched filter (numba/numpy CPU, cupy GPU)."""
 
     name: ClassVar[str] = "TLS"
     objective_sense: ClassVar[Literal["max", "min"]] = "max"
@@ -446,9 +554,10 @@ class TLSMethod(PeriodogramMethod):
     natural_domain: ClassVar[Domain] = Domain.FLUX
     settings_cls: ClassVar[type] = TLSSettings
     cpu_backend: ClassVar[str] = "numpy"
+    fast_cpu_backend: ClassVar[str | None] = "numba"
     gpu_backend: ClassVar[str | None] = "cupy"
     portable_gpu_backend: ClassVar[str | None] = "torch"
-    all_backends: ClassVar[tuple[str, ...]] = ("numpy", "cupy", "torch")
+    all_backends: ClassVar[tuple[str, ...]] = ("numba", "numpy", "cupy", "torch")
 
     def default_grid(self, lc: LightCurve, settings: TLSSettings) -> GridSpec:  # type: ignore[override]
         finite = lc.finite()

From dc75b4f93783b0318bb603edce075e228245b1b2 Mon Sep 17 00:00:00 2001
From: Tharindu Jayasinghe <tharindujay96@gmail.com>
Date: Wed, 1 Jul 2026 14:30:27 -0700
Subject: [PATCH 6/9] perf(cuda): precision-templated RawKernels, CE int
 atomics, BLS device cache

- BLS/PDM/CE/TLS kernel sources are templated on the working precision
  (REAL) and honor settings.precision on the cupy backends: float64 stays
  the default; float32 is an opt-in that runs the FLOP-bound box/matched-
  filter scans at full consumer-GPU rate (Blackwell GeForce: BLS 610 ->
  71 ms, TLS 34 -> 4 ms, ~8.6x). PDM keeps double accumulators (the
  ssd = sq - sum^2/cnt form is cancellation-prone) and CE counts are now
  exact int shared atomics (half the shared memory) with the entropy
  reduction spread across the block instead of one thread.
- TLS caches its concatenated templates in shared memory.
- BLS gains a caller-owned per-light-curve device cache so the segmented
  search uploads tau/yw/ivar once per run (cupy + torch), computes t_min
  on the host (no per-segment device sync), and returns all seven outputs
  in one stacked D2H copy. Multiband keeps one cache per band.
- numba BLS zeroes only the n_bins+1 bin prefix actually used per period.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
---
 src/cuperiod/methods/_bls_core.py           | 219 +++++++++++++-------
 src/cuperiod/methods/bls.py                 |  13 +-
 src/cuperiod/methods/conditional_entropy.py |  83 +++++---
 src/cuperiod/methods/pdm.py                 |  42 ++--
 src/cuperiod/methods/tls.py                 |  76 ++++---
 src/cuperiod/multiband/bls_mb.py            |   4 +-
 6 files changed, 272 insertions(+), 165 deletions(-)

diff --git a/src/cuperiod/methods/_bls_core.py b/src/cuperiod/methods/_bls_core.py
index 5594129..c7c83f1 100644
--- a/src/cuperiod/methods/_bls_core.py
+++ b/src/cuperiod/methods/_bls_core.py
@@ -270,8 +270,14 @@ def _kernel(tau, yw, wv, periods, dur_bins, bin_duration, oversample, width,  #
             n_bins = int(np.ceil(period / bin_duration)) + oversample
             if n_bins > width - 1:
                 n_bins = width - 1
-            my = np.zeros(width)
-            mi = np.zeros(width)
+            # Only bins [0, n_bins] are ever written or read for this period, so
+            # zero just that prefix — width is sized for the *longest* period in
+            # the batch and can be ~2x larger than n_bins for the shortest.
+            my = np.empty(width)
+            mi = np.empty(width)
+            for i in range(n_bins + 1):
+                my[i] = 0.0
+                mi[i] = 0.0
             for j in range(n_points):
                 x = tau[j]
                 w = x - period * np.floor(x / period)
@@ -370,15 +376,15 @@ def _bls_search_numba(
 
 _BLS_CUDA_SRC: Final = r"""
 extern "C" __global__ void bls_block(
-    const double* __restrict__ tau,      // (N) times - t_min
-    const double* __restrict__ yw,       // (N) y * ivar
-    const double* __restrict__ wv,       // (N) ivar
-    const double* __restrict__ periods,  // (P)
+    const REAL* __restrict__ tau,        // (N) times - t_min
+    const REAL* __restrict__ yw,         // (N) y * ivar
+    const REAL* __restrict__ wv,         // (N) ivar
+    const REAL* __restrict__ periods,    // (P)
     const int* __restrict__ dur_bins,    // (D) box widths in bins
     const int n_points, const int n_periods, const int n_dur,
-    const double bin_duration, const int oversample, const int width,
-    const double sum_y, const double sum_ivar, const double t_min,
-    const int obj_flag, const double NEG_INF,
+    const REAL bin_duration, const int oversample, const int width,
+    const REAL sum_y, const REAL sum_ivar, const double t_min,
+    const int obj_flag, const REAL NEG_INF, const REAL IVAR_EPS,
     double* o_power, double* o_depth, double* o_depth_err, double* o_depth_snr,
     double* o_duration, double* o_transit_time, double* o_loglike)
 {
@@ -386,21 +392,21 @@ def _bls_search_numba(
     if (pidx >= n_periods) return;
     const int tid = threadIdx.x;
     const int nth = blockDim.x;
-    const double period = periods[pidx];
+    const REAL period = periods[pidx];
 
     int n_bins = (int)ceil(period / bin_duration) + oversample;
     if (n_bins > width - 1) n_bins = width - 1;
 
-    extern __shared__ double sh[];
-    double* my = sh;          // weighted-y per bin
-    double* mi = sh + width;  // ivar per bin
+    extern __shared__ REAL sh[];
+    REAL* my = sh;          // weighted-y per bin
+    REAL* mi = sh + width;  // ivar per bin
 
     for (int i = tid; i < width; i += nth) { my[i] = 0.0; mi[i] = 0.0; }
     __syncthreads();
 
     for (int j = tid; j < n_points; j += nth) {
-        double x = tau[j];
-        double w = x - period * floor(x / period);
+        REAL x = tau[j];
+        REAL w = x - period * floor(x / period);
         int ind = (int)(w / bin_duration) + 1;
         if (ind < 0) ind = 0;
         if (ind > width - 1) ind = width - 1;
@@ -417,56 +423,56 @@ def _bls_search_numba(
     }
     __syncthreads();
 
-    __shared__ double cof_y[CUDA_BLOCK];
-    __shared__ double cof_i[CUDA_BLOCK];
+    __shared__ REAL cof_y[CUDA_BLOCK];
+    __shared__ REAL cof_i[CUDA_BLOCK];
     {
         int span = n_bins + 1;
         int chunk = (span + nth - 1) / nth;
         int lo = tid * chunk;
         int hi = lo + chunk; if (hi > span) hi = span;
-        double ay = 0.0, ai = 0.0;
+        REAL ay = 0.0, ai = 0.0;
         for (int i = lo; i < hi; ++i) {
             ay += my[i]; my[i] = ay; ai += mi[i]; mi[i] = ai;
         }
         cof_y[tid] = ay; cof_i[tid] = ai;
         __syncthreads();
         if (tid == 0) {
-            double sy = 0.0, si = 0.0;
+            REAL sy = 0.0, si = 0.0;
             for (int k = 0; k < nth; ++k) {
-                double ty = cof_y[k], ti = cof_i[k];
+                REAL ty = cof_y[k], ti = cof_i[k];
                 cof_y[k] = sy; cof_i[k] = si; sy += ty; si += ti;
             }
         }
         __syncthreads();
-        double oy = cof_y[tid], oi = cof_i[tid];
+        REAL oy = cof_y[tid], oi = cof_i[tid];
         for (int i = lo; i < hi; ++i) { my[i] += oy; mi[i] += oi; }
     }
     __syncthreads();
 
-    double loc_obj = NEG_INF;
+    REAL loc_obj = NEG_INF;
     int loc_n = 0, loc_d = 0;
     for (int di = 0; di < n_dur; ++di) {
         const int d = dur_bins[di];
         if (d < 1 || d >= n_bins) continue;
         const int n_max = n_bins - d;
         for (int n = tid; n <= n_max; n += nth) {
-            double y_in = my[n + d] - my[n];
-            double iv_in = mi[n + d] - mi[n];
-            double iv_out = sum_ivar - iv_in;
-            if (iv_in < 2.2204460492503131e-16) continue;
-            if (iv_out < 2.2204460492503131e-16) continue;
-            double yin = y_in / iv_in;
-            double yout = (sum_y - y_in) / iv_out;
+            REAL y_in = my[n + d] - my[n];
+            REAL iv_in = mi[n + d] - mi[n];
+            REAL iv_out = sum_ivar - iv_in;
+            if (iv_in < IVAR_EPS) continue;
+            if (iv_out < IVAR_EPS) continue;
+            REAL yin = y_in / iv_in;
+            REAL yout = (sum_y - y_in) / iv_out;
             if (yout < yin) continue;
-            double depth = yout - yin;
-            double obj = (obj_flag == 0)
-                ? depth / sqrt(1.0 / iv_in + 1.0 / iv_out)
-                : 0.5 * iv_in * depth * depth;
+            REAL depth = yout - yin;
+            REAL obj = (obj_flag == 0)
+                ? depth / sqrt((REAL)1.0 / iv_in + (REAL)1.0 / iv_out)
+                : (REAL)0.5 * iv_in * depth * depth;
             if (obj > loc_obj) { loc_obj = obj; loc_n = n; loc_d = d; }
         }
     }
 
-    __shared__ double r_obj[CUDA_BLOCK];
+    __shared__ REAL r_obj[CUDA_BLOCK];
     __shared__ int r_n[CUDA_BLOCK];
     __shared__ int r_d[CUDA_BLOCK];
     r_obj[tid] = loc_obj; r_n[tid] = loc_n; r_d[tid] = loc_d;
@@ -481,7 +487,7 @@ def _bls_search_numba(
     }
 
     if (tid == 0) {
-        double best = r_obj[0];
+        REAL best = r_obj[0];
         if (best <= NEG_INF) {
             o_power[pidx] = NEG_INF; o_depth[pidx] = 0.0; o_depth_err[pidx] = 0.0;
             o_depth_snr[pidx] = 0.0; o_duration[pidx] = 0.0;
@@ -489,40 +495,42 @@ def _bls_search_numba(
             return;
         }
         int bn = r_n[0], bd = r_d[0];
-        double y_in = my[bn + bd] - my[bn];
-        double iv_in = mi[bn + bd] - mi[bn];
-        double iv_out = sum_ivar - iv_in;
+        double y_in = (double)my[bn + bd] - (double)my[bn];
+        double iv_in = (double)mi[bn + bd] - (double)mi[bn];
+        double iv_out = (double)sum_ivar - iv_in;
         double yin = y_in / iv_in;
-        double yout = (sum_y - y_in) / iv_out;
+        double yout = ((double)sum_y - y_in) / iv_out;
         double depth = yout - yin;
         double derr = sqrt(1.0 / iv_in + 1.0 / iv_out);
         double dsnr = depth / derr;
         double loglike = 0.5 * iv_in * depth * depth;
-        double dur = (double)bd * bin_duration;
+        double dur = (double)bd * (double)bin_duration;
         o_power[pidx] = (obj_flag == 0) ? dsnr : loglike;
         o_depth[pidx] = depth;
         o_depth_err[pidx] = derr;
         o_depth_snr[pidx] = dsnr;
         o_duration[pidx] = dur;
         o_transit_time[pidx] =
-            fmod((double)bn * bin_duration + 0.5 * dur, period) + t_min;
+            fmod((double)bn * (double)bin_duration + 0.5 * dur, (double)period) + t_min;
         o_loglike[pidx] = loglike;
     }
 }
 """
 
-_cuda_kernel_cache: dict[int, Any] = {}
+_cuda_kernel_cache: dict[tuple[int, str], Any] = {}
 
 
-def _cuda_kernel(block: int) -> Any:
-    """Compile (once) and cache the BLS RawKernel for a given block size."""
-    kernel = _cuda_kernel_cache.get(block)
+def _cuda_kernel(block: int, real: str) -> Any:
+    """Compile (once) and cache the BLS RawKernel for a block size and precision."""
+    kernel = _cuda_kernel_cache.get((block, real))
     if kernel is None:
         import cupy
 
-        src = _BLS_CUDA_SRC.replace("CUDA_BLOCK", str(block))
+        src = _BLS_CUDA_SRC.replace("CUDA_BLOCK", str(block)).replace(
+            "REAL", "float" if real == "float32" else "double"
+        )
         kernel = cupy.RawKernel(src, "bls_block")
-        _cuda_kernel_cache[block] = kernel
+        _cuda_kernel_cache[(block, real)] = kernel
     return kernel
 
 
@@ -537,20 +545,43 @@ def _bls_search_cuda(
     oversample: int,
     width: int,
     obj_flag: int,
+    precision: str = "float64",
+    device_cache: dict[str, Any] | None = None,
     block: int = CUDA_BLOCK,
 ) -> dict[str, Any]:
-    """One-block-per-period CUDA box search; returns cupy arrays per output."""
+    """One-block-per-period CUDA box search; returns host float64 arrays per output.
+
+    ``device_cache`` (caller-owned, one per light curve) keeps the uploaded
+    ``tau``/``yw``/``ivar`` device arrays and their float64 sums across calls, so the
+    segmented search uploads each light curve once instead of once per segment. The
+    time origin is subtracted on the host in float64 before any (possibly float32)
+    device cast; the kernel adds ``t_min`` back into ``transit_time`` in double.
+    """
     import cupy as cp
 
-    t_d = cp.asarray(np.ascontiguousarray(t, dtype=np.float64))
-    ivar_d = cp.asarray(np.ascontiguousarray(ivar, dtype=np.float64))
-    y_d = cp.asarray(np.ascontiguousarray(y, dtype=np.float64))
-    t_min = float(t_d.min())
-    tau = t_d - t_min
-    yw = y_d * ivar_d
-    sum_y = float(yw.sum())
-    sum_ivar = float(ivar_d.sum())
-    periods_d = cp.asarray(np.ascontiguousarray(periods, dtype=np.float64))
+    rdtype = np.float32 if precision == "float32" else np.float64
+    data: dict[str, Any] | None = None
+    if device_cache is not None:
+        data = device_cache.get("cuda")
+        if data is not None and (data["n"] != t.shape[0] or data["dtype"] != rdtype):
+            data = None
+    if data is None:
+        t_min = float(t.min())
+        yw_h = y * ivar
+        data = {
+            "n": int(t.shape[0]),
+            "dtype": rdtype,
+            "t_min": t_min,
+            "tau": cp.asarray(np.ascontiguousarray(t - t_min, dtype=rdtype)),
+            "yw": cp.asarray(np.ascontiguousarray(yw_h, dtype=rdtype)),
+            "ivar": cp.asarray(np.ascontiguousarray(ivar, dtype=rdtype)),
+            "sum_y": float(yw_h.sum()),
+            "sum_ivar": float(ivar.sum()),
+        }
+        if device_cache is not None:
+            device_cache["cuda"] = data
+
+    periods_d = cp.asarray(np.ascontiguousarray(periods, dtype=rdtype))
     n_periods = int(periods_d.shape[0])
     dur_d = cp.asarray(np.asarray(dur_bins, dtype=np.int32))
 
@@ -563,15 +594,17 @@ def _bls_search_cuda(
         "transit_time",
         "log_likelihood",
     )
-    out = {name: cp.empty(n_periods, dtype=cp.float64) for name in names}
     if n_periods == 0 or len(dur_bins) == 0:
-        out["power"][...] = -np.inf
-        return out
+        out_empty = {name: np.zeros(n_periods, dtype=np.float64) for name in names}
+        out_empty["power"][...] = -np.inf
+        return out_empty
+    out = {name: cp.empty(n_periods, dtype=cp.float64) for name in names}
 
     from cuperiod.core.backend import ensure_shared_memory
 
-    smem = 2 * width * 8
-    kernel = _cuda_kernel(block)
+    real_size = 4 if precision == "float32" else 8
+    smem = 2 * width * real_size
+    kernel = _cuda_kernel(block, precision)
     ensure_shared_memory(
         kernel, smem, method="BLS", hint="the period range (max_period_days)"
     )
@@ -579,22 +612,23 @@ def _bls_search_cuda(
         (n_periods,),
         (block,),
         (
-            tau,
-            yw,
-            ivar_d,
+            data["tau"],
+            data["yw"],
+            data["ivar"],
             periods_d,
             dur_d,
-            np.int32(t_d.shape[0]),
+            np.int32(data["n"]),
             np.int32(n_periods),
             np.int32(dur_d.shape[0]),
-            np.float64(bin_duration),
+            rdtype(bin_duration),
             np.int32(oversample),
             np.int32(width),
-            np.float64(sum_y),
-            np.float64(sum_ivar),
-            np.float64(t_min),
+            rdtype(data["sum_y"]),
+            rdtype(data["sum_ivar"]),
+            np.float64(data["t_min"]),
             np.int32(obj_flag),
-            np.float64(-np.inf),
+            rdtype(-np.inf),
+            rdtype(_IVAR_EPS),
             out["power"],
             out["depth"],
             out["depth_err"],
@@ -605,7 +639,9 @@ def _bls_search_cuda(
         ),
         shared_mem=smem,
     )
-    return out
+    # One stacked D2H copy instead of seven independent transfers.
+    host = cp.asnumpy(cp.stack([out[name] for name in names]))
+    return {name: host[i] for i, name in enumerate(names)}
 
 
 def bls_power(
@@ -620,6 +656,7 @@ def bls_power(
     backend: str = "numpy",
     batch: int = DEFAULT_BATCH,
     precision: str = "auto",
+    device_cache: dict[str, Any] | None = None,
 ) -> BLSPower:
     """BLS box search over ``periods`` via numpy/torch (portable), cupy, or numba.
 
@@ -646,7 +683,12 @@ def bls_power(
     batch : int, default 2048
         Trial periods per vectorized batch (numpy/torch backends).
     precision : {"auto", "float64", "float32"}, default "auto"
-        Device-side compute precision for the torch backend (float64 except on MPS).
+        Device-side compute precision for the torch and cupy backends (float64
+        except on MPS; ``"float32"`` is an opt-in speedup on consumer GPUs).
+    device_cache : dict or None
+        Caller-owned scratch dict reused across calls that share the *same*
+        ``t``/``y``/``dy`` (the segmented search): device uploads happen once per
+        light curve instead of once per call. Pass a fresh ``{}`` per light curve.
 
     Returns
     -------
@@ -688,6 +730,8 @@ def bls_power(
             t_host, y_host, ivar_host, periods_host,
             bin_duration=bin_duration, dur_bins=dur_bins,
             oversample=oversample, width=width, obj_flag=obj_flag,
+            precision=resolve_precision(precision, "cuda"),
+            device_cache=device_cache,
         )
     elif backend == "numba":
         out = _bls_search_numba(
@@ -709,13 +753,32 @@ def bls_power(
         # precision on Apple MPS — so a raw cast would silently corrupt the phase fold.
         # The small t-min-relative tau survives float32; ``transit_time`` is shifted
         # back to absolute on the host (float64) at the return. Mirrors ``gls._prep``.
-        t_ref = float(np.min(t_host))
-        t_d = to_device_array(t_host - t_ref, device=device, dtype=tdtype)
-        y_d = to_device_array(y_host, device=device, dtype=tdtype)
-        ivar_d = to_device_array(ivar_host, device=device, dtype=tdtype)
+        data: dict[str, Any] | None = None
+        if device_cache is not None:
+            data = device_cache.get("torch")
+            if data is not None and (
+                data["n"] != t_host.shape[0]
+                or data["dtype"] != tdtype
+                or data["device"] != device
+            ):
+                data = None
+        if data is None:
+            t_ref_val = float(np.min(t_host))
+            data = {
+                "n": int(t_host.shape[0]),
+                "dtype": tdtype,
+                "device": device,
+                "t_ref": t_ref_val,
+                "t": to_device_array(t_host - t_ref_val, device=device, dtype=tdtype),
+                "y": to_device_array(y_host, device=device, dtype=tdtype),
+                "ivar": to_device_array(ivar_host, device=device, dtype=tdtype),
+            }
+            if device_cache is not None:
+                device_cache["torch"] = data
+        t_ref = data["t_ref"]
         periods_d = to_device_array(periods_host, device=device, dtype=tdtype)
         out = _bls_search(
-            array_namespace(periods_d), t_d, y_d, ivar_d, periods_d,
+            array_namespace(periods_d), data["t"], data["y"], data["ivar"], periods_d,
             bin_duration=bin_duration, dur_bins=dur_bins,
             oversample=oversample, width=width, obj_flag=obj_flag, batch=batch,
         )
diff --git a/src/cuperiod/methods/bls.py b/src/cuperiod/methods/bls.py
index 54ab348..18509c6 100644
--- a/src/cuperiod/methods/bls.py
+++ b/src/cuperiod/methods/bls.py
@@ -78,8 +78,13 @@ def _segment_power(
     periods: FloatArray,
     durations: FloatArray,
     settings: BLSSettings,
+    device_cache: dict[str, Any] | None = None,
 ) -> dict[str, FloatArray]:
-    """One segment's per-period box maxima via the configured backend."""
+    """One segment's per-period box maxima via the configured backend.
+
+    ``device_cache`` is a caller-owned dict shared by the segments of one light
+    curve, so the GPU backends upload the curve once per run rather than per segment.
+    """
     if backend == "astropy":
         from astropy.timeseries import BoxLeastSquares
 
@@ -104,6 +109,7 @@ def _segment_power(
         backend=backend,
         batch=settings.batch_periods,
         precision=settings.precision,
+        device_cache=device_cache,
     )
     return {name: getattr(power, name) for name in _SEGMENT_FIELDS}
 
@@ -214,8 +220,11 @@ def power(  # type: ignore[override]
             )
 
         chunks: dict[str, list[FloatArray]] = {name: [] for name in _SEGMENT_FIELDS}
+        device_cache: dict[str, Any] = {}
         for periods, durations in segments:
-            seg = _segment_power(backend, jd, flux, err, periods, durations, settings)
+            seg = _segment_power(
+                backend, jd, flux, err, periods, durations, settings, device_cache
+            )
             for name in _SEGMENT_FIELDS:
                 chunks[name].append(seg[name])
         return _assemble(chunks, n, finite.baseline, backend, finite.meta)
diff --git a/src/cuperiod/methods/conditional_entropy.py b/src/cuperiod/methods/conditional_entropy.py
index 573612d..c17612d 100644
--- a/src/cuperiod/methods/conditional_entropy.py
+++ b/src/cuperiod/methods/conditional_entropy.py
@@ -94,12 +94,16 @@ def _entropy_batch(
 CUDA_BLOCK: Final = 128
 
 #: One-block-per-period CE kernel: a block folds its points into a shared-memory 2-D
-#: phase-magnitude histogram (magnitude bins precomputed on the host), then one thread
-#: reduces it to the Shannon conditional entropy H(m | phase).
+#: phase-magnitude histogram (magnitude bins precomputed on the host), then reduces it
+#: to the Shannon conditional entropy H(m | phase). Counts are **integer** shared
+#: atomics — exact, faster than double atomics, and half the shared memory — and the
+#: entropy reduction is spread across the block (phase rows striped over threads, then
+#: a tree sum). ``REAL`` is the working precision for the fold (float64 by default;
+#: float32 opt-in via ``precision``); counts and entropy stay exact/double either way.
 _CE_CUDA_SRC: Final = r"""
 extern "C" __global__ void ce_block(
-    const double* __restrict__ tau, const int* __restrict__ mag_bin,
-    const double* __restrict__ periods,
+    const REAL* __restrict__ tau, const int* __restrict__ mag_bin,
+    const REAL* __restrict__ periods,
     const int n_points, const int n_periods, const int n_phase, const int n_mag,
     double* o_entropy)
 {
@@ -107,52 +111,60 @@ def _entropy_batch(
     if (pidx >= n_periods) return;
     const int tid = threadIdx.x;
     const int nth = blockDim.x;
-    const double period = periods[pidx];
+    const REAL period = periods[pidx];
     const int M = n_phase * n_mag;
 
-    extern __shared__ double cnt[];  // (M) point count per (phase, mag) cell
-    for (int i = tid; i < M; i += nth) cnt[i] = 0.0;
+    extern __shared__ int cnt[];  // (M) point count per (phase, mag) cell
+    for (int i = tid; i < M; i += nth) cnt[i] = 0;
     __syncthreads();
 
     for (int j = tid; j < n_points; j += nth) {
-        double x = tau[j];
-        double ph = (x - period * floor(x / period)) / period;  // mod(tau/period, 1)
-        int pb = (int)(ph * n_phase);
+        REAL x = tau[j];
+        REAL ph = (x - period * floor(x / period)) / period;  // mod(tau/period, 1)
+        int pb = (int)(ph * (REAL)n_phase);
         if (pb >= n_phase) pb = n_phase - 1;
         if (pb < 0) pb = 0;
-        atomicAdd(&cnt[pb * n_mag + mag_bin[j]], 1.0);
+        atomicAdd(&cnt[pb * n_mag + mag_bin[j]], 1);
     }
     __syncthreads();
 
-    if (tid == 0) {
-        double h = 0.0;
-        for (int p = 0; p < n_phase; ++p) {
-            double ci = 0.0;
-            for (int m = 0; m < n_mag; ++m) ci += cnt[p * n_mag + m];
-            if (ci > 0.0) {
-                double lci = log(ci);
-                for (int m = 0; m < n_mag; ++m) {
-                    double c = cnt[p * n_mag + m];
-                    if (c > 0.0) h += c * (lci - log(c));
-                }
+    double h = 0.0;
+    for (int p = tid; p < n_phase; p += nth) {
+        double ci = 0.0;
+        for (int m = 0; m < n_mag; ++m) ci += (double)cnt[p * n_mag + m];
+        if (ci > 0.0) {
+            double lci = log(ci);
+            for (int m = 0; m < n_mag; ++m) {
+                double c = (double)cnt[p * n_mag + m];
+                if (c > 0.0) h += c * (lci - log(c));
             }
         }
-        o_entropy[pidx] = h / (double)n_points;
     }
+    __shared__ double r_h[CUDA_BLOCK];
+    r_h[tid] = h;
+    __syncthreads();
+    for (int s = blockDim.x >> 1; s > 0; s >>= 1) {
+        if (tid < s) r_h[tid] += r_h[tid + s];
+        __syncthreads();
+    }
+    if (tid == 0) o_entropy[pidx] = r_h[0] / (double)n_points;
 }
 """
 
-_ce_kernel_cache: dict[int, Any] = {}
+_ce_kernel_cache: dict[tuple[int, str], Any] = {}
 
 
-def _ce_kernel(block: int) -> Any:
-    """Compile (once) and cache the CE RawKernel for a given block size."""
-    kernel = _ce_kernel_cache.get(block)
+def _ce_kernel(block: int, real: str) -> Any:
+    """Compile (once) and cache the CE RawKernel for a block size and precision."""
+    kernel = _ce_kernel_cache.get((block, real))
     if kernel is None:
         import cupy
 
-        kernel = cupy.RawKernel(_CE_CUDA_SRC, "ce_block")
-        _ce_kernel_cache[block] = kernel
+        src = _CE_CUDA_SRC.replace("CUDA_BLOCK", str(block)).replace(
+            "REAL", "float" if real == "float32" else "double"
+        )
+        kernel = cupy.RawKernel(src, "ce_block")
+        _ce_kernel_cache[(block, real)] = kernel
     return kernel
 
 
@@ -163,22 +175,24 @@ def _ce_cuda(
     *,
     n_phase: int,
     n_mag: int,
+    precision: str = "float64",
     block: int = CUDA_BLOCK,
 ) -> FloatArray:
     """One-block-per-period CUDA conditional entropy; returns a host float64 array."""
     import cupy as cp
 
-    tau_d = cp.asarray(np.ascontiguousarray(tau, dtype=np.float64))
+    rdtype = np.float32 if precision == "float32" else np.float64
+    tau_d = cp.asarray(np.ascontiguousarray(tau, dtype=rdtype))
     mag_d = cp.asarray(np.ascontiguousarray(mag_bin, dtype=np.int32))
-    per_d = cp.asarray(np.ascontiguousarray(periods, dtype=np.float64))
+    per_d = cp.asarray(np.ascontiguousarray(periods, dtype=rdtype))
     n_periods = int(per_d.size)
     if n_periods == 0:
         return np.zeros(0, dtype=np.float64)
     out = cp.empty(n_periods, dtype=cp.float64)
     from cuperiod.core.backend import ensure_shared_memory
 
-    smem = n_phase * n_mag * 8
-    kernel = _ce_kernel(block)
+    smem = n_phase * n_mag * 4
+    kernel = _ce_kernel(block, precision)
     ensure_shared_memory(kernel, smem, method="CE", hint="n_phase_bins / n_mag_bins")
     kernel(
         (n_periods,),
@@ -293,7 +307,8 @@ def conditional_entropy(
     if backend == "cupy":
         ensure_cuda_dll_path()
         return _ce_cuda(
-            tau, mag_bin, periods_host, n_phase=n_phase_bins, n_mag=n_mag_bins
+            tau, mag_bin, periods_host, n_phase=n_phase_bins, n_mag=n_mag_bins,
+            precision=resolve_precision(precision, "cuda"),
         )
     if backend == "numba":
         kernel = _numba_ce_kernel()
diff --git a/src/cuperiod/methods/pdm.py b/src/cuperiod/methods/pdm.py
index 1cad0ed..c5c4fd6 100644
--- a/src/cuperiod/methods/pdm.py
+++ b/src/cuperiod/methods/pdm.py
@@ -137,8 +137,8 @@ def _theta_batch(
 #: dominant atomic work per point from ``3*n_covers`` to 3.
 _PDM_CUDA_SRC: Final = r"""
 extern "C" __global__ void pdm_block(
-    const double* __restrict__ tau, const double* __restrict__ y,
-    const double* __restrict__ periods,
+    const REAL* __restrict__ tau, const REAL* __restrict__ y,
+    const REAL* __restrict__ periods,
     const int n_points, const int n_periods, const int n_bins, const int n_covers,
     const double sigma2, double* o_theta)
 {
@@ -146,9 +146,11 @@ def _theta_batch(
     if (pidx >= n_periods) return;
     const int tid = threadIdx.x;
     const int nth = blockDim.x;
-    const double period = periods[pidx];
+    const REAL period = periods[pidx];
     const int M = n_bins * n_covers;
 
+    // Accumulators stay double even at REAL=float: the within-bin sum of squared
+    // deviations (sq - sum^2/cnt) is cancellation-prone in float32.
     extern __shared__ double sh[];
     double* s_sum = sh;          // (M) sum of y per fine bin
     double* s_sq = sh + M;       // (M) sum of y^2 per fine bin
@@ -159,10 +161,10 @@ def _theta_batch(
     __syncthreads();
 
     for (int j = tid; j < n_points; j += nth) {
-        double x = tau[j];
-        double ph = (x - period * floor(x / period)) / period;  // mod(tau/period, 1)
-        double yj = y[j];
-        int f = (int)(ph * (double)M);
+        REAL x = tau[j];
+        REAL ph = (x - period * floor(x / period)) / period;  // mod(tau/period, 1)
+        double yj = (double)y[j];
+        int f = (int)(ph * (REAL)M);
         if (f >= M) f = M - 1;
         if (f < 0) f = 0;
         atomicAdd(&s_sum[f], yj);
@@ -197,17 +199,18 @@ def _theta_batch(
 }
 """
 
-_pdm_kernel_cache: dict[int, Any] = {}
+_pdm_kernel_cache: dict[tuple[int, str], Any] = {}
 
 
-def _pdm_kernel(block: int) -> Any:
-    """Compile (once) and cache the PDM RawKernel for a given block size."""
-    kernel = _pdm_kernel_cache.get(block)
+def _pdm_kernel(block: int, real: str) -> Any:
+    """Compile (once) and cache the PDM RawKernel for a block size and precision."""
+    kernel = _pdm_kernel_cache.get((block, real))
     if kernel is None:
         import cupy
 
-        kernel = cupy.RawKernel(_PDM_CUDA_SRC, "pdm_block")
-        _pdm_kernel_cache[block] = kernel
+        src = _PDM_CUDA_SRC.replace("REAL", "float" if real == "float32" else "double")
+        kernel = cupy.RawKernel(src, "pdm_block")
+        _pdm_kernel_cache[(block, real)] = kernel
     return kernel
 
 
@@ -219,14 +222,16 @@ def _pdm_cuda(
     n_bins: int,
     n_covers: int,
     sigma2: float,
+    precision: str = "float64",
     block: int = CUDA_BLOCK,
 ) -> FloatArray:
     """One-block-per-period CUDA PDM Theta; returns a host float64 array."""
     import cupy as cp
 
-    tau_d = cp.asarray(np.ascontiguousarray(tau, dtype=np.float64))
-    y_d = cp.asarray(np.ascontiguousarray(y, dtype=np.float64))
-    per_d = cp.asarray(np.ascontiguousarray(periods, dtype=np.float64))
+    rdtype = np.float32 if precision == "float32" else np.float64
+    tau_d = cp.asarray(np.ascontiguousarray(tau, dtype=rdtype))
+    y_d = cp.asarray(np.ascontiguousarray(y, dtype=rdtype))
+    per_d = cp.asarray(np.ascontiguousarray(periods, dtype=rdtype))
     n_periods = int(per_d.size)
     if n_periods == 0:
         return np.zeros(0, dtype=np.float64)
@@ -234,7 +239,7 @@ def _pdm_cuda(
     from cuperiod.core.backend import ensure_shared_memory
 
     smem = 3 * n_bins * n_covers * 8
-    kernel = _pdm_kernel(block)
+    kernel = _pdm_kernel(block, precision)
     ensure_shared_memory(kernel, smem, method="PDM", hint="n_bins / n_covers")
     kernel(
         (n_periods,),
@@ -369,7 +374,8 @@ def pdm_theta(
     if backend == "cupy":
         ensure_cuda_dll_path()
         return _pdm_cuda(
-            tau, y, periods_host, n_bins=n_bins, n_covers=n_covers, sigma2=sigma2
+            tau, y, periods_host, n_bins=n_bins, n_covers=n_covers, sigma2=sigma2,
+            precision=resolve_precision(precision, "cuda"),
         )
     if backend == "numba":
         kernel = _numba_pdm_kernel()
diff --git a/src/cuperiod/methods/tls.py b/src/cuperiod/methods/tls.py
index 7456344..d41cd33 100644
--- a/src/cuperiod/methods/tls.py
+++ b/src/cuperiod/methods/tls.py
@@ -198,29 +198,31 @@ def _matched_filter(
 #: block-reduce to the best signal residue. No (period, point) global intermediates.
 _TLS_CUDA_SRC: Final = r"""
 extern "C" __global__ void tls_block(
-    const double* __restrict__ tau, const double* __restrict__ yw,
-    const double* __restrict__ wv, const double* __restrict__ periods,
+    const REAL* __restrict__ tau, const REAL* __restrict__ yw,
+    const REAL* __restrict__ wv, const REAL* __restrict__ periods,
     const int* __restrict__ dur_bins, const int* __restrict__ tmpl_off,
-    const double* __restrict__ tmpl, const int n_dur,
-    const int n_points, const int n_periods, const int n_bins, const double W_EPS,
+    const REAL* __restrict__ tmpl, const int n_dur, const int tmpl_len,
+    const int n_points, const int n_periods, const int n_bins, const REAL W_EPS,
     double* o_sr, double* o_depth, double* o_duration, double* o_t0)
 {
     const int pidx = blockIdx.x;
     if (pidx >= n_periods) return;
     const int tid = threadIdx.x;
     const int nth = blockDim.x;
-    const double period = periods[pidx];
+    const REAL period = periods[pidx];
 
-    extern __shared__ double sh[];
-    double* a = sh;             // (n_bins) sum w*y' per bin
-    double* b = sh + n_bins;    // (n_bins) sum w per bin
+    extern __shared__ REAL sh[];
+    REAL* a = sh;                        // (n_bins) sum w*y' per bin
+    REAL* b = sh + n_bins;               // (n_bins) sum w per bin
+    REAL* s_tmpl = sh + 2 * n_bins;      // (tmpl_len) all templates, concatenated
     for (int i = tid; i < n_bins; i += nth) { a[i] = 0.0; b[i] = 0.0; }
+    for (int i = tid; i < tmpl_len; i += nth) s_tmpl[i] = tmpl[i];
     __syncthreads();
 
     for (int j = tid; j < n_points; j += nth) {
-        double x = tau[j];
-        double ph = (x - period * floor(x / period)) / period;
-        int bb = (int)(ph * n_bins);
+        REAL x = tau[j];
+        REAL ph = (x - period * floor(x / period)) / period;
+        int bb = (int)(ph * (REAL)n_bins);
         if (bb >= n_bins) bb = n_bins - 1;
         if (bb < 0) bb = 0;
         atomicAdd(&a[bb], yw[j]);
@@ -228,23 +230,23 @@ def _matched_filter(
     }
     __syncthreads();
 
-    double loc_sr = 0.0, loc_depth = 0.0;
+    REAL loc_sr = 0.0, loc_depth = 0.0;
     int loc_start = 0, loc_width = 0;
     for (int idx = tid; idx < n_dur * n_bins; idx += nth) {
         const int di = idx / n_bins;
         const int start = idx % n_bins;
         const int width = dur_bins[di];
         const int off = tmpl_off[di];
-        double num = 0.0, den = 0.0;
+        REAL num = 0.0, den = 0.0;
         for (int k = 0; k < width; ++k) {
             int bb = start + k;
             if (bb >= n_bins) bb -= n_bins;  // circular wrap
-            double gk = tmpl[off + k];
+            REAL gk = s_tmpl[off + k];
             num += gk * a[bb];
             den += gk * gk * b[bb];
         }
         if (den > W_EPS && num < 0.0) {
-            double sr = num * num / den;
+            REAL sr = num * num / den;
             if (sr > loc_sr) {
                 loc_sr = sr; loc_start = start; loc_width = width;
                 loc_depth = -num / den;
@@ -252,8 +254,8 @@ def _matched_filter(
         }
     }
 
-    __shared__ double r_sr[CUDA_BLOCK];
-    __shared__ double r_depth[CUDA_BLOCK];
+    __shared__ REAL r_sr[CUDA_BLOCK];
+    __shared__ REAL r_depth[CUDA_BLOCK];
     __shared__ int r_start[CUDA_BLOCK];
     __shared__ int r_width[CUDA_BLOCK];
     r_sr[tid] = loc_sr; r_depth[tid] = loc_depth;
@@ -278,18 +280,20 @@ def _matched_filter(
 }
 """
 
-_tls_kernel_cache: dict[int, Any] = {}
+_tls_kernel_cache: dict[tuple[int, str], Any] = {}
 
 
-def _tls_kernel(block: int) -> Any:
-    """Compile (once) and cache the TLS RawKernel for a given block size."""
-    kernel = _tls_kernel_cache.get(block)
+def _tls_kernel(block: int, real: str) -> Any:
+    """Compile (once) and cache the TLS RawKernel for a block size and precision."""
+    kernel = _tls_kernel_cache.get((block, real))
     if kernel is None:
         import cupy
 
-        src = _TLS_CUDA_SRC.replace("CUDA_BLOCK", str(block))
+        src = _TLS_CUDA_SRC.replace("CUDA_BLOCK", str(block)).replace(
+            "REAL", "float" if real == "float32" else "double"
+        )
         kernel = cupy.RawKernel(src, "tls_block")
-        _tls_kernel_cache[block] = kernel
+        _tls_kernel_cache[(block, real)] = kernel
     return kernel
 
 
@@ -302,20 +306,23 @@ def _tls_cuda(
     n_bins: int,
     dur_bins: list[int],
     templates: dict[int, FloatArray],
+    precision: str = "float64",
     block: int = CUDA_BLOCK,
 ) -> dict[str, FloatArray]:
     """One-block-per-period CUDA matched filter; returns host arrays per output."""
     import cupy as cp
 
-    flat = np.concatenate([templates[wd] for wd in dur_bins]).astype(np.float64)
+    rdtype = np.float32 if precision == "float32" else np.float64
+    real_size = 4 if precision == "float32" else 8
+    flat = np.concatenate([templates[wd] for wd in dur_bins]).astype(rdtype)
     widths = [len(templates[wd]) for wd in dur_bins]
     offsets = np.zeros(len(dur_bins), dtype=np.int32)
     if len(widths) > 1:
         offsets[1:] = np.cumsum(widths[:-1])
-    tau_d = cp.asarray(np.ascontiguousarray(tau, dtype=np.float64))
-    yw_d = cp.asarray(np.ascontiguousarray(yw, dtype=np.float64))
-    w_d = cp.asarray(np.ascontiguousarray(w, dtype=np.float64))
-    per_d = cp.asarray(np.ascontiguousarray(periods, dtype=np.float64))
+    tau_d = cp.asarray(np.ascontiguousarray(tau, dtype=rdtype))
+    yw_d = cp.asarray(np.ascontiguousarray(yw, dtype=rdtype))
+    w_d = cp.asarray(np.ascontiguousarray(w, dtype=rdtype))
+    per_d = cp.asarray(np.ascontiguousarray(periods, dtype=rdtype))
     dur_d = cp.asarray(np.asarray(dur_bins, dtype=np.int32))
     off_d = cp.asarray(offsets)
     tmpl_d = cp.asarray(flat)
@@ -326,16 +333,20 @@ def _tls_cuda(
         return {name: np.zeros(0, dtype=np.float64) for name in names}
     from cuperiod.core.backend import ensure_shared_memory
 
-    smem = 2 * n_bins * 8
-    kernel = _tls_kernel(block)
+    # bins (a, b) plus the concatenated templates all live in shared memory.
+    smem = (2 * n_bins + int(flat.size)) * real_size
+    kernel = _tls_kernel(block, precision)
     ensure_shared_memory(kernel, smem, method="TLS", hint="n_phase_bins")
     kernel(
         (n_periods,),
         (block,),
         (
             tau_d, yw_d, w_d, per_d, dur_d, off_d, tmpl_d,
-            np.int32(len(dur_bins)), np.int32(tau_d.size), np.int32(n_periods),
-            np.int32(n_bins), np.float64(_W_EPS),
+            np.int32(len(dur_bins)), np.int32(flat.size),
+            np.int32(tau_d.size), np.int32(n_periods),
+            np.int32(n_bins),
+            # _W_EPS (1e-300) underflows to 0 in float32; use a float32-scale floor.
+            rdtype(1e-30) if precision == "float32" else np.float64(_W_EPS),
             out["sr"], out["depth"], out["duration"], out["t0"],
         ),
         shared_mem=smem,
@@ -506,6 +517,7 @@ def tls_power(
         res = _tls_cuda(
             tau, yw, w, periods_host,
             n_bins=n_bins, dur_bins=dur_bins, templates=templates,
+            precision=resolve_precision(settings.precision, "cuda"),
         )
     elif backend == "numba":
         res = _tls_numba(
diff --git a/src/cuperiod/multiband/bls_mb.py b/src/cuperiod/multiband/bls_mb.py
index 13e785e..24ed3f6 100644
--- a/src/cuperiod/multiband/bls_mb.py
+++ b/src/cuperiod/multiband/bls_mb.py
@@ -81,9 +81,10 @@ def bls_multiband_power(
     dsnr_parts: list[FloatArray] = []
     dur_parts: list[FloatArray] = []
     t0_parts: list[FloatArray] = []
+    band_caches: list[dict[str, object]] = [{} for _ in bands]
     for periods, durations in segments:
         per_band = []
-        for lc in bands:
+        for lc, device_cache in zip(bands, band_caches):
             err = lc.error if lc.error is not None else np.ones_like(lc.value)
             seg = _segment_power(
                 backend,  # type: ignore[arg-type]
@@ -93,6 +94,7 @@ def bls_multiband_power(
                 periods,
                 durations,
                 settings,
+                device_cache,  # type: ignore[arg-type]
             )
             per_band.append(seg)
         dsnr = np.vstack([np.clip(s["depth_snr"], 0.0, None) for s in per_band])

From 27a1a09fce4c3322db93a7bc4debeb36156be48e Mon Sep 17 00:00:00 2001
From: Tharindu Jayasinghe <tharindujay96@gmail.com>
Date: Wed, 1 Jul 2026 14:34:37 -0700
Subject: [PATCH 7/9] perf(string-length): in-block CUDA sort kernel +
 torch.sort fast path

- New one-block-per-period RawKernel: folds, bitonic-sorts the (phase,
  original-index) pairs in shared memory (index tie-break = the same
  stable order as the array-API paths), and reduces the string length
  with no (P, N) intermediates. Curves beyond the shared-memory
  capacity (~8k points on Blackwell) fall back to the vectorized path.
  3k points x 50k periods: 123 ms vs 9.3 s numpy, parity 1e-13 with
  forced phase ties.
- torch path uses torch.sort(stable=True) to fuse the argsort +
  row-gather pair (57 ms on torch:cuda).
- ensure_shared_memory now opts in whenever dynamic+static shared
  exceeds the 48 KB default; dynamic == 48 KB alongside static arrays
  previously failed with CUDA_ERROR_INVALID_VALUE.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
---
 src/cuperiod/core/backend.py          |   5 +-
 src/cuperiod/methods/string_length.py | 161 +++++++++++++++++++++++++-
 2 files changed, 162 insertions(+), 4 deletions(-)

diff --git a/src/cuperiod/core/backend.py b/src/cuperiod/core/backend.py
index 5305901..423a796 100644
--- a/src/cuperiod/core/backend.py
+++ b/src/cuperiod/core/backend.py
@@ -193,7 +193,10 @@ def ensure_shared_memory(
             f"the device allows at most {optin // 1024} KB; reduce {hint}, or use "
             "backend='cpu'."
         )
-    if int(dynamic_bytes) > _DEFAULT_SHARED_MEM:
+    # The opt-in is needed whenever the *total* (dynamic + static) exceeds the
+    # default cap — a launch with 48 KB dynamic still fails if the kernel also has
+    # static shared arrays.
+    if needed > _DEFAULT_SHARED_MEM:
         kernel.max_dynamic_shared_size_bytes = int(dynamic_bytes)
 
 
diff --git a/src/cuperiod/methods/string_length.py b/src/cuperiod/methods/string_length.py
index cec36d5..c72d359 100644
--- a/src/cuperiod/methods/string_length.py
+++ b/src/cuperiod/methods/string_length.py
@@ -20,6 +20,7 @@
 from cuperiod.core._arrayapi import (
     array_namespace,
     device_ref,
+    is_torch_array,
     resolve_precision,
     resolve_torch_device,
     to_device_array,
@@ -63,14 +64,22 @@ def _length_batch(
     dev = device_ref(periods)
     n_periods = int(periods.shape[0])
     length = xp.empty(n_periods, dtype=periods.dtype, device=dev)
+    torch_input = is_torch_array(periods)
     for start in range(0, n_periods, batch):
         stop = min(start + batch, n_periods)
         pb = periods[start:stop]
         n_p = int(pb.shape[0])
         phase = xp.remainder(tau[None, :] / pb[:, None], 1.0)  # (P, N)
-        order = xp.argsort(phase, axis=1)
-        rows = xp.arange(n_p, dtype=idtype, device=dev)[:, None]
-        ph = phase[rows, order]            # phase sorted per row
+        if torch_input:
+            # torch.sort returns sorted values and the (stable) order in one kernel,
+            # replacing the argsort + row-gather pair.
+            import torch
+
+            ph, order = torch.sort(phase, dim=1, stable=True)
+        else:
+            order = xp.argsort(phase, axis=1)
+            rows = xp.arange(n_p, dtype=idtype, device=dev)[:, None]
+            ph = phase[rows, order]        # phase sorted per row
         mm = m_scaled[order]               # magnitudes gathered in the same order
         dphi = ph[:, 1:] - ph[:, :-1]
         dmag = mm[:, 1:] - mm[:, :-1]
@@ -82,6 +91,148 @@ def _length_batch(
     return length
 
 
+# --- GPU fast path: one CUDA block per trial period ---------------------------
+
+#: CUDA threads per block for the sort kernel (one block per trial period).
+CUDA_BLOCK: Final = 256
+
+#: Shared-memory bytes per point in the sort kernel: a float64 phase + int32 index.
+_SL_BYTES_PER_POINT: Final = 12
+
+#: One-block-per-period string-length kernel: a block folds its points, bitonic-sorts
+#: the (phase, original index) pairs in shared memory — comparing the index on equal
+#: phases makes the order *stable*, matching the array-API paths' stable argsort — and
+#: accumulates the string length from the sorted neighbours. No (P, N) intermediates
+#: and no global sort scratch. Fits light curves up to the shared-memory capacity
+#: (~4096 points at the default 48 KB); larger curves take the vectorized path.
+_SL_CUDA_SRC: Final = r"""
+extern "C" __global__ void sl_block(
+    const double* __restrict__ tau, const double* __restrict__ mag,
+    const double* __restrict__ periods,
+    const int n_points, const int n_pad, const int n_periods,
+    double* o_length)
+{
+    const int pidx = blockIdx.x;
+    if (pidx >= n_periods) return;
+    const int tid = threadIdx.x;
+    const int nth = blockDim.x;
+    const double period = periods[pidx];
+
+    extern __shared__ double sh[];
+    double* ph = sh;                    // (n_pad) folded phases
+    int* idx = (int*)(sh + n_pad);      // (n_pad) original indices
+
+    for (int i = tid; i < n_pad; i += nth) {
+        if (i < n_points) {
+            double q = tau[i] / period;
+            ph[i] = q - floor(q);       // mod(tau/period, 1), as the CPU paths
+            idx[i] = i;
+        } else {
+            ph[i] = 2.0;                // pad above any phase; sorts to the end
+            idx[i] = 0x7fffffff;
+        }
+    }
+    __syncthreads();
+
+    for (int k = 2; k <= n_pad; k <<= 1) {
+        for (int j = k >> 1; j > 0; j >>= 1) {
+            for (int i = tid; i < n_pad; i += nth) {
+                int ixj = i ^ j;
+                if (ixj > i) {
+                    bool up = ((i & k) == 0);
+                    double pa = ph[i], pb = ph[ixj];
+                    int ia = idx[i], ib = idx[ixj];
+                    bool greater = (pa > pb) || (pa == pb && ia > ib);
+                    if (greater == up) {
+                        ph[i] = pb; ph[ixj] = pa;
+                        idx[i] = ib; idx[ixj] = ia;
+                    }
+                }
+            }
+            __syncthreads();
+        }
+    }
+
+    double total = 0.0;
+    for (int i = tid; i < n_points - 1; i += nth) {
+        double dphi = ph[i + 1] - ph[i];
+        double dmag = mag[idx[i + 1]] - mag[idx[i]];
+        total += sqrt(dphi * dphi + dmag * dmag);
+    }
+    if (tid == 0) {
+        double dphi = (ph[0] + 1.0) - ph[n_points - 1];
+        double dmag = mag[idx[0]] - mag[idx[n_points - 1]];
+        total += sqrt(dphi * dphi + dmag * dmag);
+    }
+    __shared__ double r_t[CUDA_BLOCK];
+    r_t[tid] = total;
+    __syncthreads();
+    for (int s = blockDim.x >> 1; s > 0; s >>= 1) {
+        if (tid < s) r_t[tid] += r_t[tid + s];
+        __syncthreads();
+    }
+    if (tid == 0) o_length[pidx] = r_t[0];
+}
+"""
+
+_sl_kernel_cache: dict[int, Any] = {}
+
+
+def _sl_kernel(block: int) -> Any:
+    """Compile (once) and cache the string-length RawKernel for a block size."""
+    kernel = _sl_kernel_cache.get(block)
+    if kernel is None:
+        import cupy
+
+        src = _SL_CUDA_SRC.replace("CUDA_BLOCK", str(block))
+        kernel = cupy.RawKernel(src, "sl_block")
+        _sl_kernel_cache[block] = kernel
+    return kernel
+
+
+def _sl_cuda_capacity() -> int:
+    """Largest light curve the sort kernel can hold in opt-in shared memory."""
+    import cupy
+
+    optin = int(
+        cupy.cuda.Device().attributes.get("MaxSharedMemoryPerBlockOptin", 48 * 1024)
+    )
+    return (optin - CUDA_BLOCK * 8) // _SL_BYTES_PER_POINT  # static reduce buffer
+
+
+def _sl_cuda(
+    tau: FloatArray,
+    m_scaled: FloatArray,
+    periods: FloatArray,
+    *,
+    block: int = CUDA_BLOCK,
+) -> FloatArray:
+    """One-block-per-period CUDA string length; returns a host float64 array."""
+    import cupy as cp
+
+    n = int(tau.size)
+    n_pad = 1
+    while n_pad < n:
+        n_pad <<= 1
+    tau_d = cp.asarray(np.ascontiguousarray(tau, dtype=np.float64))
+    mag_d = cp.asarray(np.ascontiguousarray(m_scaled, dtype=np.float64))
+    per_d = cp.asarray(np.ascontiguousarray(periods, dtype=np.float64))
+    n_periods = int(per_d.size)
+    out = cp.empty(n_periods, dtype=cp.float64)
+    from cuperiod.core.backend import ensure_shared_memory
+
+    smem = n_pad * _SL_BYTES_PER_POINT
+    kernel = _sl_kernel(block)
+    ensure_shared_memory(kernel, smem, method="STRINGLENGTH", hint="n (light curve)")
+    kernel(
+        (n_periods,),
+        (block,),
+        (tau_d, mag_d, per_d, np.int32(n), np.int32(n_pad), np.int32(n_periods), out),
+        shared_mem=smem,
+    )
+    return np.asarray(cp.asnumpy(out), dtype=np.float64)
+
+
 # --- CPU fast path: numba-parallel, one loop-iteration per trial period -------
 
 _NUMBA_SL_KERNEL: Any = None
@@ -172,6 +323,10 @@ def string_length(
 
     if backend == "cupy":
         ensure_cuda_dll_path()
+        # The in-block sort kernel needs the whole curve in shared memory; longer
+        # curves fall back to the vectorized sort-based path.
+        if t.size <= _sl_cuda_capacity():
+            return _sl_cuda(tau, m_scaled, periods_host)
         import cupy as cp
 
         per_cp = cp.asarray(periods_host)

From 6153d76fda6b65cde529c2184422b7e8ba61659c Mon Sep 17 00:00:00 2001
From: Tharindu Jayasinghe <tharindujay96@gmail.com>
Date: Wed, 1 Jul 2026 14:43:26 -0700
Subject: [PATCH 8/9] test+docs: fast-backend parity suite, changelog, backend
 docs

- tests/test_fast_backends.py: numba-vs-numpy parity for PDM/CE/SL/
  MHAOV/TLS (SL with forced phase ties), float32 CUDA kernel checks
  (BLS tolerates isolated near-tie box flips but requires the same
  detected period), scatter-shim unit tests (numpy + torch), BLS
  device-cache reuse equivalence, and cpu->numba resolution for all six
  methods; new requires_numba marker.
- CHANGELOG Performance section; backends/installation/methods/index
  docs updated for the numba tier and float32 kernels; mypy override
  for cupyx.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
---
 CHANGELOG.md                      |  45 +++++++
 docs/guide/backends.md            |  34 ++---
 docs/guide/methods.md             |   3 +-
 docs/index.md                     |   2 +-
 docs/installation.md              |  18 +--
 pyproject.toml                    |   2 +
 src/cuperiod/methods/_bls_core.py |  10 +-
 src/cuperiod/methods/base.py      |   2 +-
 src/cuperiod/methods/tls.py       |   2 +-
 src/cuperiod/multiband/bls_mb.py  |   2 +-
 tests/conftest.py                 |   6 +
 tests/test_fast_backends.py       | 203 ++++++++++++++++++++++++++++++
 12 files changed, 301 insertions(+), 28 deletions(-)
 create mode 100644 tests/test_fast_backends.py

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 012dc95..d8c7c1c 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -6,6 +6,51 @@ All notable changes to cuPeriod are documented here. The format is based on
 
 ## [Unreleased]
 
+### Performance
+
+- **Multicore `numba` CPU kernels for PDM, CE, String-Length, MHAOV, and TLS** (BLS
+  already had one). With the `[fast]` extra installed, `backend="cpu"`/`"auto"` now
+  resolve to them; warm-kernel speedups over the vectorized numpy paths on a 3k-point
+  curve: PDM ~300x, CE ~135x, MHAOV ~57x, TLS ~53x, String-Length ~23x, with parity to
+  the numpy results at or below ~1e-11.
+- **MHAOV is rewritten around harmonic trig sums.** Every entry of the normal equations
+  is analytically a trig sum, so the Gram matrix is now assembled from `C_m`/`S_m`
+  sums computed with the Chebyshev recurrence (one `cos`/`sin` evaluation regardless of
+  the harmonic order) instead of materializing the `(F, N, 2H+1)` design tensor:
+  ~2.7-3.2x faster on numpy/cupy/torch alike, an order of magnitude less transient
+  memory (no more multi-GB tensors at 1e5 points), and gemm-free by construction —
+  the Blackwell cuBLAS workaround branch is gone because nothing calls gemm anymore.
+- **The portable GLS direct path evaluates `cos`/`sin` once instead of six times** —
+  the base-grid sums for the weights and the weighted data share one angle matrix and
+  the doubled-frequency sums follow from the double-angle identities: ~1.8x on
+  torch:cpu and torch:cuda; the frequency chunk is auto-capped by the light-curve
+  length so device memory stays bounded regardless of `N`.
+- **PDM bins each point once** into `n_bins*n_covers` fine bins and regroups every
+  cover exactly from that histogram (vectorized, CUDA, and numba paths): ~2.8x on the
+  numpy path, one third of the shared-memory atomics in the CUDA kernel.
+- **Opt-in `precision="float32"` now reaches the CUDA `RawKernel`s** (BLS, PDM, CE,
+  TLS): on consumer GPUs, whose float64 throughput is 1/64 of float32, the
+  FLOP-bound searches speed up ~8.6x (BLS full segmented run 610 -> 71 ms; TLS
+  34 -> 4 ms on an RTX Blackwell card). float64 stays the default; PDM keeps float64
+  accumulators and CE counts are exact integers at any precision. The float32 BLS
+  guards empty box windows with an `ivar` floor scaled to the total inverse variance
+  (float32 cumsum noise would otherwise fabricate boxes).
+- **String-Length gains a real CUDA kernel**: one block per period bitonic-sorts the
+  (phase, index) pairs in shared memory — stable, matching the CPU tie order — with
+  no `(P, N)` intermediates (~75x over numpy; curves longer than the shared-memory
+  capacity fall back to the vectorized path). The torch path fuses its sort+gather
+  via `torch.sort(stable=True)`.
+- **Faster scatter shims everywhere**: numpy binning goes through buffered
+  `bincount` instead of unbuffered `np.add.at`, cupy through `cupyx.scatter_add`,
+  and the batch kernels scatter row-wise so torch reads a stride-0 view instead of
+  materializing `(P, N)` broadcast copies and flat indices.
+- **BLS uploads each light curve to the GPU once per run** (cupy and torch) via a
+  per-curve device cache shared by the period segments — previously every segment
+  re-uploaded the arrays and synced on a device-side `t.min()` — and returns all
+  seven per-period outputs in one stacked device-to-host copy. The GLS NUFFT paths
+  batch the base-grid pair (`w`, `w*y`) into a single `n_trans=2` transform, and the
+  one-shot cufinufft path assembles the power on the GPU (194 -> 290 curves/s).
+
 ### Added
 
 - **Multi-vendor GPU support via PyTorch and the Python array API.** All seven
diff --git a/docs/guide/backends.md b/docs/guide/backends.md
index fa5c63d..103bb31 100644
--- a/docs/guide/backends.md
+++ b/docs/guide/backends.md
@@ -39,25 +39,25 @@ The selectors:
   - `cupy`
   - `numba` if installed, else `astropy`
 * - PDM
-  - `numpy`
+  - `numba` *(with `[fast]`)*, `numpy`
   - `cupy`
-  - `numpy`
+  - `numba` if installed, else `numpy`
 * - CE
-  - `numpy`
+  - `numba` *(with `[fast]`)*, `numpy`
   - `cupy`
-  - `numpy`
+  - `numba` if installed, else `numpy`
 * - String-Length
-  - `numpy`
+  - `numba` *(with `[fast]`)*, `numpy`
   - `cupy`
-  - `numpy`
+  - `numba` if installed, else `numpy`
 * - MHAOV
-  - `numpy`
+  - `numba` *(with `[fast]`)*, `numpy`
   - `cupy`
-  - `numpy`
+  - `numba` if installed, else `numpy`
 * - TLS
-  - `numpy`
+  - `numba` *(with `[fast]`)*, `numpy`
   - `cupy`
-  - `numpy`
+  - `numba` if installed, else `numpy`
 ```
 
 † BLS's `numpy` backend is a GPU-parity *reference*, not the product path — it shares one
@@ -172,11 +172,15 @@ How to read this:
   search). They're so fast on the CPU that the GPU's marginal gain on a single curve is
   modest — but the GPU still wins decisively for **large grids and big catalogs**
   ({doc}`batch`).
-- **PDM, CE, String-Length, MHAOV, TLS run on numpy on the CPU**, so the GPU's
-  data-parallelism delivers roughly 50–190× on a single curve. If you're searching many
-  periods or many stars with these, the GPU is a large win. The portable **`torch:cuda`**
-  column tracks alongside — competitive for these methods, and the same code reaches
-  AMD/Intel/Apple GPUs (see the portable-backend note above).
+- The **CPU times in the table are the vectorized numpy paths**. With the `[fast]` extra
+  installed, PDM, CE, String-Length, MHAOV, and TLS now run multicore **numba** kernels
+  on `"cpu"`/`"auto"` — one to two orders of magnitude faster than the numpy column
+  (e.g. PDM ~300×, CE ~135×, MHAOV ~57× on a 3k-point curve) — which puts a warm CPU
+  within reach of the GPU for a *single* curve. The GPU remains the throughput champion
+  for large grids and catalogs.
+- On **consumer NVIDIA cards** (GeForce), whose float64 throughput is 1/64 of float32,
+  the opt-in `precision="float32"` runs the BLS/TLS CUDA kernels ~8-9× faster at
+  detection-grade accuracy; float64 stays the default.
 - cuPeriod's **CPU** path already beats the established reference tools it was checked
   against (GLS ~3× astropy, PDM ~4× PyAstronomy, BLS ~18× astropy's `BoxLeastSquares`).
 
diff --git a/docs/guide/methods.md b/docs/guide/methods.md
index a75860a..2f7a489 100644
--- a/docs/guide/methods.md
+++ b/docs/guide/methods.md
@@ -83,7 +83,8 @@ for peak in pg.best_periods(5, alias_diverse=True):     # alias-aware for box se
 ```
 
 With the `[fast]` extra, BLS's CPU backend is a multicore `numba` search ~20× faster than
-astropy. Key settings ({class}`~cuperiod.BLSSettings`): `min_period_days` /
+astropy (PDM, CE, String-Length, MHAOV, and TLS gain `numba` CPU kernels too — see
+{doc}`backends`). Key settings ({class}`~cuperiod.BLSSettings`): `min_period_days` /
 `max_period_days`, `duration_min_frac` / `duration_max_frac`, `n_durations`, `objective`.
 
 ### TLS — transit least squares
diff --git a/docs/index.md b/docs/index.md
index ab68900..43e1b19 100644
--- a/docs/index.md
+++ b/docs/index.md
@@ -84,7 +84,7 @@ All seven share one API, one CLI, and the full single/batch machinery. See
 pip install cuperiod            # CPU (numpy, scipy, astropy, finufft)
 pip install "cuperiod[gpu]"     # + CUDA 12 GPU backends (cupy, cufinufft)
 pip install "cuperiod[torch]"   # + portable PyTorch backend (AMD/Intel/Apple GPUs + CPU)
-pip install "cuperiod[fast]"    # + numba multicore box search (~20× astropy BLS on CPU)
+pip install "cuperiod[fast]"    # + numba multicore CPU kernels (all methods, 20-300×)
 pip install "cuperiod[gui]"     # + interactive desktop GUI (cuperiod-gui)
 ```
 
diff --git a/docs/installation.md b/docs/installation.md
index 984ba47..2235d74 100644
--- a/docs/installation.md
+++ b/docs/installation.md
@@ -19,23 +19,27 @@ CPU, the command line, and batch processing over a process pool.
 | --- | --- | --- |
 | **gpu** | `pip install "cuperiod[gpu]"` | CUDA-12 GPU backends (`cupy-cuda12x`, `cufinufft`, the NVIDIA runtime wheels) |
 | **torch** | `pip install "cuperiod[torch]"` | the portable **PyTorch** backend — runs every method on AMD (ROCm), Intel (XPU), Apple (MPS), and a CPU path |
-| **fast** | `pip install "cuperiod[fast]"` | a multicore `numba` box search — BLS's CPU default, ~20× faster than astropy's compiled `BoxLeastSquares` |
+| **fast** | `pip install "cuperiod[fast]"` | multicore `numba` CPU kernels for BLS, PDM, CE, String-Length, MHAOV, and TLS — the CPU default when installed, one to two orders of magnitude faster than the fallback CPU paths |
 | **gui** | `pip install "cuperiod[gui]"` | the interactive desktop GUI, `cuperiod-gui` (PySide6 + pyqtgraph) — see {doc}`guide/gui` |
 | **pandas** | `pip install "cuperiod[pandas]"` | pandas `DataFrame` ingestion |
 
 Extras combine, e.g. `pip install "cuperiod[gpu,fast]"` or `"cuperiod[gui,torch]"`.
 
 :::{tip}
-The `[fast]` extra is worth installing even without a GPU: it makes BLS an order of
-magnitude faster on the CPU while matching astropy to floating point. When present it
-becomes BLS's default CPU backend automatically; otherwise BLS falls back to astropy.
+The `[fast]` extra is worth installing even without a GPU: it gives BLS, PDM, CE,
+String-Length, MHAOV, and TLS multicore JIT kernels that are one to two orders of
+magnitude faster than the fallback CPU paths while matching them to floating point.
+When present they become the default CPU backends automatically.
 :::
 
 ### What the `[fast]` extra changes
 
-Installing `numba` flips `backend="cpu"` for BLS from astropy's `BoxLeastSquares` to the
-in-house multicore box search. Nothing else about your code changes — the results match
-astropy to round-off; they just arrive ~20× sooner. (Note that `numba` currently caps
+Installing `numba` flips `backend="cpu"` (and the CPU fallback of `"auto"`) from the
+fallback implementations — astropy's `BoxLeastSquares` for BLS, the vectorized numpy
+kernels for PDM, CE, String-Length, MHAOV, and TLS — to in-house multicore JIT kernels.
+Nothing else about your code changes: the results match the fallbacks to round-off
+(BLS still matches astropy), they just arrive much sooner (~20× for BLS over astropy;
+~25–300× for the others over their numpy paths). (Note that `numba` currently caps
 `numpy < 2.5`, so installing it may downgrade numpy slightly.)
 
 ## GPU requirements
diff --git a/pyproject.toml b/pyproject.toml
index 4749bd0..7c6718d 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -162,6 +162,8 @@ module = [
     "cufinufft",
     "cupy",
     "cupy.*",
+    "cupyx",
+    "cupyx.*",
     "numba",
     "numba.*",
     "pandas",
diff --git a/src/cuperiod/methods/_bls_core.py b/src/cuperiod/methods/_bls_core.py
index c7c83f1..82b0b3f 100644
--- a/src/cuperiod/methods/_bls_core.py
+++ b/src/cuperiod/methods/_bls_core.py
@@ -628,7 +628,15 @@ def _bls_search_cuda(
             np.float64(data["t_min"]),
             np.int32(obj_flag),
             rdtype(-np.inf),
-            rdtype(_IVAR_EPS),
+            # In float64 an empty box window's ivar_in is exactly 0, so the absolute
+            # eps floor suffices. In float32 the cumsum difference of an empty window
+            # is cancellation *noise* of order sum_ivar*eps_f32 — an absolute 2.2e-16
+            # floor would admit garbage boxes — so the floor scales with the total.
+            rdtype(
+                max(_IVAR_EPS, data["sum_ivar"] * 8.0 * float(np.finfo(np.float32).eps))
+            )
+            if precision == "float32"
+            else np.float64(_IVAR_EPS),
             out["power"],
             out["depth"],
             out["depth_err"],
diff --git a/src/cuperiod/methods/base.py b/src/cuperiod/methods/base.py
index 65b6b4e..64c98e2 100644
--- a/src/cuperiod/methods/base.py
+++ b/src/cuperiod/methods/base.py
@@ -144,7 +144,7 @@ def resolve_backend(self, requested: str) -> str:
         return requested
 
     def _best_cpu_backend(self, available: set[str]) -> str:
-        """``fast_cpu_backend`` when its dependency is installed, else ``cpu_backend``."""
+        """``fast_cpu_backend`` when its dependency is present, else ``cpu_backend``."""
         if self.fast_cpu_backend is not None and self.fast_cpu_backend in available:
             return self.fast_cpu_backend
         return self.cpu_backend
diff --git a/src/cuperiod/methods/tls.py b/src/cuperiod/methods/tls.py
index d41cd33..25bcd45 100644
--- a/src/cuperiod/methods/tls.py
+++ b/src/cuperiod/methods/tls.py
@@ -558,7 +558,7 @@ def tls_power(
 
 
 class TLSMethod(PeriodogramMethod):
-    """Transit Least Squares — limb-darkened matched filter (numba/numpy CPU, cupy GPU)."""
+    """Transit Least Squares — limb-darkened matched filter (numba CPU, cupy GPU)."""
 
     name: ClassVar[str] = "TLS"
     objective_sense: ClassVar[Literal["max", "min"]] = "max"
diff --git a/src/cuperiod/multiband/bls_mb.py b/src/cuperiod/multiband/bls_mb.py
index 24ed3f6..0029ba8 100644
--- a/src/cuperiod/multiband/bls_mb.py
+++ b/src/cuperiod/multiband/bls_mb.py
@@ -84,7 +84,7 @@ def bls_multiband_power(
     band_caches: list[dict[str, object]] = [{} for _ in bands]
     for periods, durations in segments:
         per_band = []
-        for lc, device_cache in zip(bands, band_caches):
+        for lc, device_cache in zip(bands, band_caches, strict=True):
             err = lc.error if lc.error is not None else np.ones_like(lc.value)
             seg = _segment_power(
                 backend,  # type: ignore[arg-type]
diff --git a/tests/conftest.py b/tests/conftest.py
index 2954588..ad6a268 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -38,6 +38,12 @@ def _importable(name: str) -> bool:
     reason="torch (the [torch] extra) is not installed",
 )
 
+#: Skip a test unless numba (the ``[fast]`` extra) is importable.
+requires_numba = pytest.mark.skipif(
+    not _importable("numba"),
+    reason="numba (the [fast] extra) is not installed",
+)
+
 
 def _torch_gpu_available() -> bool:
     if not _importable("torch"):
diff --git a/tests/test_fast_backends.py b/tests/test_fast_backends.py
new file mode 100644
index 0000000..79be5f1
--- /dev/null
+++ b/tests/test_fast_backends.py
@@ -0,0 +1,203 @@
+"""Backend parity for the performance tiers: numba kernels, float32 CUDA, caches.
+
+The per-method test modules pin each statistic against a transparent reference; this
+module pins the *fast* backends against the vectorized numpy path — the numba CPU
+kernels for all five ported methods, the float32-templated CUDA kernels, the shared
+scatter shims, and the BLS per-light-curve device cache.
+"""
+
+from __future__ import annotations
+
+import numpy as np
+
+import cuperiod as cup
+from conftest import requires_gpu, requires_numba, requires_torch
+from cuperiod.core._arrayapi import scatter_add_rows, scatter_counts_rows
+from cuperiod.methods._bls_core import bls_power
+from cuperiod.methods.conditional_entropy import conditional_entropy
+from cuperiod.methods.mhaov import aov_power
+from cuperiod.methods.pdm import pdm_theta
+from cuperiod.methods.string_length import string_length
+from cuperiod.methods.tls import tls_power
+from synth import synthetic_eclipser, synthetic_sine
+
+PERIODS = np.linspace(0.3, 3.0, 400)
+FREQS = np.linspace(0.4, 3.0, 400)
+
+
+def _tied_sine(n: int = 300) -> tuple[np.ndarray, np.ndarray]:
+    """A sine curve with duplicated times, so phase folds tie exactly."""
+    t, mag, _ = synthetic_sine(n=n)
+    t[50:60] = t[49]
+    return t, mag
+
+
+# --- numba CPU kernels vs the vectorized numpy path ---------------------------
+
+
+@requires_numba
+def test_numba_pdm_matches_numpy() -> None:
+    t, mag, _ = synthetic_sine(n=300)
+    cpu = pdm_theta(t, mag, PERIODS, backend="numpy")
+    fast = pdm_theta(t, mag, PERIODS, backend="numba")
+    assert np.max(np.abs(cpu - fast)) < 1e-12
+
+
+@requires_numba
+def test_numba_ce_matches_numpy() -> None:
+    t, mag, _ = synthetic_sine(n=300)
+    cpu = conditional_entropy(t, mag, PERIODS, backend="numpy")
+    fast = conditional_entropy(t, mag, PERIODS, backend="numba")
+    assert np.max(np.abs(cpu - fast)) < 1e-12
+
+
+@requires_numba
+def test_numba_string_length_matches_numpy_with_ties() -> None:
+    # Duplicated times force equal phases: the numba mergesort must break ties in
+    # the same (stable) order as the array-API paths' stable argsort.
+    t, mag = _tied_sine()
+    cpu = string_length(t, mag, PERIODS, backend="numpy")
+    fast = string_length(t, mag, PERIODS, backend="numba")
+    assert np.max(np.abs(cpu - fast)) < 1e-9
+
+
+@requires_numba
+def test_numba_mhaov_matches_numpy() -> None:
+    t, mag, _ = synthetic_sine(n=300)
+    cpu = aov_power(t, mag, FREQS, n_harmonics=3, backend="numpy")
+    fast = aov_power(t, mag, FREQS, n_harmonics=3, backend="numba")
+    assert np.allclose(cpu, fast, rtol=1e-8, atol=1e-8)
+
+
+@requires_numba
+def test_numba_tls_matches_numpy() -> None:
+    t, flux, err = synthetic_eclipser(period=2.5)
+    periods = np.linspace(1.5, 4.0, 300)
+    settings = cup.TLSSettings()
+    cpu = tls_power(t, flux, err, periods, settings=settings, backend="numpy")
+    fast = tls_power(t, flux, err, periods, settings=settings, backend="numba")
+    for key in ("sde", "sr", "depth", "duration", "t0"):
+        assert np.allclose(cpu[key], fast[key], rtol=1e-10, atol=1e-10), key
+
+
+@requires_numba
+def test_cpu_request_resolves_to_numba() -> None:
+    for name in ("PDM", "CE", "STRINGLENGTH", "MHAOV", "TLS", "BLS"):
+        assert cup.get_method(name).resolve_backend("cpu") == "numba"
+
+
+# --- float32 CUDA kernels ------------------------------------------------------
+
+
+@requires_gpu
+def test_bls_cupy_float32_close_to_float64() -> None:
+    t, flux, err = synthetic_eclipser(period=2.5)
+    periods = np.linspace(1.5, 4.0, 500)
+    durations = np.asarray([0.05, 0.1])
+    p64 = bls_power(t, flux, err, periods, durations, 10, backend="cupy",
+                    precision="float64")
+    p32 = bls_power(t, flux, err, periods, durations, 10, backend="cupy",
+                    precision="float32")
+    # Detection-grade agreement: float32 rounding may flip which of two near-tied
+    # boxes wins at isolated periods, so allow a small fraction of outliers but
+    # require the same detected period.
+    close = np.isclose(p64.power, p32.power, rtol=1e-3, atol=1e-3)
+    assert close.mean() > 0.99
+    assert int(np.argmax(p64.power)) == int(np.argmax(p32.power))
+
+
+@requires_gpu
+def test_pdm_ce_tls_cupy_float32_close_to_float64() -> None:
+    t, mag, _ = synthetic_sine(n=400)
+    theta64 = pdm_theta(t, mag, PERIODS, backend="cupy", precision="float64")
+    theta32 = pdm_theta(t, mag, PERIODS, backend="cupy", precision="float32")
+    assert np.allclose(theta64, theta32, rtol=5e-3, atol=5e-3)
+
+    h64 = conditional_entropy(t, mag, PERIODS, backend="cupy", precision="float64")
+    h32 = conditional_entropy(t, mag, PERIODS, backend="cupy", precision="float32")
+    assert np.allclose(h64, h32, rtol=5e-3, atol=5e-3)
+
+    tt, flux, err = synthetic_eclipser(period=2.5)
+    periods = np.linspace(1.5, 4.0, 300)
+    r64 = tls_power(tt, flux, err, periods,
+                    settings=cup.TLSSettings(precision="float64"), backend="cupy")
+    r32 = tls_power(tt, flux, err, periods,
+                    settings=cup.TLSSettings(precision="float32"), backend="cupy")
+    assert int(np.argmax(r64["sr"])) == int(np.argmax(r32["sr"]))
+
+
+@requires_gpu
+def test_string_length_gpu_kernel_handles_ties() -> None:
+    # The in-block bitonic sort breaks equal phases by original index — the same
+    # stable order as the CPU paths.
+    t, mag = _tied_sine()
+    cpu = string_length(t, mag, PERIODS, backend="numpy")
+    gpu = string_length(t, mag, PERIODS, backend="cupy")
+    assert np.max(np.abs(cpu - gpu)) < 1e-9
+
+
+# --- scatter shims ---------------------------------------------------------------
+
+
+def test_scatter_add_rows_numpy_matches_manual() -> None:
+    rng = np.random.default_rng(0)
+    index = rng.integers(0, 7, size=(4, 50))
+    values = rng.normal(size=50)
+    target = np.zeros((4, 7))
+    scatter_add_rows(target, index, values)
+    expected = np.zeros((4, 7))
+    for p in range(4):
+        np.add.at(expected[p], index[p], values)
+    assert np.allclose(target, expected)
+
+    counts = np.zeros((4, 7))
+    scatter_counts_rows(counts, index)
+    expected_counts = np.zeros((4, 7))
+    for p in range(4):
+        np.add.at(expected_counts[p], index[p], 1.0)
+    assert np.allclose(counts, expected_counts)
+
+
+@requires_torch
+def test_scatter_add_rows_torch_matches_numpy() -> None:
+    import torch
+
+    rng = np.random.default_rng(1)
+    index = rng.integers(0, 9, size=(3, 40))
+    values = rng.normal(size=40)
+    np_target = np.zeros((3, 9))
+    scatter_add_rows(np_target, index, values)
+    t_target = torch.zeros((3, 9), dtype=torch.float64)
+    scatter_add_rows(t_target, torch.as_tensor(index), torch.as_tensor(values))
+    assert np.allclose(np_target, t_target.numpy())
+
+    np_counts = np.zeros((3, 9))
+    scatter_counts_rows(np_counts, index)
+    t_counts = torch.zeros((3, 9), dtype=torch.float64)
+    scatter_counts_rows(t_counts, torch.as_tensor(index))
+    assert np.allclose(np_counts, t_counts.numpy())
+
+
+# --- BLS device cache -------------------------------------------------------------
+
+
+@requires_torch
+def test_bls_device_cache_reuse_matches_fresh() -> None:
+    t, flux, err = synthetic_eclipser(period=2.5)
+    seg1 = np.linspace(1.5, 2.4, 200)
+    seg2 = np.linspace(2.4, 4.0, 200)
+    durations = np.asarray([0.05, 0.1])
+    cache: dict[str, object] = {}
+    with_cache = [
+        bls_power(t, flux, err, seg, durations, 10, backend="torch:cpu",
+                  device_cache=cache)
+        for seg in (seg1, seg2)
+    ]
+    fresh = [
+        bls_power(t, flux, err, seg, durations, 10, backend="torch:cpu")
+        for seg in (seg1, seg2)
+    ]
+    assert "torch" in cache  # the upload happened once and was kept
+    for got, want in zip(with_cache, fresh, strict=True):
+        assert np.allclose(got.power, want.power)
+        assert np.allclose(got.transit_time, want.transit_time)

From 51c92c60319418f8394662f19d7522d701d06861 Mon Sep 17 00:00:00 2001
From: Tharindu Jayasinghe <tharindujay96@gmail.com>
Date: Wed, 1 Jul 2026 14:47:05 -0700
Subject: [PATCH 9/9] ci: install EGL/OpenGL system libraries for the Linux GUI
 tests

The ubuntu-latest image stopped shipping libEGL.so.1, so pytest-qt now
fails at configure importing PySide6 QtGui (INTERNALERROR before any
test runs). Unrelated to the code changes on this branch.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
---
 .github/workflows/ci.yml | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 1ad1a8f..df51c29 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -21,6 +21,11 @@ jobs:
         python: ["3.11", "3.12"]
     steps:
       - uses: actions/checkout@v4
+      # PySide6 (the GUI tests, run offscreen) needs the EGL/OpenGL system
+      # libraries, which the ubuntu-latest image no longer preinstalls.
+      - name: Install Qt system libraries (Linux)
+        if: runner.os == 'Linux'
+        run: sudo apt-get update && sudo apt-get install -y libegl1 libopengl0 libgl1
       - uses: astral-sh/setup-uv@v5
       - name: Create environment
         run: uv venv --python ${{ matrix.python }}