Bluestein: 3-smooth padding, embedded call graph, cutoff 47 - #136
Bluestein: 3-smooth padding, embedded call graph, cutoff 47#136pankgeorg wants to merge 7 commits into
Conversation
Twiddle factors were regenerated on every execution: every kernel seeded Singleton's recurrence with a sincospi call (per output row of the O(n^2) DFT leaf, per j1 in the composite step, per 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. CallGraph now carries, per node, a twiddle table in the layout its kernel reads sequentially (DFT: w^k; composite: the (j1, k2) block; radix-4/3: per-level interleaved triplets/pairs addressed by a flat offset), a BluesteinScratch (chirp, its pre-scaled transform, work arrays, pow2 tables for the padded length) per Bluestein node, and the direction the tables were built for. All tables derive from one unit_roots table per node that evaluates sincospi on the first octant only when 8 | N. Planned execution is allocation-free for every size. Tables are correctly rounded, so the Float32 error no longer grows with n (~1.5 ulp at 2^22 instead of ~1000); the accuracy test grid is extended to 2^22. The old kernel signatures remain as wrappers that build tables on the fly.
The padded convolution length is now the smallest 2^a 3^b >= 2N-1 when that is enough smaller than the next power of two to be cheaper (a length with factors of 3 costs about 1.6x per n log n compared to a power of two in FFTA), e.g. 4099 -> 8748 instead of 16384 and 65537 -> 139968 instead of 262144. The padded transforms run through a CallGraph embedded in the BluesteinScratch instead of the bare power-of-two kernel. With twiddle tables the O(N^2) DFT leaf and Bluestein cross over near N = 45 (measured for ComplexF64), so DEFAULT_BLUESTEIN_CUTOFF drops from 73 to 47.
… powers of two below 2048 Measured against the power-of-two alternative, the 3-smooth pad was slower for small transforms (73 -> 162 took 6.7 us instead of 5.1 us with 256): the composite step's fixed overhead dominates there, and the 1.6x per-element factor was optimistic; 1.9x fits the measured pairs (162/256, 8748/16384, 17496/32768, 139968/262144).
|
Follow-up pushed (071712b): the before/after table above was generated with the first cost factor (1.6), and a direct check of its regressions showed the two large ones (49757, 120779) to be host noise (E is 7% and 5% faster there when measured back-to-back), but the small-prime ones real: 73 and 79 were padded to 162 = 2·81, which costs 6.7 µs against 5.1 µs for 256 — the composite step's fixed overhead dominates at that size. The factor is now 1.9 (fits the four measured pairs 162/256, 8748/16384, 17496/32768, 139968/262144) and 3-smooth pads are only considered from 2048 up. 73/79 are back to 5.0 µs; 8443 and 65537 keep their 1.5×/1.3× gains. I'll refresh the table in the description once the re-run completes. |
|
Caveat worth stating in the PR itself: the two constants in |
Codecov Report✅ All modified and coverable lines are covered by tests. Additional details and impacted files@@ Coverage Diff @@
## main #136 +/- ##
==========================================
+ Coverage 98.80% 98.82% +0.02%
==========================================
Files 5 5
Lines 585 681 +96
==========================================
+ Hits 578 673 +95
- Misses 7 8 +1 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
…measurements on two architectures On x86-64 (AVX2) a 3-smooth length costs 2.0-3.1x per n log n relative to a power of two, against 1.3-2.3x on aarch64 (NEON); with the previous 1.9 threshold the 3-smooth pad was a 1.58x win at N = 8443 on aarch64 and a 1.19x loss on x86-64. A threshold of 2.1 never loses on either machine. The DFT/Bluestein crossover is n = 23 on x86-64 and about 29 on aarch64, so DEFAULT_BLUESTEIN_CUTOFF drops from 47 to 29.
…currently never selects a 3-smooth length
1 participant
Item E of the plan in #130, stacked on #134 (twiddle tables; the diff includes that commit until it lands).
Cutoff. With tables the O(N²) DFT leaf and Bluestein cross over at N ≈ 23 on x86-64 (AVX2, Core Ultra 7 165H: Bluestein is 1.7–2.2× faster than the DFT leaf for N = 41–47) and at N ≈ 29 on aarch64 (Neoverse-N1, for the 64-point pad; the DFT leaf is briefly up to 1.45× faster again at 37–43 when the pad grows to 128, and Bluestein is ahead from 47).
DEFAULT_BLUESTEIN_CUTOFFgoes from 73 to 29, which never loses much on either machine. Both datasets are in the source comment.Padding. Bluestein's algorithm computes a length-
Ntransform as a convolution of length≥ 2N-1. FFTA pads to the next power of two; FFTW pads to a smooth length (1009 → 2025, 4099 → 8640). This PR addsbluestein_pad_length(N; factor), which takes the smallest2^a·3^b ≥ 2N-1whenfactor·m·log m < p·log pfor the next power of twop, and runs the padded transforms through aCallGraphembedded in theBluesteinScratch(so they benefit from every kernel improvement) instead of the bare power-of-two kernel. The factor turned out to be architecture dependent: measured pern log n, a 3-smooth length costs 1.3–2.3× a power of two on aarch64 but 2.0–3.1× on x86-64, and with the first value tried (1.9) the 3-smooth pad was a 1.58× win at N = 8443 on aarch64 and a 1.19× loss on x86-64 (measured by an independent run on an AVX2 machine; details below). A general-purpose library cannot know its host, so the default is 2.1, at which the chooser never loses on either machine — and, since the admissible length is always more than half the power of two, currently never fires. The mechanism stays, documented with both measurements, for when the composite/radix-3 path becomes competitive with the radix-4 kernel; the practical content of this PR is the cutoff and the embedded graph.Rader's algorithm (which FFTW uses for 61, 73, 65537, …) is not part of this PR; it would bring primes with smooth
N-1down further and is noted inbenchmark/ANALYSIS.mdas follow-up work.x86-64 calibration data (Intel Core Ultra 7 165H, AVX2, Julia 1.12.6, ComplexF64, planned execution, min of samples; same method as the aarch64 numbers): cost per
n log nof a 3-smooth length relative to a power of two — 288 vs 256: 2.74× (aarch64 2.26×), 9216 vs 16384: 2.01× (1.28×), 8748 vs 16384: 3.11× (1.45×). End-to-end Bluestein with the 3-smooth pad vs the power of two: 4099: 1.19× slower, 8443: 1.19× slower (aarch64: 1.58× faster), 65537: 1.01× slower (aarch64: 1.17× faster). DFT leaf vs Bluestein: 19: 0.71/0.82 µs, 23: 1.19/0.91, 31: 2.19/0.90, 41: 3.88/2.29, 47: 5.13/2.31, 61: 8.70/2.34.Before/after (aarch64 Neoverse-N1, Julia 1.12.6,
benchmark/suite.jl --only 1d --kinds fft, planned execution, single thread,ComplexF64; "before" is #134):Primes and prime×small sizes below 10^5,
ComplexF64, quiet host (Float32rows follow the same pattern). Noise reference: the pow2/smooth rows, which this PR does not touch, come out at 1.02× (range 0.71–1.26) in the same run. The calibration data behind the constants is inbenchmark/x86-64/calib_bluestein_x86.md,calib_cutoff_x86.mdandcalib_padchoice_x86.mdon the branch of #128.Tests:
bluestein_pad_lengthis pinned for representative sizes, the Bluestein scratch of a composite with a prime factor is checked, and the allocation-count bound of the one-shotfft(x)regression test is raised (primes ≥ 47 now build their scratch at plan time).