[uoj 34] 多项式乘法

阐述见快速傅立叶变换及其逆变换. C++自带的complex类比手写的慢不少, 但是我觉得每次都手写一堆重载有些无聊 QAQ

使用合并DFT的技巧, complex的常数基本被抵消了......所以以下代码的效率跟手写复数类+两次DFT和一次iDFT差不多.

如果只进行两三次FFT (包括iFFT), 预处理对效率的优化不是很明显.

松爷教导我们不要把数组大小开成2的幂.

#include <bits/stdc++.h>
#define Rep(i, a, b) for (int i = (a), i##_ = (b); i < i##_; ++i)
#define For(i, a, b) for (int i = (a), i##_ = (b); i <= i##_; ++i)

using namespace std;
const double pi = acos(-1);
const int MAX_N = 1<<18;
typedef complex<double> Comp;

namespace Trans {
	int inv[MAX_N + 1];
	void init(int n)
	{
		int s = 1;
		while ((1<<s) < n) ++s;
		--s;
		Rep (i, 1, n) inv[i] = (inv[i>>1]>>1) | ((i&1)<<s);
	}

	void fft(Comp a[], int n, int d=1)
	{
		Rep (i, 0, n) if (inv[i] < i) swap(a[inv[i]], a[i]);
		
		for (int m = 1; m < n; m *= 2) {
			Comp _w(cos(pi/m), d * sin(pi/m));
			for (int i = 0; i < n; i += 2*m) {
				Comp w(1);
				Rep (j, 0, m) {
					Comp t = w * a[i+m+j];
					a[i+m+j] = a[i+j] - t;
					a[i+j] += t;
					w *= _w;
				}
			}
		}
		
		if (d == -1)
			Rep (i, 0, n) a[i] = a[i].real() / n; // Comp(a[i].real()/n, a[i].imag()/n)
	}

	void convol(Comp a[], Comp b[], int n)
	{
		fft(a, n);
		Rep (i, 0, n) {
			int j = (n-i) & (n-1);
			b[i] = (conj(a[j]*a[j]) - a[i]*a[i]) * Comp(0, 0.25);
		}
		fft(b, n, -1);
	}
}

Comp a[MAX_N + 1], b[MAX_N + 1];

int main()
{
	int n, m, x, l = 1;
	scanf("%d%d", &n, &m);
	while (l < n+m+1) l *= 2;
	Trans::init(l);

	For (i, 0, n) scanf("%d", &x), a[i] = x;
	For (i, 0, m) scanf("%d", &x), a[i] += Comp(0, x);

	Trans::convol(a, b, l);

	For (i, 0, n+m) printf("%d%c", (int)floor(real(b[i]) + 0.5), " \n"[i == n+m]);
	return 0;
}