Oversample instead of Bluestein for non-nice sizes

README.md

@@ -84,8 +84,8 @@ yd = FFTL.fft(MtlArray(randn(ComplexF32, 1024, 1024)))
 **Fourier and discrete cosine transforms written in Julia, with wide support for number types and CPU and GPU execution.**
 
 - Complex, real-input and all eleven real-to-real kinds of FFTW, along any subset of
-  dimensions, for any length (the part of a length made of prime factors above 31 through
-  Bluestein's algorithm).
+  dimensions, for any length (the part of a length made of prime factors above 31 resampled onto
+  a smooth grid).
 - `Float32` and `Float64` on the fast path (SIMD codelets, threads). `BigFloat`, `Double64`,
   `Float128` and other floating-point types are transformed in their own precision by a generic
   path, accurate to a few ulps of the type.
@@ -130,7 +130,8 @@ figures there are regenerated from the committed data on every documentation bui
   stride S is a chain started with stride S, the S neighbouring columns forming the vectorised
   index. No transposition.
 - **Lengths with prime factors above 31** split off their smooth part in one Cooley–Tukey step
-  and use Bluestein's algorithm only for the rest.
+  and resample only the rest onto a smooth grid (one transform of twice its length and a short
+  interpolation).
 - **Threads** split each dimension over slices or blocks of columns; long 1D transforms use the
   four-step algorithm.
 

docs/AGENTS.md

@@ -34,7 +34,7 @@ rules there apply to everything written here as well. Everything in this folder
   plans `RPlan`, `RRFFTPlan`, `RR2RPlan`), contribute.
 - **Explanation** (`src/explanation/`): the DFT and its conventions, FFTW and FFTL (history, papers, what
   FFTL takes from FFTW; written from papers and documentation only), algorithms (Stockham chains,
-  codelets, Bluestein, multidimensional without transposition, four-step, real and real-to-real
+  codelets, resampling, multidimensional without transposition, four-step, real and real-to-real
   reductions), implementation (type specialisation, SIMD, kernel cache, cache sizes, the two chain forms,
   why libfftl needs `RChain`), design decisions (what was decided and why, with the evidence), package
   size and load time, performance (measurement protocol).
@@ -51,9 +51,9 @@ Design decisions page. These concern the manual itself:
 
 - **Interface in the text:** the manual writes the exported names without prefix; next to `using FFTL` it
   uses `import AbstractFFTs`, because `using AbstractFFTs` would make `fft` ambiguous.
-- **Feature/platform table:** built from the README, the code and the libfftl header. The Bluestein
-  convolution-length rule is taken from the code (aarch64: smallest 5-smooth M ≥ 2n−1; x86: cost model,
-  7-smooth M in [2n−1, 1.4·(2n−1)]); the FFTW-padded real layout of in-place r2c/c2r in libfftl comes
+- **Feature/platform table:** built from the README, the code and the libfftl header. The grid-length
+  rule of resampled lengths is taken from the code (aarch64: smallest 5-smooth m ≥ 2n; x86: cost model,
+  7-smooth m in [2n, 2.8·n]); the FFTW-padded real layout of in-place r2c/c2r in libfftl comes
   from the header.
 - **Figures:** CairoMakie in the tutorials. Literate runs with `execute = false` (its default for
   markdown), so the scripts execute as Documenter `@example` blocks during `makedocs`; a figure cell ends

docs/src/explanation/algorithms.md

@@ -19,8 +19,8 @@ The radices are chosen by a fixed rule (`src/factor.jl`). The rule is the result
 searches over the decompositions of many lengths on the development machines, not of a runtime
 search. It prefers radix 4 and radix 8 for powers of two. On aarch64 it merges a factor 3 with a 4
 into a 12; the x86 rules differ and are in the same file. Every other prime up to 31 is a radix of
-its own. Larger primes are not radices; the part of the length they make up goes to Bluestein's
-algorithm (below).
+its own. Larger primes are not radices; the part of the length they make up is resampled onto a
+smooth grid (below).
 
 ## Codelets
 
@@ -38,26 +38,37 @@ The generated code runs on complex SIMD vectors (`src/complexvec.jl`), with the
 the element type. On aarch64 the complex add-and-rotate and the complex multiply are single
 instructions (`FCADD`, `FCMLA`).
 
