bzoj2194 快速傅立叶之二

题目链接

problem

给出两个长度为n的数列a，b。求一个数列c满足：$$c[k] = sumlimits_{i = k} ^ na[i]b[i - k]$$

(nle 10^5)

solution

长得和卷积很像，观察一下卷积的形式：(c[k]=sumlimits_{i=0}^ia[i]b[k-i])

所以先把b数组翻转过来。

然后所求的式子就变成了(c[k]=sumlimits_{i=k}^na[i]b[n+k-i])

将(n+k)看做整体，那么形式就和卷积比较像了,令(d[n+k]=c[k]=sumlimits_{i=k}^na[i]b[n + k - i])。因为a的(n+1)到(n+k)位置均为0，b的(0)到(k-1)位置均为0.

所以上式也可以写为(d[n+k]=sumlimits_{i=0}^{n+k}a[i]b[n+k-i])。

这样就可以借助FFT计算了。

注意最后的答案是问c[k]，也就是d[k+n]

code

/*
* @Author: wxyww
* @Date:   2020-01-22 08:04:15
* @Last Modified time: 2020-01-22 08:24:04
*/
#include<cstdio>
#include<iostream>
#include<cstdlib>
#include<cstring>
#include<cmath>
#include<algorithm>
#include<queue>
#include<vector>
#include<ctime>
using namespace std;
typedef long long ll;
const int N = 400010;
const double pi = 3.1415926535897931;
ll read() {
	ll x=0,f=1;char c=getchar();
	while(c<'0'||c>'9') {
		if(c=='-') f=-1;
		c=getchar();
	}
	while(c>='0'&&c<='9') {
		x=x*10+c-'0';
		c=getchar();
	}
	return x*f;
}
struct complex {
	double x,y;
}A[N],B[N];

complex operator * (const complex &A,const complex &B) {
	return (complex){A.x * B.x - A.y * B.y,A.x * B.y + A.y * B.x};
}
complex operator + (const complex &A,const complex &B) {
	return (complex){A.x + B.x,A.y + B.y};
}
complex operator - (const complex &A,const complex &B) {
	return (complex){A.x - B.x,A.y - B.y};
}
int rev[N];
void FFT(complex *a,int n,int xs) {
	for(int i = 0;i <= n;++i) if(rev[i] > i) swap(a[rev[i]],a[i]);

	for(int m = 2;m <= n;m <<= 1) {
		complex w1 = (complex){cos(2 * pi / m),xs * sin(2 * pi / m)};
		for(int i = 0;i < n;i += m) {
			complex w = (complex){1,0};
			for(int k = 0;k < (m >> 1);++k) {
				complex u = a[i + k];
				complex t = w * a[i + k + (m >> 1)];
				a[i + k] = u + t;
				a[i + k + (m >> 1)] = u - t;
				w = w * w1;
			}
		}
	}

}

int main() {
	int n = read();
	for(int i = 0;i < n;++i) {
		A[i].x = read();B[n - i - 1].x = read();
	}
	--n;
	int tot = 1;
	while(tot <= (n << 1)) tot <<= 1;

	for(int i = 0;i <= tot;++i) rev[i] = (rev[i >> 1] >> 1) | (i & 1 ? tot >> 1 : 0);


	FFT(A,tot,1);FFT(B,tot,1);

	for(int i = 0;i <= tot;++i) A[i] = A[i] * B[i];
	FFT(A,tot,-1);
	
	for(int i = 0;i <= n;++i) printf("%d
",(int)round(A[i + n].x / tot));



	return 0;
}