【BZOJ】3527: [Zjoi2014]力

FFT卷积算法解析

最新推荐文章于 2018-03-09 19:09:25 发布

weixin_30568715

最新推荐文章于 2018-03-09 19:09:25 发布

阅读量55

点赞数

CC 4.0 BY-SA版权

文章标签： php

原文链接：http://www.cnblogs.com/Dragon-Light/p/6491619.html

本文介绍了一种基于快速傅立叶变换（FFT）的卷积算法实现，通过具体代码示例详细展示了如何处理输入数据，并利用FFT进行高效计算。特别针对特定问题场景，通过消去$f_i$并换元来简化计算过程。

题目链接：http://www.lydsy.com/JudgeOnline/problem.php?id=3527

把${f_i}$消去之后换元就是卷积的形式，直接算就可以了。

 1 #include<iostream>
 2 #include<cstdio>
 3 #include<algorithm>
 4 #include<vector>
 5 #include<cstdlib>
 6 #include<cmath>
 7 #include<complex>
 8 #include<cstring>
 9 using namespace std;
10 #define maxn 1001000
11 #define llg long long 
12 #define Pi acos(-1.0)
13 #define yyj(a) freopen(a".in","r",stdin),freopen(a".out","w",stdout);
14 typedef complex<double> XU;
15 
16 struct FFT
17 {
18     llg rev[maxn],L,n;
19     double c[maxn];
20     void fft(XU *z,llg f)
21     {
22         for (llg i=0;i<n;i++) if (i<rev[i]) swap(z[i],z[rev[i]]);
23         for (llg i=1;i<n;i*=2)
24         {
25             XU wn(cos(Pi/i),sin(Pi/i));
26             for (llg p=i*2,j=0;j<n;j+=p)
27             {
28                 XU w(1,0);
29                 for (llg k=0;k<i;k++,w*=wn)
30                 {
31                     XU x=z[j+k],y=z[j+k+i]*w;
32                     z[j+k]=x+y,z[j+k+i]=x-y;
33                 } 
34             }
35         }
36         if (f==-1) reverse(z+1,z+n);
37     }
38 
39     XU a[maxn],b[maxn];
40     
41     void work()
42     {
43         llg m=(n-1)*3;
44         for (n=1;n<=m;n*=2) L++;
45         for (llg i=0;i<n;i++) rev[i]=(rev[i/2]/2)|((i&1)<<(L-1));
46         fft(a,1),fft(b,1);
47         for (llg i=0;i<n;i++) a[i]*=b[i];
48         fft(a,-1);
49         for (llg i=0;i<n;i++) c[i]=(a[i].real()/n);
50     }
51 }F1,F2;
52 
53 double x[maxn];
54 
55 int main()
56 {
57     yyj("li");
58     llg N;
59     cin>>N;
60     F1.n=N;
61     for (llg i=0;i<N;i++) scanf("%lf",&x[i]),F1.a[i]=x[i];
62     for (llg i=0;i<=N-2;i++) F1.b[i]=(double)-1.0/(double)(N-1-i)/(double)(N-1-i);
63     for (llg i=N;i<=2*N-2;i++) F1.b[i]=(double)1.0/(double)(N-1-i)/(double)(N-1-i);
64     F1.work();
65     for (llg i=N-1;i<2*N-1;i++) printf("%.9lf\n",F1.c[i]) ;
66     return 0;
67 }