extern crate strided;
extern crate num;

use std::num::{Int, Float};
use num::complex::{Complex, Complex64};
use strided::{MutStrided, Strided};

// naive decimation-in-time radix-2 Cooley-Tukey fast Fourier
// transform.
fn fft(input: Strided<Complex64>, mut output: MutStrided<Complex64>, forward: bool) {
    /*
    assert!(input.len() == output.len() && input.len().count_ones() == 1); // is a power of two

    // base case.
    if input.len() == 1 {
        output[0] = input[0];
        return
    }

    let (evens, odds) = input.substrides2();
    let (mut start, mut end) = output.split_at(input.len() / 2);

    // recursively perform FFTs, writing the results into the first
    // and second half of the output array respectively.
    fft(evens, start.reborrow(), forward);
    fft(odds, end.reborrow(), forward);

    let two_pi: f64 = Float::two_pi();
    let expn = two_pi / input.len() as f64;
    // exp(±2π/N)
    let twiddle = Complex::from_polar(&1.0, &if forward {-expn} else {expn});

    let mut factor = Complex::one();

    // combine the subFFTs with the relations:
    //   X_k = E_k + exp(-2πk/N) * O_k
    //   X_{k+N/2} = E_k - exp(-2πk/N) * O_k
    // (or exp(2πk/N) for a forward FFT.)
    for (even, odd) in start.iter_mut().zip(end.iter_mut()) {
        let twiddled = factor * *odd;
        let e = *even;

        *even = e + twiddled;
        *odd = e - twiddled;
        factor = factor * twiddle;
    }
    */
}

fn main() {
    let (one, zero) = (Complex::one(), Complex::zero());
    let mut a = [one + one, one, one + one, one];
    let mut b = [zero, zero, zero, zero];

    fft(Strided::new(&a), MutStrided::new(&mut b), true);
    println!("{} -> {}", a.as_slice(), b.as_slice());

    fft(Strided::new(&b), MutStrided::new(&mut a), false);
    println!("{} -> {}", b.as_slice(), a.as_slice());
}