2N点的实数傅里叶变换

2n-fft.c

#include <cstdio>
#include <cmath>
#include <complex>
#include <algorithm>
 
const int MaxL = 18, MaxN = 1 << MaxL;
typedef std::complex<double> complex_t;
complex_t eps[MaxL], eps_inv[MaxL], A[MaxN], B[MaxN];
 
void init_eps()
{
	double pi = acos(-1), angle = 2.0 * pi / (1 << MaxL);
	eps[MaxL - 1] = complex_t(cos(angle), sin(angle));
	for(int i = MaxL - 2; i >= 0; --i)
		eps[i] = eps[i + 1] * eps[i + 1];
	for(int i = 0; i != MaxL; ++i)
		eps_inv[i] = conj(eps[i]);
}
 
void transform(int n, complex_t *z, complex_t *w)
{
	for(int i = 0, j = 0; i != n; ++i)
	{
		if(i > j) std::swap(z[i], z[j]);
		for(int l = n >> 1; (j ^= l) < l; l >>= 1);
	}
 
	for(int i = 2; i <= n; i <<= 1, ++w)
	{
		int m = i >> 1;
		complex_t *e = new complex_t[m];
		e[0] = 1.0;
		for(int j = 1; j != m; ++j)
			e[j] = e[j - 1] * *w;
		for(int j = 0; j < n; j += i)
		{
			for(int p = 0; p != m; ++p)
			{
				complex_t x = z[j + m + p] * e[p];
				z[j + m + p] = z[j + p] - x;
				z[j + p] += x;
			}
		}
		delete[] e;
	}
}
 
int main()
{
	int n, m;
	std::scanf("%d %d", &n, &m);
	init_eps();
	for(int i = 0; i <= n; ++i)
	{
		double a;
		std::scanf("%lf", &a);
		A[i] = a;
	}
 
	for(int i = 0; i <= m; ++i)
	{
		double a;
		std::scanf("%lf", &a);
		A[i] += std::complex<double>{0, a};
	}
 
	int p = 1;
	while(p < n + m + 1) p <<= 1;
	transform(p, A, eps);
	for(int i = 1; i != p; ++i)
	{
		double x1 = A[i].real(), y1 = A[i].imag();
		double x2 = A[p - i].real(), y2 = A[p - i].imag();
		std::complex<double> a = { (x1 + x2) * 0.5, (y1 - y2) * 0.5 };
		std::complex<double> b = { (y1 + y2) * 0.5, -(x1 - x2) * 0.5 };
		B[i] = a * b;
	}
 
	B[0] = A[0].imag() * A[0].real();
 
	transform(p, B, eps_inv);
	for(int i = 0; i <= n + m; ++i)
		std::printf("%d ", int(B[i].real() / p + 0.5));
	return 0;
}