Store twiddle tables and Bluestein data in the plan

[pankgeorg]

Item A of the plan in #130: twiddle factors and Bluestein data are computed once at plan time and stored in the CallGraph, instead of being regenerated on every execution.

Before, every kernel seeded Singleton's recurrence with singleton_params (a sincospi) — once per output row of the O(n²) DFT leaf, once per j1 in fft_composite!, and per recursion level of the radix-4/3 kernels — and fft_bluestein! allocated three pad-length buffers and recomputed the chirp and its FFT on every call. For n = 5 the trig calls were the entire cost of the transform; for n = 1000 about 1 800 sincospi per execution.

Now

• CallGraph gets twiddles::Vector{Vector{T}} (one table per node), bluestein::Vector{BluesteinScratch{T}} + blue_index, and dir::Direction (tables are direction specific; fft! throws if called with the other direction — plans always have a fixed direction).
• Table layouts (documented in src/callgraph.jl): DFT nodes store w^k; composite nodes store the (j1, k2) block in the order the kernel walks it; the radix-4 and radix-3 kernels store per-level interleaved (w^k, w^2k, w^3k) / (w^k, w^2k) triplets/pairs and address them with a flat offset passed down the recursion. All tables are derived from one unit_roots table per node, which evaluates sincospi only on the first octant when 8 | N (so planning a 2^20 transform costs ~130k trig calls, 8 ms, not 1.4 M).
• BluesteinScratch holds the chirp, its transform (pre-scaled by 1/pad_len), the two work arrays and the pow2 tables for the padded length; the convolution is written with forward transforms only (ifft(y) = conj(fft(conj(y)))/n) so a single table set suffices. Planned execution is now allocation-free for every size, including primes and composites with a Bluestein factor.
• The old kernel signatures (fft_dft!(..., d::Direction) etc.) remain as thin wrappers that build the table on the fly, so direct callers (the test suite, and the usage suggested in README benchmark result cannot be reproduced #119) keep working.
• Accuracy: tables are correctly rounded, so the Float32 error no longer grows with n: at n = 2^22 the relative error vs a Float64 reference drops from ~1000 ulp to 1.5 ulp (and from 10–28 ulp to 1.2–2.2 ulp across the whole 2^16…2^22 and 3^9…3^11 grid — FFTW is at ~1.5 ulp) (the accuracy test grid is extended to 2^20 and 2^22).

Cost: plan creation now builds the tables — ~26 µs at 4096, ~0.4 ms at 65536, ~8 ms at 2^20, ~70 ms at 2^22 (ComplexF64); memory per plan grows by roughly one table of n entries per node level. One-shot fft(x) (plan + execute) is still faster than before at every size in the sweep.

Before/after (aarch64 Neoverse-N1, Julia 1.12.6, benchmark/suite.jl, planned execution, single thread; FFTA/FFTW is 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	2.37× (2.08–2.62)	7.45× → 3.20×	48 MiB → 0	0.6 µs → 1.2 ms
1d/awkward	fft	Float64	29	2.26× (1.91–2.57)	5.63× → 2.55×	96 MiB → 0	0.7 µs → 1.5 ms
1d/pow2	fft	Float32	20	1.63× (1.39–1.85)	3.77× → 2.32×	0 → 0	0.1 µs → 35.0 µs
1d/pow2	fft	Float64	40	1.60× (1.22–1.88)	3.08× → 1.96×	0 → 0	0.1 µs → 48.1 µs
1d/prime	fft	Float32	19	2.28× (1.41–3.10)	7.54× → 3.32×	96 MiB → 0	0.1 µs → 109.5 µs
1d/prime	fft	Float64	19	2.51× (1.98–3.66)	5.81× → 2.37×	192 MiB → 0	0.1 µs → 135.2 µs
1d/smooth	fft	Float32	24	2.12× (1.38–2.76)	10.18× → 4.97×	0 → 0	0.8 µs → 70.3 µs
1d/smooth	fft	Float64	24	2.34× (1.38–3.62)	7.33× → 3.23×	0 → 0	0.9 µs → 88.9 µs
2d	fft	Float32	9	1.47× (1.38–1.59)	6.06× → 4.16×	32 KiB → 32 KiB	0.2 µs → 2.2 µs
2d	fft	Float64	14	1.83× (1.42–2.90)	5.18× → 2.86×	9 MiB → 64 KiB	0.2 µs → 6.9 µs
3d	fft	Float32	5	1.43× (1.37–1.50)	11.63× → 8.17×	2 KiB → 2 KiB	0.3 µs → 1.3 µs
3d	fft	Float64	5	1.49× (1.43–1.54)	7.56× → 5.39×	4 KiB → 4 KiB	0.3 µs → 1.6 µs
batched_dim1	fft	Float32	6	1.57× (1.38–1.65)	5.07× → 3.25×	0 → 0	0.1 µs → 10.0 µs
batched_dim1	fft	Float64	6	1.50× (1.32–1.66)	3.08× → 2.04×	0 → 0	0.1 µs → 12.6 µs
batched_dim2	fft	Float32	6	1.34× (1.14–1.55)	2.90× → 2.14×	0 → 0	0.1 µs → 9.9 µs
batched_dim2	fft	Float64	6	1.31× (1.07–1.60)	2.24× → 1.77×	0 → 0	0.1 µs → 12.2 µs

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

• 1.07× — fft Float64 64×65536 dims=(2,): 324.29 ms → 302.69 ms (FFTW 269.61 ms)
• 1.14× — fft Float32 64×16384 dims=(2,): 49.58 ms → 43.64 ms (FFTW 28.96 ms)
• 1.15× — fft Float64 64×16384 dims=(2,): 62.99 ms → 54.93 ms (FFTW 34.18 ms)
• 1.16× — fft Float32 64×65536 dims=(2,): 279.49 ms → 241.41 ms (FFTW 194.20 ms)
• 3.10× — fft Float32 101 dims=(1,): 15.8 µs → 5.1 µs (FFTW 1.9 µs)
• 3.11× — fft Float64 101 dims=(1,): 16.6 µs → 5.3 µs (FFTW 2.4 µs)
• 3.19× — fft Float64 1000 dims=(1,): 138.0 µs → 43.3 µs (FFTW 10.3 µs)
• 3.30× — fft Float64 21875 dims=(1,): 5.62 ms → 1.70 ms (FFTW 506.0 µs)
• 3.62× — fft Float64 25 dims=(1,): 2.0 µs → 0.6 µs (FFTW 0.1 µs)
• 3.66× — fft Float64 7 dims=(1,): 0.3 µs → 0.1 µs (FFTW 0.0 µs)

261 matched cases; geometric-mean speedup 1.94×; 0 cases slower by >5%.

(Real-transform rows to follow in a comment: the suite run mis-detected mul! for real plans on this branch; fixed in #128.)

Tests: new test/twiddles.jl pins the table layouts and accuracy, the direction check, the on-the-fly wrappers and zero allocations for planned execution; existing suite passes unchanged apart from the extended accuracy grid.

[codecov]

Codecov Report

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

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #134      +/-   ##
==========================================
- Coverage   98.80%   98.79%   -0.01%     
==========================================
  Files           5        5              
  Lines         585      666      +81     
==========================================
+ Hits          578      658      +80     
- 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.

[pankgeorg]

Real-transform rows for the table in the description (same machine and settings; "before" is main; the 10 errors are the unsupported 3D rfft, as on main):

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	rfft	Float32	29	2.34× (2.05–2.99)	7.73× → 3.34×	33 MiB → 9 MiB	0.4 µs → 1.9 ms
1d/awkward	rfft	Float64	29	2.21× (1.91–3.03)	6.17× → 2.81×	66 MiB → 18 MiB	0.4 µs → 2.2 ms
1d/pow2	rfft	Float32	20	1.43× (0.95–1.74)	4.47× → 3.21×	16 MiB → 16 MiB	0.1 µs → 369.0 µs
1d/pow2	rfft	Float64	20	1.42× (0.95–1.88)	3.29× → 2.39×	32 MiB → 32 MiB	0.1 µs → 423.6 µs
1d/prime	rfft	Float32	19	2.74× (2.13–3.47)	6.57× → 2.43×	116 MiB → 20 MiB	0.1 µs → 980.5 µs
1d/prime	rfft	Float64	19	2.52× (2.06–3.04)	6.25× → 2.52×	232 MiB → 40 MiB	0.1 µs → 1.1 ms
1d/smooth	rfft	Float32	24	1.99× (1.09–2.81)	10.47× → 5.36×	16 MiB → 16 MiB	0.8 µs → 506.2 µs
1d/smooth	rfft	Float64	24	2.19× (1.05–3.29)	8.66× → 4.16×	32 MiB → 32 MiB	0.8 µs → 571.6 µs
2d	rfft	Float32	9	1.44× (1.32–1.54)	13.23× → 9.28×	80 MiB → 80 MiB	0.2 µs → 162.9 µs
2d	rfft	Float64	14	1.81× (1.32–2.89)	11.83× → 6.68×	160 MiB → 160 MiB	0.2 µs → 180.6 µs
batched_dim1	rfft	Float32	6	1.47× (1.37–1.52)	6.39× → 4.33×	32 MiB → 32 MiB	0.1 µs → 193.2 µs
batched_dim1	rfft	Float64	9	1.49× (1.42–1.55)	4.90× → 3.29×	65 MiB → 65 MiB	0.1 µs → 189.1 µs
batched_dim2	rfft	Float32	6	1.38× (1.31–1.42)	3.54× → 2.61×	32 MiB → 32 MiB	0.1 µs → 193.7 µs
batched_dim2	rfft	Float64	6	1.37× (1.29–1.43)	3.05× → 2.28×	65 MiB → 65 MiB	0.1 µs → 199.6 µs

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

• 0.95× — rfft Float32 8 dims=(1,): 0.1 µs → 0.1 µs (FFTW 0.0 µs)
• 0.95× — rfft Float64 8 dims=(1,): 0.1 µs → 0.1 µs (FFTW 0.0 µs)
• 1.05× — rfft Float64 12 dims=(1,): 0.3 µs → 0.3 µs (FFTW 0.0 µs)
• 1.09× — rfft Float32 12 dims=(1,): 0.3 µs → 0.2 µs (FFTW 0.0 µs)
• 3.15× — rfft Float64 25 dims=(1,): 1.6 µs → 0.5 µs (FFTW 0.1 µs)
• 3.29× — rfft Float64 21875 dims=(1,): 5.11 ms → 1.56 ms (FFTW 206.9 µs)
• 3.33× — rfft Float32 17 dims=(1,): 0.7 µs → 0.2 µs (FFTW 0.1 µs)
• 3.33× — rfft Float32 71 dims=(1,): 10.6 µs → 3.2 µs (FFTW 1.4 µs)
• 3.38× — rfft Float32 61 dims=(1,): 7.9 µs → 2.3 µs (FFTW 1.0 µs)
• 3.47× — rfft Float32 41 dims=(1,): 3.7 µs → 1.1 µs (FFTW 0.5 µs)

234 matched cases; geometric-mean speedup 1.94×; 1 cases slower by >5%.

[wheeheee]

Made a few patches for easy performance wins :) wh_patch1.patch
There are still a few more I spotted that you could ask Claude to do (that mostly benefits unplanned FFTs) e.g.

1. inlined sincospi in unit_roots and twiddle
2. there are 7 extra stores when N is a multiple of 8, which could maybe be eliminated?
3. it might be more cache-friendly (for large N) to fill in an octant or two first (say, W[k + 1] and W[q + k + 1], then read from them and negate/swap appropriately to 2 other pairs, and so on