[oscardssmith]

Stacked on #30.

What

• One kernel for all large prime radices. A stage of a prime radix from 17 on runs prime_stage! (src/primestage.jl): the radix and a table of cosines and sines are arguments, and the kernel is compiled once per vector type. It covers the stages vectorised over the lanes, the last stage, lane blocks and the grouped first stage. The radices up to 16 keep their straight-line codelets.
• MAX_RADIX_PRIME is 73 (was 31). The KernelAbstractions kernels and the generic number types keep 31; factorize takes the limit as an argument.
• The lazy machinery is gone: the callable tables per stage kind and the lazy_radices preference.
• No @generated stage kernels. Every combination the chain can select has a plain method. One table (STAGE_KERNELS) lists the kernels and their arguments; the methods and the dispatch switches are made from it. stage! and stage_range! share one body builder, and the scalar tail of the grouped first stage is the grouped kernel on a vector of one complex value. The generated bodies of the radices 2–16 are unchanged (compared for every kind, element type, sign, vector type and radix).

Why

A straight-line codelet of a prime r has about r²/2 statements with as many constants (525 statements for 31) and needed 82 kernel methods per radix on AVX-512, which is why 17–31 were compiled at first use. With the loop a prime costs a table.

Measurements

AMD EPYC 9354 (AVX-512), Julia 1.13.3, one thread, plan_fft and mul!, rounds in alternating processes, median time over the median time on #30:

set	cases	geometric mean	above 1.05
lengths and shapes with a factor 17–31	40	0.82	1 (Float64 100 × 31: 1.06)
long lengths with such a factor (31·2⁶ … 17·2²⁰, times primes above 1000)	36	0.94	1 (Float64 19·2¹⁷: 1.08)
the primes 37–73, against Bluestein before	100	0.29	0
general cases without such factors	88	1.00	3 (Float64 lines of 32, 48, 60: 1.07–1.15)

• The three general cases were at 1.04, 0.83 and 1.05 in a second run of 6 rounds on 20 short lines (geometric mean 1.00). The stage kernels of those lengths compile to the same machine code as on One chain form: every smooth length runs on RChain #30; I did not find a cause in the callers.
• The primes 37–73 alone: 0.43–0.79 (Float64), 0.48–0.90 (Float32); as a factor of a longer length: 0.13–0.53.
• Error against a 200-bit reference, Float64, lengths 841, 6256, 6859: 1.4–1.8 ε, with the codelets 3.5–4.6 ε.

Sizes on the same machine:

	#30	this PR
libfftl.so (x86-64-v3, --trim=safe)	49.5 MB	26.8 MB
code in libfftl.so	18.6 MB	8.8 MB
package image	58.6 MB	66.2 MB

The package image grows because the kernel and its tables are in the image and the lazy codelets were not; 1.4 MB of that is typed code of an inlined helper (see below).

The limit

Measured with the limit at 127: as a factor of a longer length a prime from 79 on is still 1.3–4× faster than Bluestein, but the line of the prime alone is slower from 79 on for Float32 (1.05; 1.5 at 127) and from 97 on for Float64. 73 is the largest prime for which every measured case is faster.

What made the loop as fast as unrolled code

• Up to eight output pairs per pass over the inputs (four with 16 vector registers).
• 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; 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, was up to 1.4× slower for a whole transform. The recursion is what stores the 1.4 MB of typed code.
• The calls of the kernel sit behind a function of their own, outside the dispatch switches; inside, they cost every case about 1 ns.

Tests

• Quick selection with 2 threads: 10956 pass, 10 broken as before. New cases for the primes 37–73 (alone, first, last, twice, with a scalar remainder) and for relinking a chain with a looped radix.
• Tests that used 37, 41 and 43 as examples of primes without a radix use 79, 83 and 89, so that the Bluestein and composite paths keep their coverage.
• The radix table test keeps its set of lengths (prime factors up to 31), so the per-platform files are unchanged.
• libfftl builds; fftw3_ctest.sh passes (3940 comparisons, 0 failed).

Not done

• aarch64 is not verified, and the GPU path was not run (its limit is unchanged).
• The loop for 11 and 13: an earlier per-radix form of it was at 0.98 (Float64) and 1.05–1.06 (Float32) of the codelets; not measured with the single kernel.
• A larger limit for lengths that are not a prime alone.
• The full weekly test lists and the documentation build were not run.

Decision and evidence: docs/src/explanation/design-decisions.md ("The prime radices from 17 to 73 share one looped kernel"). AGENTS.md no longer says that @generated methods stay as a fallback.

🤖 Generated with Claude Code

[tknopp]

Checked on aarch64 (Apple M1 Max, same setup as in #30, time over main). The quick selection passes with 2 threads (10956 pass, 10 broken as before), so the PRIME_BLOCK = 8 path on aarch64 is correct.

case	#31 / main
Float64 17 / 31	0.71 / 0.75
Float64 100 × 31	1.03
Float32 19 · 512	1.11
Float64 37 / 67 / 73, Float32 41	0.36 / 0.52 / 0.54 / 0.52
Float64 37 · 64, Float32 71 · 128	0.36 / 0.43
general cases (16 … 4096, planes)	0.89–1.13, as in #30

The primes 37–73 confirm your numbers. The one case against the x86 measurement is Float32 19 · 512 at 1.11; perhaps PRIME_BLOCK = 8 is not the best choice on NEON (32 registers, but only 128 bit wide). Not a blocker; worth a look when someone tunes aarch64.

The image growth to 66 MB is a fair price for no compilation on first use and half the libfftl.so. Please make sure image-size.md says so explicitly, since it described the 32 MB of the lazy radices.

[oscardssmith]

I think I may have more image size improvements. (right now we compile all eagerly and precompile maximally and neither likely makes sense)