rustfft/algorithm/
good_thomas_algorithm.rs

1use std::cmp::max;
2use std::sync::Arc;
3
4use num_complex::Complex;
5use num_integer::Integer;
6use strength_reduce::StrengthReducedUsize;
7use transpose;
8
9use crate::array_utils;
10use crate::{common::FftNum, FftDirection};
11use crate::{Direction, Fft, Length};
12
13/// Implementation of the [Good-Thomas Algorithm (AKA Prime Factor Algorithm)](https://en.wikipedia.org/wiki/Prime-factor_FFT_algorithm)
14///
15/// This algorithm factors a size n FFT into n1 * n2, where GCD(n1, n2) == 1
16///
17/// Conceptually, this algorithm is very similar to the Mixed-Radix, except because GCD(n1, n2) == 1 we can do some
18/// number theory trickery to reduce the number of floating-point multiplications and additions. Additionally, It can
19/// be faster than Mixed-Radix at sizes below 10,000 or so.
20///
21/// ~~~
22/// // Computes a forward FFT of size 1200, using the Good-Thomas Algorithm
23/// use rustfft::algorithm::GoodThomasAlgorithm;
24/// use rustfft::{Fft, FftPlanner};
25/// use rustfft::num_complex::Complex;
26/// use rustfft::num_traits::Zero;
27///
28/// let mut buffer = vec![Complex{ re: 0.0f32, im: 0.0f32 }; 1200];
29///
30/// // we need to find an n1 and n2 such that n1 * n2 == 1200 and GCD(n1, n2) == 1
31/// // n1 = 48 and n2 = 25 satisfies this
32/// let mut planner = FftPlanner::new();
33/// let inner_fft_n1 = planner.plan_fft_forward(48);
34/// let inner_fft_n2 = planner.plan_fft_forward(25);
35///
36/// // the good-thomas FFT length will be inner_fft_n1.len() * inner_fft_n2.len() = 1200
37/// let fft = GoodThomasAlgorithm::new(inner_fft_n1, inner_fft_n2);
38/// fft.process(&mut buffer);
39/// ~~~
40pub struct GoodThomasAlgorithm<T> {
41    width: usize,
42    width_size_fft: Arc<dyn Fft<T>>,
43
44    height: usize,
45    height_size_fft: Arc<dyn Fft<T>>,
46
47    reduced_width: StrengthReducedUsize,
48    reduced_width_plus_one: StrengthReducedUsize,
49
50    inplace_scratch_len: usize,
51    outofplace_scratch_len: usize,
52    immut_scratch_len: usize,
53
54    len: usize,
55    direction: FftDirection,
56}
57
58impl<T: FftNum> GoodThomasAlgorithm<T> {
59    /// Creates a FFT instance which will process inputs/outputs of size `width_fft.len() * height_fft.len()`
60    ///
61    /// `GCD(width_fft.len(), height_fft.len())` must be equal to 1
62    pub fn new(mut width_fft: Arc<dyn Fft<T>>, mut height_fft: Arc<dyn Fft<T>>) -> Self {
63        assert_eq!(
64            width_fft.fft_direction(), height_fft.fft_direction(),
65            "width_fft and height_fft must have the same direction. got width direction={}, height direction={}",
66            width_fft.fft_direction(), height_fft.fft_direction());
67
68        let mut width = width_fft.len();
69        let mut height = height_fft.len();
70        let direction = width_fft.fft_direction();
71
72        // This algorithm doesn't work if width and height aren't coprime
73        let gcd = num_integer::gcd(width as i64, height as i64);
74        assert!(gcd == 1,
75                "Invalid width and height for Good-Thomas Algorithm (width={}, height={}): Inputs must be coprime",
76                width,
77                height);
78
79        // The trick we're using for our index remapping will only work if width < height, so just swap them if it isn't
80        if width > height {
81            std::mem::swap(&mut width, &mut height);
82            std::mem::swap(&mut width_fft, &mut height_fft);
83        }
84
85        let len = width * height;
86
87        // Collect some data about what kind of scratch space our inner FFTs need
88        let width_inplace_scratch = width_fft.get_inplace_scratch_len();
89        let height_inplace_scratch = height_fft.get_inplace_scratch_len();
90        let height_outofplace_scratch = height_fft.get_outofplace_scratch_len();
91
92        // Computing the scratch we'll require is a somewhat confusing process.
93        // When we compute an out-of-place FFT, both of our inner FFTs are in-place
94        // When we compute an inplace FFT, our inner width FFT will be inplace, and our height FFT will be out-of-place
95        // For the out-of-place FFT, one of 2 things can happen regarding scratch:
96        //      - If the required scratch of both FFTs is <= self.len(), then we can use the input or output buffer as scratch, and so we need 0 extra scratch
97        //      - If either of the inner FFTs require more, then we'll have to request an entire scratch buffer for the inner FFTs,
98        //          whose size is the max of the two inner FFTs' required scratch
99        let max_inner_inplace_scratch = max(height_inplace_scratch, width_inplace_scratch);
100        let outofplace_scratch_len = if max_inner_inplace_scratch > len {
101            max_inner_inplace_scratch
102        } else {
103            0
104        };
105
106        // For the in-place FFT, again the best case is that we can just bounce data around between internal buffers, and the only inplace scratch we need is self.len()
107        // If our height fft's OOP FFT requires any scratch, then we can tack that on the end of our own scratch, and use split_at_mut to separate our own from our internal FFT's
108        // Likewise, if our width inplace FFT requires more inplace scracth than self.len(), we can tack that on to the end of our own inplace scratch.
109        // Thus, the total inplace scratch is our own length plus the max of what the two inner FFTs will need
110        let inplace_scratch_len = len
111            + max(
112                if width_inplace_scratch > len {
113                    width_inplace_scratch
114                } else {
115                    0
116                },
117                height_outofplace_scratch,
118            );
119
120        let immut_scratch_len = max(
121            width_fft.get_inplace_scratch_len(),
122            len + height_fft.get_inplace_scratch_len(),
123        );
124
125        Self {
126            width,
127            width_size_fft: width_fft,
128
129            height,
130            height_size_fft: height_fft,
131
132            reduced_width: StrengthReducedUsize::new(width),
133            reduced_width_plus_one: StrengthReducedUsize::new(width + 1),
134
135            inplace_scratch_len,
136            outofplace_scratch_len,
137            immut_scratch_len,
138
139            len,
140            direction,
141        }
142    }
143
144    fn reindex_input(&self, source: &[Complex<T>], destination: &mut [Complex<T>]) {
145        // A critical part of the good-thomas algorithm is re-indexing the inputs and outputs.
146        // To remap the inputs, we will use the CRT mapping, paired with the normal transpose we'd do for mixed radix.
147        //
148        // The algorithm for the CRT mapping will work like this:
149        // 1: Keep an output index, initialized to 0
150        // 2: The output index will be incremented by width + 1
151        // 3: At the start of the row, compute if we will increment output_index past self.len()
152        //      3a: If we will, then compute exactly how many increments it will take,
153        //      3b: Increment however many times as we scan over the input row, copying each element to the output index
154        //      3c: Subtract self.len() from output_index
155        // 4: Scan over the rest of the row, incrementing output_index, and copying each element to output_index, thne incrementing output_index
156        // 5: The first index of each row will be the final index of the previous row plus one, but because of our incrementing (width+1) inside the loop, we overshot, so at the end of the row, subtract width from output_index
157        //
158        // This ends up producing the same result as computing the multiplicative inverse of width mod height and etc by the CRT mapping, but with only one integer division per row, instead of one per element.
159        let mut destination_index = 0;
160        for mut source_row in source.chunks_exact(self.width) {
161            let increments_until_cycle =
162                1 + (self.len() - destination_index) / self.reduced_width_plus_one;
163
164            // If we will have to rollover output_index on this row, do it in a separate loop
165            if increments_until_cycle < self.width {
166                let (pre_cycle_row, post_cycle_row) = source_row.split_at(increments_until_cycle);
167
168                for input_element in pre_cycle_row {
169                    destination[destination_index] = *input_element;
170                    destination_index += self.reduced_width_plus_one.get();
171                }
172
173                // Store the split slice back to input_row, os that outside the loop, we can finish the job of iterating the row
174                source_row = post_cycle_row;
175                destination_index -= self.len();
176            }
177
178            // Loop over the entire row (if we did not roll over) or what's left of the row (if we did) and keep incrementing output_row
179            for input_element in source_row {
180                destination[destination_index] = *input_element;
181                destination_index += self.reduced_width_plus_one.get();
182            }
183
184            // The first index of the next will be the final index this row, plus one.
185            // But because of our incrementing (width+1) inside the loop above, we overshot, so subtract width, and we'll get (width + 1) - width = 1
186            destination_index -= self.width;
187        }
188    }
189
190    fn reindex_output(&self, source: &[Complex<T>], destination: &mut [Complex<T>]) {
191        // A critical part of the good-thomas algorithm is re-indexing the inputs and outputs.
192        // To remap the outputs, we will use the ruritanian mapping, paired with the normal transpose we'd do for mixed radix.
193        //
194        // The algorithm for the ruritanian mapping will work like this:
195        // 1: At the start of every row, compute the output index = (y * self.height) % self.width
196        // 2: We will increment this output index by self.width for every element
197        // 3: Compute where in the row the output index will wrap around
198        // 4: Instead of starting copying from the beginning of the row, start copying from after the rollover point
199        // 5: When we hit the end of the row, continue from the beginning of the row, continuing to increment the output index by self.width
200        //
201        // This achieves the same result as the modular arithmetic ofthe ruritanian mapping, but with only one integer divison per row, instead of one per element
202        for (y, source_chunk) in source.chunks_exact(self.height).enumerate() {
203            let (quotient, remainder) =
204                StrengthReducedUsize::div_rem(y * self.height, self.reduced_width);
205
206            // Compute our base index and starting point in the row
207            let mut destination_index = remainder;
208            let start_x = self.height - quotient;
209
210            // Process the first part of the row
211            for x in start_x..self.height {
212                destination[destination_index] = source_chunk[x];
213                destination_index += self.width;
214            }
215
216            // Wrap back around to the beginning of the row and keep incrementing
217            for x in 0..start_x {
218                destination[destination_index] = source_chunk[x];
219                destination_index += self.width;
220            }
221        }
222    }
223
224    fn perform_fft_inplace(&self, buffer: &mut [Complex<T>], scratch: &mut [Complex<T>]) {
225        let (scratch, inner_scratch) = scratch.split_at_mut(self.len());
226
227        // Re-index the input, copying from the buffer to the scratch in the process
228        self.reindex_input(buffer, scratch);
229
230        // run FFTs of size `width`
231        let width_scratch = if inner_scratch.len() > buffer.len() {
232            &mut inner_scratch[..]
233        } else {
234            &mut buffer[..]
235        };
236        self.width_size_fft
237            .process_with_scratch(scratch, width_scratch);
238
239        // transpose
240        transpose::transpose(scratch, buffer, self.width, self.height);
241
242        // run FFTs of size 'height'
243        self.height_size_fft
244            .process_outofplace_with_scratch(buffer, scratch, inner_scratch);
245
246        // Re-index the output, copying from the scratch to the buffer in the process
247        self.reindex_output(scratch, buffer);
248    }
249
250    fn perform_fft_immut(
251        &self,
252        input: &[Complex<T>],
253        output: &mut [Complex<T>],
254        scratch: &mut [Complex<T>],
255    ) {
256        // Re-index the input, copying from the input to the output in the process
257        self.reindex_input(input, output);
258
259        // run FFTs of size `width`
260        self.width_size_fft.process_with_scratch(output, scratch);
261
262        let (scratch, inner_scratch) = scratch.split_at_mut(self.len());
263
264        // transpose
265        transpose::transpose(output, scratch, self.width, self.height);
266
267        // run FFTs of size 'height'
268        self.height_size_fft
269            .process_with_scratch(scratch, inner_scratch);
270
271        // Re-index the output, copying from the input to the output in the process
272        self.reindex_output(scratch, output);
273    }
274
275    fn perform_fft_out_of_place(
276        &self,
277        input: &mut [Complex<T>],
278        output: &mut [Complex<T>],
279        scratch: &mut [Complex<T>],
280    ) {
281        // Re-index the input, copying from the input to the output in the process
282        self.reindex_input(input, output);
283
284        // run FFTs of size `width`
285        let width_scratch = if scratch.len() > input.len() {
286            &mut scratch[..]
287        } else {
288            &mut input[..]
289        };
290        self.width_size_fft
291            .process_with_scratch(output, width_scratch);
292
293        // transpose
294        transpose::transpose(output, input, self.width, self.height);
295
296        // run FFTs of size 'height'
297        let height_scratch = if scratch.len() > output.len() {
298            &mut scratch[..]
299        } else {
300            &mut output[..]
301        };
302        self.height_size_fft
303            .process_with_scratch(input, height_scratch);
304
305        // Re-index the output, copying from the input to the output in the process
306        self.reindex_output(input, output);
307    }
308}
309boilerplate_fft!(
310    GoodThomasAlgorithm,
311    |this: &GoodThomasAlgorithm<_>| this.len,
312    |this: &GoodThomasAlgorithm<_>| this.inplace_scratch_len,
313    |this: &GoodThomasAlgorithm<_>| this.outofplace_scratch_len,
314    |this: &GoodThomasAlgorithm<_>| this.immut_scratch_len
315);
316
317/// Implementation of the Good-Thomas Algorithm, specialized for smaller input sizes
318///
319/// This algorithm factors a size n FFT into n1 * n2, where GCD(n1, n2) == 1
320///
321/// Conceptually, this algorithm is very similar to MixedRadix, except because GCD(n1, n2) == 1 we can do some
322/// number theory trickery to reduce the number of floating point operations. It typically performs
323/// better than MixedRadixSmall, especially at the smallest sizes.
324///
325/// ~~~
326/// // Computes a forward FFT of size 56 using GoodThomasAlgorithmSmall
327/// use std::sync::Arc;
328/// use rustfft::algorithm::GoodThomasAlgorithmSmall;
329/// use rustfft::algorithm::butterflies::{Butterfly7, Butterfly8};
330/// use rustfft::{Fft, FftDirection};
331/// use rustfft::num_complex::Complex;
332///
333/// let mut buffer = vec![Complex{ re: 0.0f32, im: 0.0f32 }; 56];
334///
335/// // we need to find an n1 and n2 such that n1 * n2 == 56 and GCD(n1, n2) == 1
336/// // n1 = 7 and n2 = 8 satisfies this
337/// let inner_fft_n1 = Arc::new(Butterfly7::new(FftDirection::Forward));
338/// let inner_fft_n2 = Arc::new(Butterfly8::new(FftDirection::Forward));
339///
340/// // the good-thomas FFT length will be inner_fft_n1.len() * inner_fft_n2.len() = 56
341/// let fft = GoodThomasAlgorithmSmall::new(inner_fft_n1, inner_fft_n2);
342/// fft.process(&mut buffer);
343/// ~~~
344pub struct GoodThomasAlgorithmSmall<T> {
345    width: usize,
346    width_size_fft: Arc<dyn Fft<T>>,
347
348    height: usize,
349    height_size_fft: Arc<dyn Fft<T>>,
350
351    input_output_map: Box<[usize]>,
352
353    direction: FftDirection,
354}
355
356impl<T: FftNum> GoodThomasAlgorithmSmall<T> {
357    /// Creates a FFT instance which will process inputs/outputs of size `width_fft.len() * height_fft.len()`
358    ///
359    /// `GCD(width_fft.len(), height_fft.len())` must be equal to 1
360    pub fn new(width_fft: Arc<dyn Fft<T>>, height_fft: Arc<dyn Fft<T>>) -> Self {
361        assert_eq!(
362            width_fft.fft_direction(), height_fft.fft_direction(),
363            "n1_fft and height_fft must have the same direction. got width direction={}, height direction={}",
364            width_fft.fft_direction(), height_fft.fft_direction());
365
366        let width = width_fft.len();
367        let height = height_fft.len();
368        let len = width * height;
369
370        assert_eq!(width_fft.get_outofplace_scratch_len(), 0, "GoodThomasAlgorithmSmall should only be used with algorithms that require 0 out-of-place scratch. Width FFT (len={}) requires {}, should require 0", width, width_fft.get_outofplace_scratch_len());
371        assert_eq!(height_fft.get_outofplace_scratch_len(), 0, "GoodThomasAlgorithmSmall should only be used with algorithms that require 0 out-of-place scratch. Height FFT (len={}) requires {}, should require 0", height, height_fft.get_outofplace_scratch_len());
372
373        assert!(width_fft.get_inplace_scratch_len() <= width, "GoodThomasAlgorithmSmall should only be used with algorithms that require little inplace scratch. Width FFT (len={}) requires {}, should require {} or less", width, width_fft.get_inplace_scratch_len(), width);
374        assert!(height_fft.get_inplace_scratch_len() <= height, "GoodThomasAlgorithmSmall should only be used with algorithms that require little inplace scratch. Height FFT (len={}) requires {}, should require {} or less", height, height_fft.get_inplace_scratch_len(), height);
375
376        // compute the multiplicative inverse of width mod height and vice versa. x will be width mod height, and y will be height mod width
377        let gcd_data = i64::extended_gcd(&(width as i64), &(height as i64));
378        assert!(gcd_data.gcd == 1,
379                "Invalid input width and height to Good-Thomas Algorithm: ({},{}): Inputs must be coprime",
380                width,
381                height);
382
383        // width_inverse or height_inverse might be negative, make it positive by wrapping
384        let width_inverse = if gcd_data.x >= 0 {
385            gcd_data.x
386        } else {
387            gcd_data.x + height as i64
388        } as usize;
389        let height_inverse = if gcd_data.y >= 0 {
390            gcd_data.y
391        } else {
392            gcd_data.y + width as i64
393        } as usize;
394
395        // NOTE: we are precomputing the input and output reordering indexes, because benchmarking shows that it's 10-20% faster
396        // If we wanted to optimize for memory use or setup time instead of multiple-FFT speed, we could compute these on the fly in the perform_fft() method
397        let input_iter = (0..len)
398            .map(|i| (i % width, i / width))
399            .map(|(x, y)| (x * height + y * width) % len);
400        let output_iter = (0..len).map(|i| (i % height, i / height)).map(|(y, x)| {
401            (x * height * height_inverse as usize + y * width * width_inverse as usize) % len
402        });
403
404        let input_output_map: Vec<usize> = input_iter.chain(output_iter).collect();
405
406        Self {
407            direction: width_fft.fft_direction(),
408
409            width,
410            width_size_fft: width_fft,
411
412            height,
413            height_size_fft: height_fft,
414
415            input_output_map: input_output_map.into_boxed_slice(),
416        }
417    }
418
419    fn perform_fft_immut(
420        &self,
421        input: &[Complex<T>],
422        output: &mut [Complex<T>],
423        scratch: &mut [Complex<T>],
424    ) {
425        // These asserts are for the unsafe blocks down below. we're relying on the optimizer to get rid of this assert
426        assert_eq!(self.len(), input.len());
427        assert_eq!(self.len(), output.len());
428
429        let (input_map, output_map) = self.input_output_map.split_at(self.len());
430
431        // copy the input using our reordering mapping
432        for (output_element, &input_index) in output.iter_mut().zip(input_map.iter()) {
433            *output_element = input[input_index];
434        }
435
436        // run FFTs of size `width`
437        self.width_size_fft.process_with_scratch(output, scratch);
438
439        // transpose
440        unsafe { array_utils::transpose_small(self.width, self.height, output, scratch) };
441
442        // run FFTs of size 'height'
443        self.height_size_fft.process_with_scratch(scratch, output);
444
445        // copy to the output, using our output redordeing mapping
446        for (input_element, &output_index) in scratch.iter().zip(output_map.iter()) {
447            output[output_index] = *input_element;
448        }
449    }
450
451    fn perform_fft_out_of_place(
452        &self,
453        input: &mut [Complex<T>],
454        output: &mut [Complex<T>],
455        _scratch: &mut [Complex<T>],
456    ) {
457        // These asserts are for the unsafe blocks down below. we're relying on the optimizer to get rid of this assert
458        assert_eq!(self.len(), input.len());
459        assert_eq!(self.len(), output.len());
460
461        let (input_map, output_map) = self.input_output_map.split_at(self.len());
462
463        // copy the input using our reordering mapping
464        for (output_element, &input_index) in output.iter_mut().zip(input_map.iter()) {
465            *output_element = input[input_index];
466        }
467
468        // run FFTs of size `width`
469        self.width_size_fft.process_with_scratch(output, input);
470
471        // transpose
472        unsafe { array_utils::transpose_small(self.width, self.height, output, input) };
473
474        // run FFTs of size 'height'
475        self.height_size_fft.process_with_scratch(input, output);
476
477        // copy to the output, using our output redordeing mapping
478        for (input_element, &output_index) in input.iter().zip(output_map.iter()) {
479            output[output_index] = *input_element;
480        }
481    }
482
483    fn perform_fft_inplace(&self, buffer: &mut [Complex<T>], scratch: &mut [Complex<T>]) {
484        // These asserts are for the unsafe blocks down below. we're relying on the optimizer to get rid of this assert
485        assert_eq!(self.len(), buffer.len());
486        assert_eq!(self.len(), scratch.len());
487
488        let (input_map, output_map) = self.input_output_map.split_at(self.len());
489
490        // copy the input using our reordering mapping
491        for (output_element, &input_index) in scratch.iter_mut().zip(input_map.iter()) {
492            *output_element = buffer[input_index];
493        }
494
495        // run FFTs of size `width`
496        self.width_size_fft.process_with_scratch(scratch, buffer);
497
498        // transpose
499        unsafe { array_utils::transpose_small(self.width, self.height, scratch, buffer) };
500
501        // run FFTs of size 'height'
502        self.height_size_fft
503            .process_outofplace_with_scratch(buffer, scratch, &mut []);
504
505        // copy to the output, using our output redordeing mapping
506        for (input_element, &output_index) in scratch.iter().zip(output_map.iter()) {
507            buffer[output_index] = *input_element;
508        }
509    }
510}
511boilerplate_fft!(
512    GoodThomasAlgorithmSmall,
513    |this: &GoodThomasAlgorithmSmall<_>| this.width * this.height,
514    |this: &GoodThomasAlgorithmSmall<_>| this.len(),
515    |_| 0,
516    |this: &GoodThomasAlgorithmSmall<_>| this.len()
517);
518
519#[cfg(test)]
520mod unit_tests {
521    use super::*;
522    use crate::test_utils::check_fft_algorithm;
523    use crate::{algorithm::Dft, test_utils::BigScratchAlgorithm};
524    use num_integer::gcd;
525    use num_traits::Zero;
526    use std::sync::Arc;
527
528    #[test]
529    fn test_good_thomas() {
530        for width in 1..12 {
531            for height in 1..12 {
532                if gcd(width, height) == 1 {
533                    test_good_thomas_with_lengths(width, height, FftDirection::Forward);
534                    test_good_thomas_with_lengths(width, height, FftDirection::Inverse);
535                }
536            }
537        }
538    }
539
540    #[test]
541    fn test_good_thomas_small() {
542        let butterfly_sizes = [2, 3, 4, 5, 6, 7, 8, 16];
543        for width in &butterfly_sizes {
544            for height in &butterfly_sizes {
545                if gcd(*width, *height) == 1 {
546                    test_good_thomas_small_with_lengths(*width, *height, FftDirection::Forward);
547                    test_good_thomas_small_with_lengths(*width, *height, FftDirection::Inverse);
548                }
549            }
550        }
551    }
552
553    fn test_good_thomas_with_lengths(width: usize, height: usize, direction: FftDirection) {
554        let width_fft = Arc::new(Dft::new(width, direction)) as Arc<dyn Fft<f32>>;
555        let height_fft = Arc::new(Dft::new(height, direction)) as Arc<dyn Fft<f32>>;
556
557        let fft = GoodThomasAlgorithm::new(width_fft, height_fft);
558
559        check_fft_algorithm(&fft, width * height, direction);
560    }
561
562    fn test_good_thomas_small_with_lengths(width: usize, height: usize, direction: FftDirection) {
563        let width_fft = Arc::new(Dft::new(width, direction)) as Arc<dyn Fft<f32>>;
564        let height_fft = Arc::new(Dft::new(height, direction)) as Arc<dyn Fft<f32>>;
565
566        let fft = GoodThomasAlgorithmSmall::new(width_fft, height_fft);
567
568        check_fft_algorithm(&fft, width * height, direction);
569    }
570
571    #[test]
572    fn test_output_mapping() {
573        let width = 15;
574        for height in 3..width {
575            if gcd(width, height) == 1 {
576                let width_fft =
577                    Arc::new(Dft::new(width, FftDirection::Forward)) as Arc<dyn Fft<f32>>;
578                let height_fft =
579                    Arc::new(Dft::new(height, FftDirection::Forward)) as Arc<dyn Fft<f32>>;
580
581                let fft = GoodThomasAlgorithm::new(width_fft, height_fft);
582
583                let mut buffer = vec![Complex { re: 0.0, im: 0.0 }; fft.len()];
584
585                fft.process(&mut buffer);
586            }
587        }
588    }
589
590    // Verify that the Good-Thomas algorithm correctly provides scratch space to inner FFTs
591    #[test]
592    fn test_good_thomas_inner_scratch() {
593        let scratch_lengths = [1, 5, 24];
594
595        let mut inner_ffts = Vec::new();
596
597        for &len in &scratch_lengths {
598            for &inplace_scratch in &scratch_lengths {
599                for &outofplace_scratch in &scratch_lengths {
600                    for &immut_scratch in &scratch_lengths {
601                        inner_ffts.push(Arc::new(BigScratchAlgorithm {
602                            len,
603                            inplace_scratch,
604                            outofplace_scratch,
605                            immut_scratch,
606                            direction: FftDirection::Forward,
607                        }) as Arc<dyn Fft<f32>>);
608                    }
609                }
610            }
611        }
612
613        for width_fft in inner_ffts.iter() {
614            for height_fft in inner_ffts.iter() {
615                if width_fft.len() == height_fft.len() {
616                    continue;
617                }
618
619                let fft = GoodThomasAlgorithm::new(Arc::clone(width_fft), Arc::clone(height_fft));
620
621                let mut inplace_buffer = vec![Complex::zero(); fft.len()];
622                let mut inplace_scratch = vec![Complex::zero(); fft.get_inplace_scratch_len()];
623
624                fft.process_with_scratch(&mut inplace_buffer, &mut inplace_scratch);
625
626                let mut outofplace_input = vec![Complex::zero(); fft.len()];
627                let mut outofplace_output = vec![Complex::zero(); fft.len()];
628                let mut outofplace_scratch =
629                    vec![Complex::zero(); fft.get_outofplace_scratch_len()];
630
631                fft.process_outofplace_with_scratch(
632                    &mut outofplace_input,
633                    &mut outofplace_output,
634                    &mut outofplace_scratch,
635                );
636
637                let immut_input = vec![Complex::zero(); fft.len()];
638                let mut immut_output = vec![Complex::zero(); fft.len()];
639                let mut immut_scratch = vec![Complex::zero(); fft.get_immutable_scratch_len()];
640
641                fft.process_immutable_with_scratch(
642                    &immut_input,
643                    &mut immut_output,
644                    &mut immut_scratch,
645                );
646            }
647        }
648    }
649}
rustfft/algorithm/good_thomas_algorithm.rs

rustfft/algorithm/
good_thomas_algorithm.rs
