One looped kernel for the prime radices 17 to 73

AGENTS.md

@@ -26,9 +26,10 @@ Read this file before changing the package. The documentation holds the rest:
   before it is written down, or it is marked as not verified. A change that claims to be faster comes
   with an interleaved A/B measurement on one machine (medians, and the spread of the same code against
   itself); see `docs/src/explanation/performance.md`.
-- **Kernels of a fixed set are plain methods, not `@generated`.** A package image keeps the complete
-  optimised IR of every precompiled `@generated` specialisation. `@generated` methods stay as the
-  fallback for combinations outside the precompiled set. See `docs/src/explanation/image-size.md`.
+- **Kernels are plain methods, not `@generated`.** A package image keeps the complete optimised IR
+  of every precompiled `@generated` specialisation. The stage kernels are defined with `@eval` for
+  the fixed set of radices, signs and vector types; there is no `@generated` fallback. See
+  `docs/src/explanation/image-size.md`.
 - **Nothing private in the repository.** Code, documentation, data files, commit messages and pull request
   descriptions contain no host names or accounts of machines, no local paths, no private contact data and no
   names of people beyond the authors of the package. A machine is named by its processor ("AMD EPYC 7702",

README.md

@@ -84,7 +84,7 @@ 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
+  dimensions, for any length (the part of a length made of prime factors above 73 through
   Bluestein's algorithm).
 - `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
@@ -117,19 +117,20 @@ figures there are regenerated from the committed data on every documentation bui
 ## How it works
 
 - **Generated codelets.** A small symbolic generator writes straight-line code for DFTs of
-  size 2–31 (composite sizes split in registers, odd primes in the symmetric form). Codelets are
-  generated with `@generated` for the radix, element type and SIMD width actually used.
+  size 2–16 (composite sizes split in registers, odd primes in the symmetric form), one method
+  per radix, element type and SIMD width. The prime radices 17–73 share one kernel that loops
+  over a table of cosines and sines.
 - **Complex numbers in SIMD registers.** On aarch64 CPUs with the ARMv8.3 complex-number
   extension (all Apple silicon), additions with a factor ±i and complex multiplications use the
   `FCADD` and `FCMLA` instructions; elsewhere the same operations are built from shuffles.
 - **Stockham chains.** A transform is a sequence of radix stages (mostly 4 and 8) without bit
-  reversal. The sequence is part of the plan's type, so it compiles into one function. The radices
+  reversal. The sequence is data of the plan; nothing is compiled per length. The radices
   are chosen by a fixed rule found in offline searches over the decompositions of many lengths;
   there is no run-time measuring step.
 - **Multidimensional transforms** run the same chains along every dimension: a dimension with
   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
+- **Lengths with prime factors above 73** split off their smooth part in one Cooley–Tukey step
   and use Bluestein's algorithm only for the rest.
 - **Threads** split each dimension over slices or blocks of columns; long 1D transforms use the
   four-step algorithm.

docs/src/explanation/algorithms.md

@@ -18,13 +18,13 @@ which keeps the data in natural order throughout:
 The radices are chosen by a fixed rule (`src/factor.jl`). The rule is the result of exhaustive
 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
+into a 12; the x86 rules differ and are in the same file. Every other prime up to 73 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).
 
 ## Codelets
 
-A stage of radix ``r`` is a small DFT applied to ``r`` values at once. For each ``r`` the codelet
+A stage of radix ``r`` is a small DFT applied to ``r`` values at once. For each ``r`` up to 16 the codelet
 is straight-line code without loops, generated by a symbolic generator (`src/codegen.jl`). The
 generator writes the DFT as a sum of terms and simplifies it:
 
@@ -38,9 +38,15 @@ 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`).
 
+The prime radices from 17 to 73 use the same symmetric form, but as loops: one kernel
+(`src/primestage.jl`) takes the radix and a table of the cosines and sines as arguments. Straight-line
+code for such a radix has about ``r^2/2`` statements with as many different constants, which no longer
+fit in registers; the loop computes up to eight output pairs in one pass over the inputs, so that the
+inputs are loaded once for all of them.
+
 ## 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``,
+Lengths with a prime factor above 73 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}.

docs/src/explanation/design-decisions.md

@@ -58,7 +58,8 @@ showed 100 `Float32` and 64 `Float64` lengths slower than 1.5×, the worst at 3.
 
 ## Lengths with large prime factors: split off the smooth part, Bluestein on the rest
 