-## Bluestein's algorithm
-
-Lengths with a prime factor above 31 go through Bluestein's identity. With ``jk = (j^2 + k^2 - (k-j)^2)/2``,
-
-```math
-X_k = c_k \sum_{j} (x_j c_j)\, \overline{c_{k-j}}, \qquad c_j = e^{\sigma \pi i j^2/n}.
-```
-
-The sum is a cyclic convolution with a kernel that depends only on ``n``. Its length ``M``
-is a smooth number with ``M \ge 2n - 1``, so the convolution is two transforms of smooth lengths
-and one pointwise product. The kernel's transform is computed once when the plan is made.
-
-Per line, Bluestein costs two transforms of length ``M`` and a pointwise product, so it is used
-only for the part of a length that the smooth codelets cannot cover. On aarch64 the convolution
-length is the smallest 5-smooth number with ``M \ge 2n - 1``. On x86 a cost model fitted to
-measured chains picks the fastest 7-smooth length in ``[2n - 1,\, 1.4\,(2n - 1)]``: a length with
-many radix 9, 5 or 3 stages transforms several times slower than a power of two a few percent
-larger (`src/plan.jl`).
-
-On x86 a convolution of 8 MiB or more (``M \ge 2^{19}`` in `Float64`) runs as a
+## Lengths with large prime factors: resampling
+
+Lengths with a prime factor above 31 are transformed by oversampled interpolation onto a smooth
+grid (a type-2 nonuniform FFT evaluated at the points ``k\,m/n``; `src/resample.jl`):
+
+1. the input, divided by the Fourier transform of an interpolation kernel ``\psi``, is
+   zero-padded to a smooth length ``m \ge 2n``,
+2. one transform of length ``m``,
+3. output ``k`` is the sum of the ``w`` grid values around ``k\,m/n``, weighted by
+   ``\psi(z) = e^{\beta(\sqrt{1 - z^2} - 1)}`` (the exponential of a semicircle).
+
+The result is not exact: its error is about ``e^{-\pi w \sqrt{1 - n/m}}`` relative to ``\lVert x \rVert``.
+FFTL takes the narrowest ``w`` (a multiple of 4) for an error of ``10\,\varepsilon``: 16 taps in
+`Float64`, 8 in `Float32`. The weights of every output are tabulated up to a table size (1 MiB
+in `Float64`, 4 MiB in `Float32`, where the two measured about even), above it they are
+evaluated from piecewise polynomials. Generic number types (`BigFloat` and others,
+`src/generic.jl`) keep Bluestein's algorithm, which is exact: there the width grows with the
+precision and the interpolation would cost more than the transform it saves.
+
+Compared with Bluestein's algorithm, which needs two transforms of a length ``\ge 2n - 1``, one
+transform is saved for ``n\,w`` multiply-adds. Measured against it on a Tiger Lake core (one
+thread, before Bluestein was removed), the time ratio was 0.51–0.82 for primes from 1009 to
+100003, at an error of 2–5·10⁻¹⁵ in `Float64` (Bluestein about 5·10⁻¹⁶) and 1–2·10⁻⁷ in
+`Float32` (Bluestein about 2.5·10⁻⁷).
+
+The grid length is a smooth number. On aarch64 it is the smallest 5-smooth ``m \ge 2n``. On x86
+a cost model fitted to measured chains picks the fastest 7-smooth length in ``[2n,\, 2.8\,n]``: a
+length with many radix 9, 5 or 3 stages transforms several times slower than a power of two a
+few percent larger (`grid_length` in `src/plan.jl`).
+
+On x86 a grid of 8 MiB or more (``m \ge 2^{19}`` in `Float64`) is transformed by a
 [four-step](@ref "Long one-dimensional transforms: four-step") transform on the calling thread
 instead of a chain. A chain of that length no longer fits in cache and makes one pass over
 memory per stage; the four-step makes two. On a Zen 2 core this is 1.06× faster at ``2^{19}``
@@ -67,7 +78,7 @@ and 1.3× from ``2^{20}`` on. On the M1 the chain is faster at every length meas
 ### Lengths with a smooth part
 
 A length ``n = s\,q`` whose large prime factors make up ``q`` and whose small ones make up
