快速傅立叶之二
Description
请计算C[k]=sigma(a[i]*b[i-k]) 其中 k < = i < n ,并且有 n < = 10 ^ 5。 a,b中的元素均为小于等于100的非负整数。
Input
第一行一个整数N,接下来N行,第i+2..i+N-1行,每行两个数,依次表示a[i],b[i] (0 < = i < N)。
Output
输出N行,每行一个整数,第i行输出C[i-1]。
Sample Input
5
3 1
2 4
1 1
2 4
1 4
Sample Output
24
12
10
6
1
Code
#include <bits/stdc++.h>
using namespace std;
template<typename T> inline void read(T &x){
x=0;T f=1;char ch=getchar();
while(!isdigit(ch)){if(ch=='-')f=-1;ch=getchar();}
while(isdigit(ch)){x=x*10+ch-'0';ch=getchar();}
x*=f;
}
const double pi=acos(-1.);
const int N=300100;
struct Comp{
double x,y;
Comp(double _x=0.,double _y=0.):x(_x),y(_y){}
}A[N],B[N];
inline Comp operator +(const Comp &a,const Comp &b){return Comp(a.x+b.x,a.y+b.y);}
inline Comp operator -(const Comp &a,const Comp &b){return Comp(a.x-b.x,a.y-b.y);}
inline Comp operator *(const Comp &a,const Comp &b){return Comp(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
int n,_n,m,R[N],L;
void FFT(Comp *a, int f){
for(int i=0;i<n;i++)if(i<R[i])swap(a[i],a[R[i]]);
for(int i=1;i<n;i<<=1){
Comp wn=Comp(cos(pi/i),f*sin(pi/i));
for(int j=0,p=i<<1;j<n;j+=p){
Comp w=Comp(1.,0.);
for(int k=0;k<i;w=w*wn,k++){
Comp A0=a[j+k],A1=w*a[j+k+i];
a[j+k]=A0+A1;a[j+k+i]=A0-A1;
}
}
}
}
int main(){
read(n); _n=n;
for(int i=0,x,y;i<n;i++){
read(x);read(y);
A[n-i-1].x=x; B[i].x=y;
}
for(m=2*n-2,n=1;n<=m;n<<=1,L++);
for(int i=0;i<n;i++)R[i]=(R[i>>1]>>1)|((i&1)<<(L-1));
FFT(A,1);FFT(B,1);for(int i=0;i<n;i++)A[i]=A[i]*B[i];FFT(A,-1);
for(int i=0;i<_n;i++) printf("%d\n", (int)round(A[_n-1-i].x/n));
return 0;
}