[BZOJ2194] 快速傅里叶之二 题解

题意：求$$C_k=sum_{k}^{n-1}a_ib_{i-k}.$$
(n leq 1e5).

考虑反转数组(a)，生成新数组(a').
那么$$C_k=sum_{i=k}^{n-1}a'{n-1-i}b{i-k},$$
考虑把(i)改成从(0)开始。那么

[C_k=sum_{i=0}^{n-k-1}a'_{n-k-1-i}b_i. ]

考虑用( ext{FFT})计算卷积的标准形式：

[C'_x=sum_{i=0}^{x}A_{x-i}B_i. ]

考虑(x=n-k-1)时的情况：

[C'_{n-k-1}=sum_{i=0}^{n-k-1}a'_{n-k-1-i}b_i. ]

发现(C')本质上就是将(C)的前(n)个反转了一下，而(C')是可以直接计算的。
于是，可以直接用( ext{FFT})计算出(a')和(b)的卷积(C')，再反转一下前(n)项，输出前(n)项即可。

#include <bits/stdc++.h>
#define REP(i , x , y) for(__typeof(y) i = x; i <= y; i++)
#define PER(i , y , x) for(__typeof(x) i = y; i >= x; i--)
typedef long long LL;
typedef long double ld;
typedef unsigned long long ULL;
/* do not give up ! try your best! Read the meaning clearly! */
template < typename T > void Input(T& t) {
	char c = getchar(); T x = 1; t = 0; while(!isdigit(c)) {if(c == '-') x = -1; c = getchar();}
	while(isdigit(c)) t = t * 10 + c - '0' , c = getchar();t *= x;
}
template < typename T , typename... Args > void Input(T& t , Args&... args) {Input(t); Input(args...);}
template < typename T > T mul(T x , T y , T _) {x %= _,y %= _;return ((x * y - (T)(((ld)x * y + 0.5) / _) * _) % _ + _) % _;}

using namespace std;

const int MAXN = 2097152 + 10;
const double Pi = acos(-1.0);
int a[MAXN] , b[MAXN] , c[MAXN] , n;
static complex < double > A[MAXN] , B[MAXN] , C[MAXN];
struct FastFourierTransform {
	complex < double > omega[MAXN] , omegaInverse[MAXN];
	void init(int n) {
		for(int i = 0; i < n; i++) {
			omega[i] = complex < double > (cos(2 * Pi / n * i) , sin(2 * Pi / n * i));
			omegaInverse[i] = conj(omega[i]);
		}
	}
	void Transform(complex < double > *a , const int n , const complex < double > *omega) {
		int k = 1;
		while((1 << k) < n) ++k;
		for(int i = 0; i < n; i++) {
			int t = 0;
			for(int j = 0;j < k; j++)
				if(i & (1 << j)) t |= (1 << (k - 1 - j));
			if(i < t) swap(a[i] , a[t]);
		}

		for(int l = 2; l <= n; l <<= 1) {
			int m = (l >> 1);
			for(complex < double > *p = a; p != a + n; p += l) {
				for(int i = 0; i < m; i++) {
					complex < double > k = omega[n / l * i] * p[m + i];
					p[m + i] = p[i] - k;
					p[i] += k;
				}
			}
		}
	}
	void DFT(complex < double > *a , int n) {
		Transform(a , n , omega);
	}
	void IDFT(complex < double > *a , int n) {
		Transform(a , n , omegaInverse);
		for(int i = 0; i < n; i++) a[i] /= n;
	}
}FFT;
int main() {
	Input(n);
	for(int i = 0; i < n; i++) Input(a[i] , b[i]);
	reverse(a , a + n);
	for(int i = 0; i < n; i++) A[i].real(a[i]) , B[i].real(b[i]);
	int N = 1;
	while(N < 2 * n) N <<= 1;
	FFT.init(N);
	FFT.DFT(A , N); FFT.DFT(B , N);
	for(int i = 0; i < N; i++) C[i] = A[i] * B[i];
	FFT.IDFT(C , N);
	for(int i = 0; i < N; i++) c[i] = static_cast < int > (C[i].real() + 0.5);
	int cnt = 0;
	for(int i = 0; i < n; i++) printf("%d
" , c[n - 1 - i]);
	return 0;
}