-``s > 1`` does not need a convolution of length ``2n``. Since ``s`` and ``q`` are coprime, one
+``s > 1`` does not need a grid of length ``2n``. Since ``s`` and ``q`` are coprime, one
 Cooley–Tukey step with ``j = j_2 + q\,j_1`` and ``k = k_1 + s\,k_2`` splits the transform:
 
 ```math
@@ -76,26 +87,30 @@ X_{k_1 + s k_2} = \sum_{j_2=0}^{q-1} e^{\sigma 2\pi i\, j_2 k_2/q}\;
 ```
 
 1. A chain of length ``s`` over the ``q`` interleaved lines ``j_2`` (the inner sum),
-2. multiplication by the twiddle factors ``e^{\sigma 2\pi i\, j_2 k_1/n}``, folded into the first
-   chirp product of the next step,
-3. ``s`` Bluestein transforms of length ``q`` (the outer sum), each written with stride ``s``.
-
-The arithmetic is about the same as one Bluestein transform of length ``n``: ``s`` convolutions
-of length about ``2q`` are again about ``2n`` points. The gain is where the data lives. A
-convolution of length ``2q`` fits in cache, one of length ``2n`` often does not. On a Zen 2 core
-(one thread, `Float64`) ``262145 = 65 \cdot 4033`` takes 6.6 ms instead of 16.8 ms, and
-``1048577 = 17 \cdot 61681`` takes 43 ms instead of 89 ms.
+2. multiplication by the twiddle factors ``e^{\sigma 2\pi i\, j_2 k_1/n}``, folded into the
+   deconvolution of the next step,
+3. ``s`` resampled transforms of length ``q`` (the outer sum), each written with stride ``s``.
+
+From 32 lanes on, the ``s`` transforms of the last step run as one batched pass (as in FFTJ):
+the grid holds rows of lanes, lane fastest, so the transform of length ``m`` is one chain over
+the lanes and every output row is ``w`` multiply-adds of whole grid rows. Lanes go in chunks
+whose grid stays in cache. A grid of length about ``2q`` fits in cache where one of length
+``2n`` often does not.
+
+With fewer than 32 lanes each lane is one resampled transform of its own, vectorised along the
+line. A batched pass vectorises across the lanes, so few lanes, or a lane count that is not a
+multiple of the vector width, leave most of the vector unused: on an EPYC 9354 (AVX-512, one
+thread, ``q`` = 37, 1009, 10007) the batched pass took 1.7–3.8 times as long as the lanes one by
+one at ``s = 2`` in `Float64` (2.1–7.3 in `Float32`), 1.1–1.6 (1.6–2.9) at ``s = 7``, 0.73–1.04
+(0.84–1.31) at ``s = 16`` and 0.56–1.04 (0.87–0.96) at ``s = 32`` (`RS_BATCH_LANES` in
+`src/resample.jl`). The rule depends on the length only, not on the number of threads, so a
+threaded plan computes the same values as a serial one.
 
 When ``q`` itself is large it is split again: ``q`` becomes the full power of its largest prime,
