2194: 快速傅立叶之二
Time Limit: 10 Sec Memory Limit: 259 MB
Submit: 1273 Solved: 745Description
请计算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 4Sample Output
24
12
10
6
1HINT
Source
【分析】
这个卷积也很容易看出来
C[k]=sigma(a[i]*b[i-k])
变形,先讲b数组倒置:
C[k]=sigma(a[i]*b[n-i+k])
令D[k]=sigma(a[i]*b[k-i])
则C[k]=D[n+k]
D数组就是标准卷积了,用FFT加速。
递归:
![](https://i-blog.csdnimg.cn/blog_migrate/8f900a89c6347c561fdf2122f13be562.gif)
![](https://i-blog.csdnimg.cn/blog_migrate/961ddebeb323a10fe0623af514929fc1.gif)
1 #include<cstdio> 2 #include<cstdlib> 3 #include<cstring> 4 #include<iostream> 5 #include<algorithm> 6 #include<cmath> 7 using namespace std; 8 #define Maxn 100010*4 9 const double pi=3.141592653; 10 11 struct P 12 { 13 double x,y; 14 P() {x=y=0;} 15 P(double x,double y):x(x),y(y){} 16 friend P operator + (P x,P y) {return P(x.x+y.x,x.y+y.y);} 17 friend P operator - (P x,P y) {return P(x.x-y.x,x.y-y.y);} 18 friend P operator * (P x,P y) {return P(x.x*y.x-x.y*y.y,x.x*y.y+x.y*y.x);} 19 }a[Maxn],b[Maxn]; 20 21 void fft(P *s,int n,int f) 22 { 23 if(n==1) return; 24 P a0[n>>1],a1[n>>1]; 25 for(int i=0;i<=n;i+=2) a0[i>>1]=s[i],a1[i>>1]=s[i+1]; 26 fft(a0,n>>1,f);fft(a1,n>>1,f); 27 P wn(cos(2*pi/n),f*sin(2*pi/n)),w(1,0); 28 for(int i=0;i<(n>>1);i++,w=w*wn) s[i]=a0[i]+w*a1[i],s[i+(n>>1)]=a0[i]-w*a1[i]; 29 } 30 31 int main() 32 { 33 int n; 34 scanf("%d",&n);n--; 35 for(int i=0;i<=n;i++) scanf("%lf%lf",&a[i].x,&b[i].x); 36 for(int i=0;i<=n/2;i++) swap(b[i].x,b[n-i].x); 37 int nn=1; 38 while(nn<=2*n) nn<<=1; 39 fft(a,nn,1);fft(b,nn,1); 40 for(int i=0;i<=nn;i++) a[i]=a[i]*b[i]; 41 fft(a,nn,-1); 42 for(int i=n;i<=n+n;i++) printf("%d\n",(int)(a[i].x/nn+0.5)); 43 return 0; 44 }
2017-04-13 10:39:50