#ifndef CA_UWATERLOO_ALUMNI_DWHARDER_FFT_ARRAY
#define CA_UWATERLOO_ALUMNI_DWHARDER_FFT_ARRAY

#include <cassert>
#include <iostream>
#include <cmath>
#include <complex>

// Author:  Douglas Wilhelm Harder
// Copyright (c) 2010 by Douglas Wilhelm Harder.  All rights reserved.

namespace Data_structures {
	/****************************************************
	 * ************************************************ *
	 * *          The Fast Fourier Transform          * *
	 * ************************************************ *
	 ****************************************************/

	class Array {
		public:
			/*
			 *  This is an implementation of the FFT which runs in 
			 *  O(n ln(n)) time with O(1) memory.  It is meant to
			 *  demonstrate an iterative implementation of the FFT.
			 *
			 *  Reference:  Maple 8, Waterloo Maple Inc., Waterloo, Ontario.
			 */

			static void FFT( std::complex<double> *array, int n ) {
				// Assert that n is a power of 2
				assert( n > 0 && (n & (~n + 1)) == n );

				// atan(1) == pi/4
				// Multiplication by 4 adds 2 to the exponent
				double const PI = 4.0*std::atan( 1.0 );

				for ( int i = n/2; i >= 1; i /= 2 ) {
					std::complex<double> w = 1.0;
					// The conjugate principal nth root of unity
					std::complex<double> wn =
						std::exp( std::complex<double>( 0.0, -PI/i ) );

					for ( int j = 0; j < i; ++j ) {
						for ( int k = j; k < n; k += 2*i ) {
							int m = k + i;

							std::complex<double> tmp = array[k] + array[m];
							array[m] = (array[k] - array[m])*w;
							array[k] = tmp;
						}

						w *= wn;
					}
				}

				for ( int i = 0, j = 0; i < (n - 1); ++i ) {
					if ( i < j ) {
						std::complex<double> tmp = array[j];
						array[j] = array[i];
						array[i] = tmp;
					}

					int k;

					for ( k = n/2; k <= j; k /= 2 ) {
						j = j - k;
					}

					j += k;
				}
			}

			/*
			 *  The inverse Fourier transform may be calculated using the
			 *  Fourier transform.
			 */

			static void iFFT( std::complex<double> *array, int n ) {
				for ( int i = 0; i < n; ++i ) {
					array[i] = std::conj( array[i] );
				}

				FFT( array, n );

				for ( int i = 0; i < n; ++i ) {
					array[i] = std::conj( array[i] )/double( n );
				}
			}
	};
}

#endif