From d446c25b54c009d8c708e8acfc25c056f9a300c2 Mon Sep 17 00:00:00 2001 From: Henrik Date: Wed, 9 Sep 2026 07:58:11 +0200 Subject: [PATCH 1/2] Precompute the Rader's permutation 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. --- src/algorithm/raders_algorithm.rs | 99 ++++++++++++++++--------------- 1 file changed, 52 insertions(+), 47 deletions(-) diff --git a/src/algorithm/raders_algorithm.rs b/src/algorithm/raders_algorithm.rs index 5bc45e10..c3e2883d 100644 --- a/src/algorithm/raders_algorithm.rs +++ b/src/algorithm/raders_algorithm.rs @@ -42,10 +42,15 @@ pub struct RadersAlgorithm { inner_fft: Arc>, inner_fft_data: Box<[Complex]>, - primitive_root: u64, - primitive_root_inverse: u64, - - len: StrengthReducedU64, + // Buffer offsets for the input reordering: inner FFT input k reads buffer position + // permutation[k], which is g^(k+1) mod len, less one. The sequence depends only on the + // primitive root and the length, so building it once here keeps a serial chain of + // strength-reduced modular multiplies out of the hot loops. avx_raders.rs precomputes its + // output mapping for the same reason. This costs 4 * (len - 1) bytes, the same order as the + // twiddles avx_raders.rs and Bluestein's already store at this size. + permutation: Box<[u32]>, + + len: usize, inplace_scratch_len: usize, outofplace_scratch_len: usize, immut_scratch_len: usize, @@ -95,6 +100,15 @@ impl RadersAlgorithm { ((twiddle_input as u64 * primitive_root_inverse) % reduced_len) as usize; } + // Precompute the input reordering. Stored already offset by one so the hot loops are a + // plain indexed load with no arithmetic. + let mut permutation = Vec::with_capacity(inner_fft_len); + let mut input_index = 1u64; + for _ in 0..inner_fft_len { + input_index = (input_index * primitive_root) % reduced_len; + permutation.push(u32::try_from(input_index - 1).unwrap()); + } + let required_inner_scratch = inner_fft.get_inplace_scratch_len(); let extra_inner_scratch = if required_inner_scratch <= inner_fft_len { 0 @@ -112,10 +126,9 @@ impl RadersAlgorithm { inner_fft, inner_fft_data: inner_fft_input.into_boxed_slice(), - primitive_root, - primitive_root_inverse, + permutation: permutation.into_boxed_slice(), - len: reduced_len, + len, inplace_scratch_len, outofplace_scratch_len: extra_inner_scratch, immut_scratch_len, @@ -123,6 +136,31 @@ impl RadersAlgorithm { } } + /// Copies `src` into `dest`, applying the input reordering. + #[inline] + fn gather_input(&self, src: &[Complex], dest: &mut [Complex]) { + for (dest_element, &input_index) in dest.iter_mut().zip(self.permutation.iter()) { + *dest_element = src[input_index as usize]; + } + } + + /// Copies `src` into `dest`, conjugating and applying the output reordering. + /// + /// The output reordering walks `g^-1` where the input reordering walks `g`. Since + /// `g^(len - 1) == 1` we have `g^-k == g^(len - 1 - k)`, so the inverse sequence is the + /// forward table read backwards and rotated by one. Peeling the rotated element off keeps + /// the loop a plain zip over two slices instead of a chained iterator. + #[inline] + fn scatter_output(&self, src: &[Complex], dest: &mut [Complex]) { + let (&last_index, head) = self.permutation.split_last().unwrap(); + let (last_element, src_head) = src.split_last().unwrap(); + + for (src_element, &output_index) in src_head.iter().zip(head.iter().rev()) { + dest[output_index as usize] = src_element.conj(); + } + dest[last_index as usize] = last_element.conj(); + } + fn perform_fft_immut( &self, input: &[Complex], @@ -135,13 +173,7 @@ impl RadersAlgorithm { let (scratch, extra_scratch) = scratch.split_at_mut(self.len() - 1); // copy the input into the scratch space, reordering as we go - let mut input_index = 1; - for output_element in scratch.iter_mut() { - input_index = ((input_index as u64 * self.primitive_root) % self.len) as usize; - - let input_element = input[input_index - 1]; - *output_element = input_element; - } + self.gather_input(input, scratch); self.inner_fft.process_with_scratch(scratch, extra_scratch); @@ -163,12 +195,7 @@ impl RadersAlgorithm { self.inner_fft.process_with_scratch(scratch, extra_scratch); // copy the final values into the output, reordering as we go - let mut output_index = 1; - for scratch_element in scratch { - output_index = - ((output_index as u64 * self.primitive_root_inverse) % self.len) as usize; - output[output_index - 1] = scratch_element.conj(); - } + self.scatter_output(scratch, output); } fn perform_fft_out_of_place( @@ -182,13 +209,7 @@ impl RadersAlgorithm { let (input_first, input) = input.split_first_mut().unwrap(); // copy the input into the output, reordering as we go. also compute a sum of all elements - let mut input_index = 1; - for output_element in output.iter_mut() { - input_index = ((input_index as u64 * self.primitive_root) % self.len) as usize; - - let input_element = input[input_index - 1]; - *output_element = input_element; - } + self.gather_input(input, output); // perform the first of two inner FFTs let inner_scratch = if scratch.len() > 0 { @@ -225,12 +246,7 @@ impl RadersAlgorithm { self.inner_fft.process_with_scratch(input, inner_scratch); // copy the final values into the output, reordering as we go - let mut output_index = 1; - for input_element in input { - output_index = - ((output_index as u64 * self.primitive_root_inverse) % self.len) as usize; - output[output_index - 1] = input_element.conj(); - } + self.scatter_output(input, output); } fn perform_fft_inplace(&self, buffer: &mut [Complex], scratch: &mut [Complex]) { // The first output element is just the sum of all the input elements, and we need to store off the first input value @@ -240,13 +256,7 @@ impl RadersAlgorithm { let (scratch, extra_scratch) = scratch.split_at_mut(self.len() - 1); // copy the buffer into the scratch, reordering as we go. also compute a sum of all elements - let mut input_index = 1; - for scratch_element in scratch.iter_mut() { - input_index = ((input_index as u64 * self.primitive_root) % self.len) as usize; - - let buffer_element = buffer[input_index - 1]; - *scratch_element = buffer_element; - } + self.gather_input(buffer, scratch); // perform the first of two inner FFTs let inner_scratch = if extra_scratch.len() > 0 { @@ -274,17 +284,12 @@ impl RadersAlgorithm { self.inner_fft.process_with_scratch(scratch, inner_scratch); // copy the final values into the output, reordering as we go - let mut output_index = 1; - for scratch_element in scratch { - output_index = - ((output_index as u64 * self.primitive_root_inverse) % self.len) as usize; - buffer[output_index - 1] = scratch_element.conj(); - } + self.scatter_output(scratch, buffer); } } boilerplate_fft!( RadersAlgorithm, - |this: &RadersAlgorithm<_>| this.len.get() as usize, + |this: &RadersAlgorithm<_>| this.len, |this: &RadersAlgorithm<_>| this.inplace_scratch_len, |this: &RadersAlgorithm<_>| this.outofplace_scratch_len, |this: &RadersAlgorithm<_>| this.immut_scratch_len From 78396b798e0fd5f720a62793930170ab2c4558e8 Mon Sep 17 00:00:00 2001 From: Henrik Date: Mon, 14 Sep 2026 21:53:44 +0200 Subject: [PATCH 2/2] Halve the Rader's twiddle precompute via conjugate symmetry len is an odd prime here, so the multiplicative group mod len is cyclic of even order len - 1, and primitive_root_inverse^((len-1)/2) == -1 mod len. That means twiddle_input at position p + (len-1)/2 equals len - twiddle_input at position p, so the second half of the twiddle array is just the complex conjugate of the first half. Only the first half needs an actual compute_twiddle call and a step through the modular-multiply chain; the rest is filled in with a conjugate. Idea from https://github.com/ejmahler/RustFFT/pull/178#discussion_r3995689422 --- src/algorithm/raders_algorithm.rs | 15 ++++++++++++++- 1 file changed, 14 insertions(+), 1 deletion(-) diff --git a/src/algorithm/raders_algorithm.rs b/src/algorithm/raders_algorithm.rs index c3e2883d..b5bab24f 100644 --- a/src/algorithm/raders_algorithm.rs +++ b/src/algorithm/raders_algorithm.rs @@ -91,10 +91,23 @@ impl RadersAlgorithm { // precompute the coefficients to use inside the process method let inner_fft_scale = T::one() / T::from_usize(inner_fft_len).unwrap(); let mut inner_fft_input = vec![Complex::zero(); inner_fft_len]; + + // primitive_root() above only returns Some for len > 2, so len is an odd prime here and + // inner_fft_len = len - 1 is even. That makes the multiplicative group mod len cyclic of + // even order, and in a cyclic group of even order, g^(order/2) is the unique element of + // order 2, which mod a prime is -1. That holds for primitive_root_inverse just as much + // as for primitive_root, so the twiddle_input sequence e_p = primitive_root_inverse^p + // mod len satisfies e_{p + (len-1)/2} == len - e_p. Since compute_twiddle(len - x, ..) + // is the complex conjugate of compute_twiddle(x, ..), the second half of this array is + // just the conjugate of the first half: only the first half needs an actual + // compute_twiddle (trig) call and a step through the modular-multiply chain. Idea from + // https://github.com/ejmahler/RustFFT/pull/178#discussion_r3995689422 + let (first_half, second_half) = inner_fft_input.split_at_mut(inner_fft_len / 2); let mut twiddle_input = 1; - for input_cell in &mut inner_fft_input { + for (input_cell, conjugate_input_cell) in first_half.iter_mut().zip(second_half) { let twiddle = twiddles::compute_twiddle(twiddle_input, len, direction); *input_cell = twiddle * inner_fft_scale; + *conjugate_input_cell = input_cell.conj(); twiddle_input = ((twiddle_input as u64 * primitive_root_inverse) % reduced_len) as usize;