-**Decision (2026-11-07).** For n = s·q, with q the part made of primes above 31: one Cooley–Tukey
+**Decision (2026-11-07).** For n = s·q, with q the part made of primes above `MAX_RADIX_PRIME` (73, at the
+time of the decision 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
@@ -115,8 +116,11 @@ for the rest), and before that the idea of removing `Chain` without changing the
 
 ## Kernels of a fixed set are plain methods, not `@generated`
 
-**Decision (2026-11-10).** The stage kernels of the codelet set that `RChain` dispatches to are plain
-methods defined in `src/kernel_methods.jl` from the bodies the `@generated` methods build. The `@generated` methods remain as the fallback.
+**Decision (2026-11-10, fallback removed 2026-11-13).** The stage kernels of the codelet set that `RChain`
+dispatches to are plain methods, defined with `@eval` in `src/kernel_methods.jl` from the functions that
+build their bodies. One table (`STAGE_KERNELS` in `src/kernels.jl`) lists the kernels and their arguments;
+the methods and the dispatch switches of `RChain` are made from it. There are no `@generated` stage
+kernels: every combination the chain can select has a plain method, so a fallback would never run.
 
 **Why.** A package image keeps the complete optimised IR of every precompiled `@generated`
 specialisation, 130 MB of the 146 MB image. See [Package size and load time](@ref "Package size and load time").
@@ -149,16 +153,64 @@ radix 16 middle stage 284 → 198, radix 8 100 → 79. One line, `plan_fft`, med
 rounds: lengths 512–16384 take 0.79–0.94 (`Float32`) and 0.82–0.92 (`Float64`) of the time, lengths
 from 32768 on are unchanged (the same code against itself differs by 1–5 %). Not verified on aarch64.
 
-## Codelets of the radices 17–31 are compiled on first use
-
-**Decision (2026-11-10).** The codelets of the radices 17, 19, 23, 29 and 31 are not part of the package
-image. The first call compiles them (70–100 ms); later calls go through a table of callable objects at a
-cost of about 0.02 µs per transform.
-
-**Why.** These radices were 54 % of the statements of the stage kernels and are rarely needed. The image
-shrinks from 50 MB to 32 MB. Making 11 and 13 lazy as well would save 4 MB more; it was not done, because
-lengths with those factors are common. The list is the preference `lazy_radices`; the trimmed C library
-sets it to empty, because `juliac` compiles nothing at run time.
+## The prime radices from 17 to 73 share one looped kernel
+
+**Decision (2026-11-13).** A stage of a prime radix from 17 on runs `prime_stage!` (`src/primestage.jl`):
+one kernel for all such radices, with the radix and a table of cosines and sines as arguments, compiled
+once per vector type. It uses the symmetric form of the odd-prime codelet as loops and computes up to
+eight output pairs (four with 16 vector registers) in one pass over the inputs. It covers the stage
+kinds that the codelets have separate kernels for: stages vectorised over the lanes, the last stage,
+lane blocks, and the grouped first stage. `MAX_RADIX_PRIME` is 73. The radices up to 16 keep their
+straight-line codelets. The KernelAbstractions kernels and the generic number types keep their limit
+of 31.
+
+**Why.** A straight-line codelet of a prime r has about r²/2 statements with as many different
+constants, which do not fit in registers (525 statements for 31), and every radix needed 82 kernel
+methods on AVX-512 (kind × sign × vector type). The codelets of 17–31 were therefore compiled at first
+use and not with the package (decision of 2026-11-10, replaced by this one). With the loop, a prime
+costs a table.
+
+Measurements: AMD EPYC 9354 (AVX-512), one thread, `plan_fft` and `mul!`, rounds in alternating
+processes, median time over the median time of the previous code.
+
+- Lengths and shapes with a factor 17–31 against the codelets (40 cases, 4 rounds): geometric mean
+  0.82, 36 cases below 0.95, one above 1.05 (`Float64` 100 × 31, 1.06: a single stage in place, which
+  copies the inputs of a butterfly first). 17 and 31 alone: 0.63–0.84; 1088 = 17·64: 0.96 and 0.89;
+  837 = 27·31 and 841 = 29²: 0.72–0.82; 62 × 62: 0.70–0.74.
+- Long lengths (36 cases, 3 rounds; 31·2⁶ to 17·2²⁰, and a factor 17–31 times a prime above 1000):
+  geometric mean 0.94, one case above 1.05 (`Float64` 19·2¹⁷, 1.08). 31·2⁶ to 31·2¹⁴: 0.74–0.91;
+  17·2²⁰: 0.99 and 1.01.
+- The primes 37–73 against Bluestein's algorithm, which they used before (100 cases, 2 rounds; p, 8p,
+  64p, 1024p, 16384p): geometric mean 0.29. The prime alone: 0.43–0.79 (`Float64`) and 0.48–0.90
+  (`Float32`); as a factor of a longer length: 0.13–0.53.
+- 92 general cases without such factors (the set of the decision above, 4 rounds): geometric mean 1.00
+  without the two cases with a dimension of 37 (0.10 and 0.16), 3 cases above 1.05 (`Float64` lines of
+  32, 48 and 60: 1.07–1.15). In 6 rounds on 20 short lines alone these three were at 1.04, 0.83 and
+  1.05 and the geometric mean 1.00; the stage kernels of these lengths compile to the same machine code
+  as before.
+- Accuracy against a 200-bit reference, `Float64`, 841, 6256 and 6859: 1.4–1.8 ε, with the codelets
+  3.5–4.6 ε (the table is built from angles reduced to one turn).
+
+What made the loop as fast as the unrolled sums: several rows per pass, so that the inputs are loaded
+once for all of them and the multiply-add chains of the rows run side by side; no buffer for the sums
+and differences of the inputs (storing whole registers was the slowest part: 1.3–1.4× per butterfly
+with the buffer); the pass over a block of rows out of line and the outputs of a block written by
+recursion over the tuples of accumulators (a loop with an index into them, or a closure over them, put
+the accumulators on the stack: up to 1.4× for a whole transform).
+
+**The limit.** With the limit at 127 for the measurement, a prime from 79 on is still 1.3–4× faster than
+Bluestein's algorithm as a factor of a longer length, but the line of that prime alone is slower from
+79 on for `Float32` (1.05, 1.5 at 127) and from 97 on for `Float64` (1.01, 1.15 at 127). 73 is the
+largest prime for which every measured case is faster.
+
+**Not done.** The loop for 11 and 13: with one kernel per radix (an earlier form of the loop) it was at
+0.98 of the codelets for `Float64` and 1.05–1.06 for `Float32`, and 1.08–1.10 for 7; not measured with
+the single kernel. A larger limit for lengths that are not a prime alone. aarch64 is not verified.
+
+**Costs.** The kernel and its tables are in the package image, which the lazy codelets were not: the
+image file is 66 MB instead of 59 MB on this machine (68.5 MB before the decision on the chain form).
+The C library contains every kernel in either case and shrinks: `libfftl.so` (x86-64-v3, `--trim=safe`)
+is 26.8 MB instead of 49.5 MB, its code 8.8 MB instead of 18.6 MB.
 
 ## Run-time pass plans for trimmed programs
 

