[bzoj2194] 快速傅立叶之二

a, b是长度为n的数列, 从0开始标号. 求c[k] = Σ a[i] b[i-k] (k≤i<n). (a[i], b[i]为不大于100的非负整数, n≤10^5) 令a'[i] = a[n-1-i], 则 c[k] = Σ a'[n-1-k-j] b[j] (0≤j<n-k-1).

#include <bits/stdc++.h>
#define rep(i, a, b) for (int i = (a), i##_end = (b); i < i##_end; ++i)
#define per(i, a, b) for (int i = (a), i##_end = (b); i >= i##_end; --i)

using namespace std;

typedef complex<double> Complex;
const int N = (1<<18) + 1;
const double pi = acos(-1);

namespace Convol {
	int n, rev[N];
	void init(int _n, int s)
	{
		n = _n;
		--s;
		rep (i, 1, n)
			rev[i] = (rev[i>>1] >> 1) | ((i & 1) << s);
	}

	void fft(Complex a[], int d=1)
	{
		rep (i, 0, n)
			if (rev[i] < i)
				swap(a[rev[i]], a[i]);
		for (int m = 1; m < n; m <<= 1) {
			Complex _w(cos(pi/m), d*sin(pi/m));
			for (int j = 0; j < n; j += m<<1) {
				Complex w(1);
				rep (k, 0, m) {
					Complex t = w * a[j+m+k];
					a[j+m+k] = a[j+k] - t;
					a[j+k] += t;
					w *= _w;
				}
			}
		}
		if (d == -1)
			rep (i, 0, n) a[i] /= n;
	}

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

Complex x[N], y[N];

int main()
{
	int n, m = 1, s = 0;
	scanf("%d", &n);
	while (m < 2*n-1) m <<= 1, ++s;
	Convol::init(m, s);
	rep (i, 0, n) {
		int a, b;
		scanf("%d%d", &a, &b);
		x[n-1-i] += a;
		x[i] += Complex(0, b);
	}
	Convol::convol(x, y);
	per (i, n-1, 0)
		printf("%d\n", (int)round(real(y[i])));
	return 0;
}