use std::sync::Arc;
use num_complex::Complex;
use strength_reduce::StrengthReducedUsize;
use transpose;
use common::{FFTnum, verify_length, verify_length_divisible};
use math_utils;
use array_utils;
use ::{Length, IsInverse, FFT};
use algorithm::butterflies::FFTButterfly;
pub struct GoodThomasAlgorithm<T> {
width: usize,
width_size_fft: Arc<FFT<T>>,
height: usize,
height_size_fft: Arc<FFT<T>>,
input_x_stride: usize,
input_y_stride: usize,
len: StrengthReducedUsize,
inverse: bool,
}
impl<T: FFTnum> GoodThomasAlgorithm<T> {
pub fn new(width_fft: Arc<FFT<T>>, height_fft: Arc<FFT<T>>) -> Self {
assert_eq!(
width_fft.is_inverse(), height_fft.is_inverse(),
"width_fft and height_fft must both be inverse, or neither. got width inverse={}, height inverse={}",
width_fft.is_inverse(), height_fft.is_inverse());
let width = width_fft.len();
let height = height_fft.len();
let is_inverse = width_fft.is_inverse();
let (gcd, mut width_inverse, mut height_inverse) =
math_utils::extended_euclidean_algorithm(width as i64, height as i64);
assert!(gcd == 1,
"Invalid input width and height to Good-Thomas Algorithm: ({},{}): Inputs must be coprime",
width,
height);
if width_inverse < 0 {
width_inverse += height as i64;
}
if height_inverse < 0 {
height_inverse += width as i64;
}
Self {
width: width,
width_size_fft: width_fft,
height: height,
height_size_fft: height_fft,
input_x_stride: height_inverse as usize * height,
input_y_stride: width_inverse as usize * width,
len: StrengthReducedUsize::new(width * height),
inverse: is_inverse,
}
}
fn perform_fft(&self, input: &mut [Complex<T>], output: &mut [Complex<T>]) {
for (y, row) in output.chunks_mut(self.width).enumerate() {
let input_base = y * self.input_y_stride;
for (x, output_cell) in row.iter_mut().enumerate() {
let input_index = (input_base + x * self.input_x_stride) % self.len;
*output_cell = input[input_index];
}
}
self.width_size_fft.process_multi(output, input);
transpose::transpose(input, output, self.width, self.height);
self.height_size_fft.process_multi(output, input);
for (x, row) in input.chunks(self.height).enumerate() {
let output_base = x * self.height;
for (y, input_cell) in row.iter().enumerate() {
let output_index = (output_base + y * self.width) % self.len;
output[output_index] = *input_cell;
}
}
}
}
impl<T: FFTnum> FFT<T> for GoodThomasAlgorithm<T> {
fn process(&self, input: &mut [Complex<T>], output: &mut [Complex<T>]) {
verify_length(input, output, self.len());
self.perform_fft(input, output);
}
fn process_multi(&self, input: &mut [Complex<T>], output: &mut [Complex<T>]) {
verify_length_divisible(input, output, self.len());
for (in_chunk, out_chunk) in input.chunks_mut(self.len()).zip(output.chunks_mut(self.len())) {
self.perform_fft(in_chunk, out_chunk);
}
}
}
impl<T> Length for GoodThomasAlgorithm<T> {
#[inline(always)]
fn len(&self) -> usize {
self.width * self.height
}
}
impl<T> IsInverse for GoodThomasAlgorithm<T> {
#[inline(always)]
fn is_inverse(&self) -> bool {
self.inverse
}
}
pub struct GoodThomasAlgorithmDoubleButterfly<T> {
width: usize,
width_size_fft: Arc<FFTButterfly<T>>,
height: usize,
height_size_fft: Arc<FFTButterfly<T>>,
input_output_map: Box<[usize]>,
inverse: bool,
}
impl<T: FFTnum> GoodThomasAlgorithmDoubleButterfly<T> {
pub fn new(width_fft: Arc<FFTButterfly<T>>, height_fft: Arc<FFTButterfly<T>>) -> Self {
assert_eq!(
width_fft.is_inverse(), height_fft.is_inverse(),
"n1_fft and height_fft must both be inverse, or neither. got width inverse={}, height inverse={}",
width_fft.is_inverse(), height_fft.is_inverse());
let width = width_fft.len();
let height = height_fft.len();
let len = width * height;
let (gcd, mut width_inverse, mut height_inverse) =
math_utils::extended_euclidean_algorithm(width as i64, height as i64);
assert!(gcd == 1,
"Invalid input n1 and height to Good-Thomas Algorithm: ({},{}): Inputs must be coprime",
width,
height);
if width_inverse < 0 {
width_inverse += height as i64;
}
if height_inverse < 0 {
height_inverse += width as i64;
}
let input_iter = (0..len)
.map(|i| (i % width, i / width))
.map(|(x, y)| (x * height + y * width) % len);
let output_iter = (0..len)
.map(|i| (i % height, i / height))
.map(|(y, x)| (x * height * height_inverse as usize + y * width * width_inverse as usize) % len);
let input_output_map: Vec<usize> = input_iter.chain(output_iter).collect();
GoodThomasAlgorithmDoubleButterfly {
inverse: width_fft.is_inverse(),
width: width,
width_size_fft: width_fft,
height: height,
height_size_fft: height_fft,
input_output_map: input_output_map.into_boxed_slice(),
}
}
unsafe fn perform_fft(&self, input: &mut [Complex<T>], output: &mut [Complex<T>]) {
let (input_map, output_map) = self.input_output_map.split_at(self.len());
for (output_element, &input_index) in output.iter_mut().zip(input_map.iter()) {
*output_element = input[input_index];
}
self.width_size_fft.process_multi_inplace(output);
array_utils::transpose_small(self.width, self.height, output, input);
self.height_size_fft.process_multi_inplace(input);
for (input_element, &output_index) in input.iter().zip(output_map.iter()) {
output[output_index] = *input_element;
}
}
}
impl<T: FFTnum> FFT<T> for GoodThomasAlgorithmDoubleButterfly<T> {
fn process(&self, input: &mut [Complex<T>], output: &mut [Complex<T>]) {
verify_length(input, output, self.len());
unsafe { self.perform_fft(input, output) };
}
fn process_multi(&self, input: &mut [Complex<T>], output: &mut [Complex<T>]) {
verify_length_divisible(input, output, self.len());
for (in_chunk, out_chunk) in input.chunks_mut(self.len()).zip(output.chunks_mut(self.len())) {
unsafe { self.perform_fft(in_chunk, out_chunk) };
}
}
}
impl<T> Length for GoodThomasAlgorithmDoubleButterfly<T> {
#[inline(always)]
fn len(&self) -> usize {
self.width * self.height
}
}
impl<T> IsInverse for GoodThomasAlgorithmDoubleButterfly<T> {
#[inline(always)]
fn is_inverse(&self) -> bool {
self.inverse
}
}
#[cfg(test)]
mod unit_tests {
use super::*;
use std::sync::Arc;
use test_utils::{check_fft_algorithm, make_butterfly};
use algorithm::DFT;
use num_integer::gcd;
#[test]
fn test_good_thomas() {
for width in 1..12 {
for height in 1..12 {
if gcd(width, height) == 1 {
test_good_thomas_with_lengths(width, height, false);
test_good_thomas_with_lengths(width, height, true);
}
}
}
}
#[test]
fn test_good_thomas_double_butterfly() {
let butterfly_sizes = [2,3,4,5,6,7,8,16];
for width in &butterfly_sizes {
for height in &butterfly_sizes {
if gcd(*width, *height) == 1 {
test_good_thomas_butterfly_with_lengths(*width, *height, false);
test_good_thomas_butterfly_with_lengths(*width, *height, true);
}
}
}
}
fn test_good_thomas_with_lengths(width: usize, height: usize, inverse: bool) {
let width_fft = Arc::new(DFT::new(width, inverse)) as Arc<FFT<f32>>;
let height_fft = Arc::new(DFT::new(height, inverse)) as Arc<FFT<f32>>;
let fft = GoodThomasAlgorithm::new(width_fft, height_fft);
check_fft_algorithm(&fft, width * height, inverse);
}
fn test_good_thomas_butterfly_with_lengths(width: usize, height: usize, inverse: bool) {
let width_fft = make_butterfly(width, inverse);
let height_fft = make_butterfly(height, inverse);
let fft = GoodThomasAlgorithmDoubleButterfly::new(width_fft, height_fft);
check_fft_algorithm(&fft, width * height, inverse);
}
}