Precompute the Rader's permutation - #178
Open
HEnquist wants to merge 1 commit into
Open
Conversation
Rader's reorders the input by walking g^k mod len, one strength-reduced modular multiply per element, in a serial chain recomputed on every call. Each step is two 64x64->128 widening multiplies. wasm has no high-multiply opcode, so LLVM emits two `call __multi3` plus a linear-memory round trip per element there, where aarch64 does it in seven instructions. That, and not the scattered access pattern, is where the backends diverge: NEON has no gather either. The sequence depends only on the primitive root and the length, so build it once at construction into a Box<[u32]> and leave the hot loops as plain indexed loads. avx_raders.rs already precomputes its output mapping for the same reason. The inverse walk is that same table read backwards and rotated by one, since g^-k == g^(len - 1 - k), so one table covers both directions. Measured on an M1, f64, over eight smooth primes. Rader's runs its inner FFT twice, so subtracting two inner FFTs isolates the permutation pass: len raders ns -> raders ns speedup permutation ns/element 1009 17027 10544 1.61x 8.09 -> 1.62 2053 39685 26789 1.48x 8.10 -> 1.72 4051 116000 90233 1.29x 8.33 -> 1.90 8101 177284 126558 1.40x 8.32 -> 1.99 15121 426733 337206 1.27x 8.17 -> 2.13 32401 803375 615661 1.30x 8.12 -> 2.31 65537 1239810 867489 1.43x 8.08 -> 2.36 100801 2988425 2416661 1.24x 8.09 -> 2.27 Under wasm, where the chain was costing 25.2-26.9 ns per element, it drops to 1.87-3.07 and Rader's comes out 1.9-3.5x faster (node 24, f64). The table costs 4 * (len - 1) bytes, about 400 KB for a prime near 100000, the same order as the twiddles avx_raders.rs and Bluestein's already store at that size.
Owner
This is what was stopping me before, I was never able to figure out how to reuse the same buffer for both reorders. This is great work, thank you. Perhaps we should reconsider the raders vs bluestein's tuning with this in mind. Not necessarily as a part of this PR. I'll review next week. |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
2 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Rader's reorders the input by walking g^k mod len, one strength-reduced modular multiply per
element, in a serial chain recomputed on every call. This builds the sequence once at construction
instead.
two
call __multi3plus a linear-memory round trip per element there, where aarch64 does it inseven instructions
Box<[u32]>covers both directions: the inverse walk is that same table read backwards androtated by one, since
g^-k == g^(len - 1 - k)avx_raders.rsand Bluestein's already storepermutation cost drops from 25.2-26.9 ns to 1.87-3.07 and Rader's comes out 1.9-3.5x faster