docs/src/explanation/image-size.md

@@ -42,35 +42,32 @@ later Julia versions, so re-run the count below after an upgrade.
 
 ## What keeps it small
 
-1. **Plain methods for the fixed set.** The stage kernels that `RChain` dispatches to (radix × sign ×
-   vector type) are defined in `src/kernel_methods.jl` with `@eval`, from the same body-building
-   functions the generators call. The `@generated` methods stay as the fallback for every other combination
-   (the lazy radices), so behaviour is unchanged. A plain method keeps its (compressed) source
-   and its machine code, not its optimised IR.
-2. **Lazy radices.** The codelets of the radices in `LAZY_RADICES` (17, 19, 23, 29, 31 by default; the
-   preference `lazy_radices`, read in `src/rchain.jl`) are 54 % of the statements of the stage kernels and are rarely needed. They are not
-   in the image. A call with such a codelet index reaches `cstage_*_lazy`, which finds a callable object
-   for the (radix, sign, vector type) in a table and creates it at the first call; Julia compiles its
-   kernel then. A first call of a length with such a factor takes 70–100 ms; later calls cost about
-   0.02 µs more than a precompiled codelet. Calling the stage kernel with run-time `Val` arguments
-   instead of the table cost 0.22 µs per transform, 10 % at n = 1088.
-
-   `juliac` compiles nothing at run time, so a program that is compiled with `--trim` sets the preference
-   `lazy_radices = []`: every codelet is then compiled into it. `libfftl/LocalPreferences.toml` does this
-   for the C library; a project that builds FFTL into its own trimmed library sets it in its
-   `LocalPreferences.toml` (the project must have FFTL as a direct dependency):
-
-   ```toml
-   [FFTL]
-   lazy_radices = []
-   ```
-
-To add a radix, put it in `RADICES`; `kernel_methods.jl` defines its plain methods. To make a radix lazy,
-add it to the default of the preference in `src/rchain.jl`.
-
-The numbers on this page were measured while the smooth lengths up to 128 had compiled chains of their
-own (583 more plain methods). Without them the image on an AMD EPYC 9354 (AVX-512, which has more vector
-types than the M1) went from 72 MB to 61 MB; the M1 numbers have not been measured again.
+1. **Plain methods for the fixed set.** The stage kernels (`STAGE_KERNELS` in `src/kernels.jl`: radix ×
+   sign × vector type) are defined in `src/kernel_methods.jl` with `@eval`, from the functions that build
+   their bodies. A plain method keeps its (compressed) source and its machine code, not its optimised IR.
+   There are no `@generated` stage kernels.
+2. **One kernel for the large prime radices.** The radices 17 to 73 have no codelets of their own: one
+   kernel (`src/primestage.jl`) takes the radix and a table as arguments and is compiled once per vector
+   type. Straight-line codelets for 17, 19, 23, 29 and 31 were 54 % of the statements of the stage kernels;
+   they used to be left out of the image and compiled at the first call (70–100 ms, the "lazy radices"
+   column of the table above).
+
+To add a radix up to 16, put it in `RADICES`; `kernel_methods.jl` defines its methods. A larger prime
+needs no code: raise `MAX_RADIX_PRIME` (see [Design decisions](@ref "Design decisions") for where the
+kernel stops being faster than Bluestein's algorithm).
+
+The numbers on this page were measured on the M1 while the smooth lengths up to 128 had compiled chains
+of their own (583 more plain methods) and the radices 17–31 were compiled at first use. On an AMD EPYC
+9354 (AVX-512, which has seven vector types where the M1 has three) the image file went from 68.5 MB
+with the compiled chains to 58.6 MB without them and to 66.2 MB with the kernel for the large primes
+(which is in the image, where the lazy codelets were not). The C library, which contains every kernel
+in either case, went from 55.6 MB to 49.5 MB and to 26.8 MB (`libfftl.so`, x86-64-v3, `--trim=safe`). The M1 numbers have not been measured again.
+
+With the workload, nearly all of the image is the stage kernels: the dispatch switches of `RChain` name
+every kernel, so the first transform compiles all of them. On the EPYC 9354 (measured with the prime
+limit at 127, which does not change the kernels) the image was 65 MB with the workload and 20.5 MB
+without, precompilation took 61 s and 13 s, `using FFTL` 185 ms and 110 ms, and the first transform
+16–37 ms with the workload and 23–35 s without.
 
 ## What is left
 
@@ -88,7 +85,6 @@ Further steps that were considered:
 | step | gain | cost |
 |:--|:--|:--|
 | methods generic in the element and vector type instead of one per vector type | about 4 MB (stored sources ÷ 3) | the compiler must fold `ncomplex` and `sizeof` of the type parameters; needs a speed check |
-| radices 11 and 13 lazy as well | 4 MB (28 MB), a few ms of load time (inside the noise) | 70–100 ms at the first call of lengths with those factors, which are common |
 | no PrecompileTools | about 25 ms of load time (Dates, TOML, Preferences) | the workload is written with `precompile` calls |
 | no `Float32` in the workload | 16 MB (at 73 MB) | 6 s at the first `Float32` call |
 

docs/src/explanation/implementation.md

@@ -68,6 +68,11 @@ is a few nanoseconds of arithmetic, so the loop is written to cost little next t
 - The switch over the codelet index is inlined into the loop: a stage is one jump and one call.
 - The rows of a plane run inside one call of the loop function.
 
+A stage of a prime radix from 17 on calls one kernel for all such radices (`prime_stage!` in
+`src/primestage.jl`), with the radix and its table as arguments; it is compiled once per vector type.
+Its calls sit behind a function of their own, outside the switch, because they take more arguments
+than fit in registers and would cost every case of the switch a larger stack frame.
+
 The addresses in the stage records are valid for the chain object they were taken for. A chain
 remembers its own address at that moment; a copy of a plan, or a plan stored in a package image and
 loaded into another process, is at another address and links its records again before it uses them.
@@ -76,7 +81,7 @@ loaded into another process, is at another address and links its records again b
 at build time and compiles nothing at run time; the fixed codelet set is all a chain of any length needs. The
 [libfftl how-to](@ref "C library (libfftl)") describes what the library is built from, and
 [Design decisions](@ref "Design decisions") gives the measurements behind the single chain form. The stage kernels of
-the fixed codelet set are plain methods, not `@generated` ones, to keep the package image small; see
+the fixed codelet set are plain methods defined with `@eval`, not `@generated` ones, to keep the package image small; see
 [Package size and load time](@ref "Package size and load time").
 
 ## Lengths and the factorisation rule