#ifndef CA_UWATERLOO_ALUMNI_DWHARDER_FFT_BASIC
#define CA_UWATERLOO_ALUMNI_DWHARDER_FFT_BASIC

#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 Basic {
		public:
			/*
			 *  This is a simple implementation of a fast Fourier transform
			 *  to demonstrate that the algorithm runs in O(n ln(n)) time.
			 *
			 *  For clarity, new memory is allocated for the two arrays;
			 *  however, it is possible to implement this algorithm with
			 *  only O(1) additional memory.
			 */
		
			static void FFT( std::complex<double> *array, int n ) {
				// Assert that n is a power of 2
				assert( n > 0 && (n & (~n + 1)) == n );
		
				if ( n == 1 ) {
					return;
				}
		
				std::complex<double> even[n/2];
				std::complex<double>  odd[n/2];
		
				for ( int i = 0; i < n/2; ++i ) {
					even[i] = array[2*i];
					 odd[i] = array[2*i + 1];
				}
		
				FFT( even, n/2 );
				FFT(  odd, n/2 );
		
				// atan(1) == pi/4
				// Multiplication by 4 adds 2 to the exponent
				double const PI = 4.0*std::atan( 1.0 );

				std::complex<double> w = 1.0;
				// The conjugate principal nth root of unity
				std::complex<double> wn =
					std::exp( std::complex<double>( 0.0, -2.0*PI/n ) );
		
				for ( int i = 0; i < n/2; ++i ) {
					array[i]       = even[i] + w*odd[i];
					array[i + n/2] = even[i] - w*odd[i];
					w = w * wn;
				}
			}
		
			/*
			 *  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