-and the first step of length ``s = n/q`` is again a split transform (or a Bluestein transform
-when ``s`` has no smooth part). ``3660001 = 73 \cdot 181 \cdot 277`` runs as Bluestein transforms of
-length 13213 and 277. Each split costs one more pass over memory, so a leftover ``q`` stays one
-Bluestein transform as long as its convolution takes at most 4 MiB. Below that, splitting
-measured 1.15–1.5× slower; above it, 1.02–1.4× faster (`composite_split` in `src/plan.jl`).
-
-Rader's algorithm, which turns a prime length ``p`` into a cyclic convolution of length
-``p - 1``, would halve the convolution where ``p - 1`` is smooth. Of 45 random lengths with a
-large prime factor between ``2^{10}`` and ``2^{22}``, five qualify; FFTL does not implement it.
+and the first step of length ``s = n/q`` is again a split transform (or a resampled transform
+when ``s`` has no smooth part). Each split costs one more pass over memory, so a leftover ``q``
+stays one resampled transform as long as its grid takes at most 4 MiB (`composite_split` in
+`src/plan.jl`; the threshold was measured with Bluestein's algorithm in place of resampling).
 
 ## Multidimensional transforms without transposition
 

docs/src/explanation/design-decisions.md

@@ -43,7 +43,7 @@ line, all buffers at fixed page offsets, time of the radix-16 sequence over the
 
 ## Rules for lengths are derived on some lengths and checked on others
 
-**Decision (2026-11-04).** The radix choice and the Bluestein length are rules on the prime
+**Decision (2026-11-04).** The radix choice and the grid length of a resampled transform are rules on the prime
 factorisation, not tables per length. A rule is derived from one part of the lengths (a brute-force
 search over all decompositions, or a cost model) and checked on lengths that were not used for the
 derivation, for example by parity of n or by a random split; only the checked result counts as
@@ -56,23 +56,44 @@ not single points.
 the M1. A dense sweep of every length 2…1024 against FFTW MEASURE (one thread, Zen 2, 2046 cases)
 showed 100 `Float32` and 64 `Float64` lengths slower than 1.5×, the worst at 3.4×.
 
-## Lengths with large prime factors: split off the smooth part, Bluestein on the rest
+## Lengths with large prime factors: split off the smooth part, resample the rest
 
 **Decision (2026-11-07).** For n = s·q, with q the part made of primes above 31: one Cooley–Tukey
-step, a chain of length s, twiddles, and s Bluestein transforms of length q (`Composite`). A leftover q
-stays one Bluestein while its convolution is at most 4 MiB; above that, q is split into its prime
-powers, recursively. Bluestein convolutions of 8 MiB and more run as a serial four-step algorithm on
-x86 (1.3× from 2²⁰ on, Zen 2) but not on aarch64, where the chain is faster up to 2²² (M1 Max). With
-threads, `CompositeMT` is used when s is at least the number of threads and is bitwise equal to the
-serial kernel; otherwise the whole length is one `BluesteinMT`.
-
-**Why.** The convolutions of length about 2q stay in cache, and the arithmetic equals one Bluestein of
-length n. On 45 non-smooth lengths (Zen 2) FFTL moved from ×0.64 to ×1.04 of the fastest other library
-with one thread, and from ×0.77 to ×1.16 with 16 threads.
+step, a chain of length s, twiddles, and s resampled transforms of length q (`Composite`): one batched
+pass from 32 lanes on, lane by lane below (2026-11-12). A leftover q stays one resampled transform while its grid is at most 4 MiB; above that,
+q is split into its prime powers, recursively. Grids of 8 MiB and more are transformed by a serial
+four-step algorithm on x86 (1.3× from 2²⁰ on, Zen 2) but not on aarch64, where the chain is faster up
+to 2²² (M1 Max). With threads, `CompositeMT` is used when s is at least the number of threads and is
+bitwise equal to the serial kernel; otherwise the whole length is one `ResampleMT`.
+
+**Why.** The grids of length about 2q stay in cache. On 45 non-smooth lengths (Zen 2) FFTL moved from
+×0.64 to ×1.04 of the fastest other library with one thread, and from ×0.77 to ×1.16 with 16 threads
+(measured with Bluestein's algorithm on the leftover). Resampling needs one transform of a length of
+at least 2q where Bluestein's algorithm needs two: it took 0.51–0.82 of Bluestein's time for primes
+from 1009 to 100003 (Tiger Lake, one thread), at an error of 2–5·10⁻¹⁵ in `Float64` (Bluestein about
+5·10⁻¹⁶). The batched pass vectorises across the lanes and loses with few of them: against the lanes one
+by one it took 1.7–3.8 times as long at s = 2, 1.1–1.6 at s = 7, 0.73–1.04 at s = 16 and 0.56–1.04 at
+s = 32 (`Float64`, q = 37, 1009, 10007, EPYC 9354, one thread; `Float32` 2.1–7.3, 1.6–2.9, 0.84–1.31,
+0.87–0.96). The switch is on the lane count and not on the number of threads, which keeps serial and
+threaded plans bitwise equal.
+
+The grid of a single resampled transform has a margin on each side that holds copies of the other end
+(2026-11-12), so every output is interpolated from one contiguous load. Without it the outputs whose
+window wraps around the grid ends (7 per `Float64` transform, 3 per `Float32`) took a scalar loop at
+about ten times the cost of the others: 250 of 559 ns at n = 37 and 290 of 1104 ns at n = 101. With the
+margin, time over the time without it: n = 37: 0.55, 101: 0.75, 127: 0.79, 1009: 0.96, 100003: 0.98,
+101³: 0.78; against `main`: 0.81, 0.81, 0.77, 0.64, 0.80, 0.83, and 0.85 for 101³ on 8 threads (1.17
+before). `Float64`, EPYC 9354, minima of interleaved runs, the same code against itself within 6 %.
+The margin is a whole page: a margin of 128 bytes moves the grid off the page offset of the work
+buffers, and the prime length 100003 then took 1.03–1.05 of the time without margins.
+
+**Replaced.** Bluestein's algorithm on the leftover and `BluesteinMT` for a whole line on threads
+(2026-11-07). Generic number types and the GPU passes keep Bluestein's algorithm.
 
 **Not done / not checked.** Rader's algorithm (5 of the 45 sweep lengths qualify). Both thresholds were
 measured on Zen 2 only, and the 4 MiB rule is checked on 12 lengths that also took part in its
-derivation.
+derivation. The 32-lane rule was measured on one AVX-512 machine only, not on AVX2 or aarch64, and the grid
+margin on x86 only (512-bit and 256-bit vectors), not on aarch64.
 
 ## Two chain forms: compiled for short smooth lengths, run-time data for the rest
 
@@ -144,7 +165,7 @@ sets it to empty, because `juliac` compiles nothing at run time.
 
 **Decision (2026-11-10).** A pass plan (`plan_pass`, `plan_r2c_pass`, `plan_c2r_pass`) built with
 `_runtime` (`plan_pass_runtime`, `plan_r2c_pass_runtime`, `plan_c2r_pass_runtime`) uses run-time kernels
-only: RChain, Bluestein and composite transforms over RChain, the four-step algorithm over RChain. Its
+only: RChain, resampled and composite transforms over RChain, the four-step algorithm over RChain. Its
 type depends on the element type and the number of dimensions, not on the length. The specialised plans
 are unchanged.
 
@@ -158,8 +179,8 @@ and matches the direct DFT to 3·10⁻¹⁴ when called from C. Planning on 64 
 specialised and the run-time plans is within noise of each other on the same views.
 
 **Costs.** The run-time kernels dispatch per stage, which costs a little on short lines (the reason the
-specialised plans stay the default). For real lengths with large prime factors, the run-time kernel is a
-Bluestein transform, without the composite split of the specialised plan.
+specialised plans stay the default). For real lengths with large prime factors, the run-time kernel is one
+resampled transform, without the composite split of the specialised plan.
 
 ## The interface: FFTW.jl's names, FFTL's own functions, AbstractFFTs as a backend
 

docs/src/explanation/implementation.md

@@ -31,8 +31,8 @@ is how the KernelAbstractions kernels use the codelet generator without a second
 ## Kernel cache
 
 Kernels are cached by kind, element type, length, sign and radix sequence. The cache has 128
-entries and drops the least recently used one when it is full. A Bluestein kernel is two entries,
-its chain and its tables. `clear_cache!` empties it. Within one plan, dimensions with the same
+entries and drops the least recently used one when it is full. A resampling kernel is several
+entries: its chain and its tables. `clear_cache!` empties it. Within one plan, dimensions with the same
 length and stride class share a kernel object.
 
 ## Cache sizes and blocks

docs/src/explanation/performance.md

@@ -50,6 +50,7 @@ application they happen at different times.
 - DCT-I and DST-I are computed through a real DFT of about twice the length, and DCT-IV of odd
   length likewise. The extra length shows in the measurements.
 - On the GPU, the real-to-real kinds are not tuned yet.
-- Bluestein transforms a line through two transforms of length ``M ≥ 2q - 1``, where ``q`` is
-  at most the part of the length made of prime factors above 31 (the smooth part is split off,
-  see [Algorithms](@ref "Lengths with a smooth part")). A prime length pays this in full.
+- A length with prime factors above 31 is resampled: one transform of a smooth length
+  ``m ≥ 2q`` and ``w`` multiply-adds per output, where ``q`` is at most the part of the length
+  made of those prime factors (the smooth part is split off, see
+  [Algorithms](@ref "Lengths with a smooth part")). A prime length pays this in full.