Fast Fourier Transform

fft.hpp

#include <complex>
#include <vector>
#include <cmath>

#define COMPLEXE(u) std::complex<double>(cos(u),-sin(u))

int bitwiseReverse(int num, int log2n)
{
	int ans = 0;
	while (log2n--) {
		ans <<= 1;
		ans |= (num & 0x1);
		num >>= 1;
	}
	return ans;
}

std::vector<std::complex<double>> bitwiseReverse(const std::vector<double> &a)
{
	const int N = a.size();
	const int LG2N = log2(N);
	std::vector<std::complex<double>> reversedA(N);
	for (int i = 0; i < N; i++) 
		reversedA[bitwiseReverse(i, LG2N)] = a[i];
	return reversedA;
}

std::vector<std::complex<double>> fft(const std::vector<double> &a)
{
	using std::complex;
	using std::vector;
	const int N = a.size();
	vector<complex<double>> y = bitwiseReverse(a);
	for (int i = 1; i < N; i *= 2) {
		complex<double> wn = COMPLEXE(M_PI/i);
		for (int j = 0; j < N; j += i*2) {
			complex<double> w = complex<double>(1);
			for (int k = 0; k < i; k++) {
				complex<double> t = w*y[j+k+i];
				complex<double> u = y[j+k];
				y[j+k] = u+t;
				y[j+k+i] = u-t;
				w = w * wn;
			}
		}
	}
	return y;
}