Straight-line codelets for the 8/16/32/64-point base cases of the radix-4 kernel

[pankgeorg]

Item C of the plan in #130, stacked on #134 (twiddle tables; the diff includes that commit until it lands).

Before, fft_pow2_radix4! recursed down to 2- and 4-point base cases, so roughly half of the work on cache-resident powers of two was call and index overhead; FFTW's advantage at these sizes is its library of generated straight-line codelets.

Now, for Complex{Float32} and Complex{Float64} the recursion stops at 64 points (or the whole transform if smaller) in a fully unrolled radix-2 DIT codelet produced by an @generated function keyed on the size and direction — every intermediate in its own SSA variable, twiddles folded to constants, trivial multiplications removed, explicit fma so rounding does not depend on which array type the codelet is specialised for. This is what FFTW's genfft emits offline; here Julia's compiler does it. Other element types (BigFloat, Float16, symbolic) keep the generic recursion, which stays the reference implementation.

Compile latency. Each (size, element type, direction) codelet costs 0.1–0.6 s of LLVM time on first use, so the package gets a PrecompileTools workload that compiles the 8/16/32/64-point codelets and the common plan/execute paths at package precompile time (new dependency, PrecompileTools = "1"). Net effect on a fresh session: using FFTA; fft(randn(ComplexF64, 4096)) goes from ~1.9 s to well under a millisecond. Standalone, the 8/16/32/64-point codelets run 3.3×/1.6×/2.2×/1.3× faster than the recursion they replace (14.6/47/125/334 ns vs FFTW's 24.8/40.6/80.2/308 ns on this machine).

Before/after (aarch64 Neoverse-N1, Julia 1.12.6, benchmark/suite.jl, planned execution, single thread; "before" is #134, FFTA/FFTW vs FFTW 3.3.11 ESTIMATE):

class	kind	type	cases	FFTA speedup geomean (min–max)	FFTA/FFTW before → after	max bytes/exec before → after	plan time geomean before → after
1d/awkward	fft	Float32	29	1.17× (0.92–1.53)	3.20× → 2.74×	0 → 0	1.2 ms → 1.1 ms
1d/awkward	fft	Float64	29	1.14× (1.02–1.48)	2.55× → 2.25×	0 → 0	1.5 ms → 1.5 ms
1d/pow2	fft	Float32	20	1.39× (1.05–2.22)	2.32× → 1.68×	0 → 0	35.0 µs → 34.3 µs
1d/pow2	fft	Float64	40	1.37× (1.00–2.42)	1.96× → 1.43×	0 → 0	48.1 µs → 44.9 µs
1d/prime	fft	Float32	19	1.21× (0.99–1.59)	3.32× → 2.71×	0 → 0	109.5 µs → 96.0 µs
1d/prime	fft	Float64	19	1.17× (0.96–1.57)	2.37× → 2.04×	0 → 0	135.2 µs → 121.4 µs
1d/smooth	fft	Float32	24	1.07× (0.97–1.26)	4.97× → 4.60×	0 → 0	70.3 µs → 68.6 µs
1d/smooth	fft	Float64	24	1.07× (0.99–1.23)	3.23× → 3.06×	0 → 0	88.9 µs → 88.7 µs
2d	fft	Float32	9	1.38× (1.07–1.77)	4.16× → 3.05×	32 KiB → 32 KiB	2.2 µs → 2.1 µs
2d	fft	Float64	14	1.32× (1.09–1.75)	2.86× → 2.19×	64 KiB → 64 KiB	6.9 µs → 6.4 µs
3d	fft	Float32	5	1.44× (1.20–1.73)	8.17× → 5.47×	2 KiB → 2 KiB	1.3 µs → 1.2 µs
3d	fft	Float64	5	1.37× (1.04–1.65)	5.39× → 3.91×	4 KiB → 4 KiB	1.6 µs → 1.5 µs
batched_dim1	fft	Float32	6	1.33× (1.22–1.55)	3.25× → 2.47×	0 → 0	10.0 µs → 9.9 µs
batched_dim1	fft	Float64	6	1.29× (1.21–1.49)	2.04× → 1.60×	0 → 0	12.6 µs → 12.7 µs
batched_dim2	fft	Float32	6	1.18× (1.06–1.41)	2.14× → 1.81×	0 → 0	9.9 µs → 9.7 µs
batched_dim2	fft	Float64	6	1.12× (0.88–1.36)	1.77× → 1.54×	0 → 0	12.2 µs → 12.5 µs

Largest slowdowns / speedups (planned execution, FFTA before → after; FFTW for reference):

• 0.88× — fft Float64 64×1024 dims=(2,): 1.51 ms → 1.71 ms (FFTW 1.12 ms)
• 0.91× — fft Float64 64×4096 dims=(2,): 8.16 ms → 9.01 ms (FFTW 5.50 ms)
• 0.92× — fft Float32 49233 dims=(1,): 9.15 ms → 9.92 ms (FFTW 3.11 ms)
• 0.96× — fft Float64 7 dims=(1,): 0.1 µs → 0.1 µs (FFTW 0.0 µs)
• 2.12× — fft Float32 32 dims=(1,): 0.3 µs → 0.1 µs (FFTW 0.1 µs)
• 2.21× — fft Float64 32 dims=(1,): 0.3 µs → 0.1 µs (FFTW 0.1 µs)
• 2.22× — fft Float32 8 dims=(1,): 0.1 µs → 0.0 µs (FFTW 0.0 µs)
• 2.29× — fft Float64 32 dims=(1,): 0.3 µs → 0.1 µs (FFTW 0.1 µs)
• 2.37× — fft Float64 8 dims=(1,): 0.1 µs → 0.0 µs (FFTW 0.0 µs)
• 2.42× — fft Float64 8 dims=(1,): 0.1 µs → 0.0 µs (FFTW 0.0 µs)

261 matched cases; geometric-mean speedup 1.22×; 3 cases slower by >5%.

Against main the combination #134 + this PR is a 2.35× geometric-mean speedup over the 495 comparable cases of the suite. Three strided-batched cases are within 10% either way (noise-level on this shared host).

Tests: the existing suite (exact-equality checks between strided and contiguous execution included) passes; a 3D-rfft-unsupported error count in the benchmark is the known limitation, not new.

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 98.90%. Comparing base (7aeb327) to head (5b02c9d).
⚠️ Report is 4 commits behind head on main.

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #135      +/-   ##
==========================================
+ Coverage   98.80%   98.90%   +0.09%     
==========================================
  Files           5        6       +1     
  Lines         585      729     +144     
==========================================
+ Hits          578      721     +143     
- Misses          7        8       +1     

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[wheeheee]

I suspect copyto! is often more optimized for large N

[wheeheee]

examples for when LLVM chooses not to contract for some array types? might have seen something similar before but cannot remember.

in general I think muladd is still the way to go, it accommodates machines with no native fma (emulated fma is slow). IMHO granular fast-math flags would be exactly what's needed here, but the issue for that on the Julia repo is still open...

[wheeheee]

maybe better to use traits?