FFT小结
先上模板
1 #include<cstdio>
2 #include<cmath>
3 const int P=(1<<21)*479+1;
4 typedef long long ll;
5 ll power(ll t,int k,int mod){ll f=1;for(;k;k>>=1,t=t*t%mod)if(k&1)f=f*t%mod;return f;}
6 int a[1<<20],b[1<<20],n,m,k,w[2][1<<20],f;
7 void FFT(int x[],int k,int v)
8 {
9 int i,j,l,tmp;
10 for(i=j=0;i<k;i++)
11 {
12 if(i>j)tmp=x[i],x[i]=x[j],x[j]=tmp;
13 for(l=k>>1;(j^=l)<l;l>>=1);
14 }
15 for(i=2;i<=k;i<<=1)
16 for(j=0;j<k;j+=i)
17 for(l=0;l<i>>1;l++)
18 {
19 tmp=1LL*x[j+l+(i>>1)]*w[v][k/i*l]%P;
20 x[j+l+(i>>1)]=(1LL*x[j+l]-tmp+P)%P;
21 x[j+l]=(1LL*x[j+l]+tmp)%P;
22 }
23 }
24 int main(){
25 scanf("%d",&n);
26 for(i=0;i<n;i++)scanf("%d%d",&a[i],&b[i]);
27 for(k=1;k<n<<1;k<<=1);
28 w[0][0]=w[0][k]=1;j=power(3,(P-1)/k,P);
29 for(i=1;i<k;i++)w[0][i]=1LL*w[0][i-1]*j%P;
30 for(i=0;i<=k;i++)w[1][i]=w[0][k-i];
31 FFT(a,k,0);FFT(b,k,0);
32 for(i=0;i<k;i++)a[i]=1LL*a[i]*b[i]%P;
33 FFT(a,k,1);j=power(k,P-2,P);
34 for(i=0;i<2*n-1;i++)printf("%d\n",1LL*a[i]*j%P);
35 }
1 #include<cstdio>
2 #include<cmath>
3 typedef long long ll;
4 typedef double ld;
5 const ld PI=2*asin(1);
6 struct P{ld x,y;};
7 P operator+(const P&a,const P&b){return (P){a.x+b.x,a.y+b.y};}
8 P operator-(const P&a,const P&b){return (P){a.x-b.x,a.y-b.y};}
9 P operator*(const P&a,const P&b){double d=a.x*b.x,e=a.y*b.y,f=(a.x+a.y)*(b.x+b.y);return (P){d-e,f-e-d};}
10
11 int a[1<<20],b[1<<20],n,m,k,f;ll c[1<<20];
12 P w[2][1<<20],x[1<<20],y[1<<20];
13 void FFT(P*x,int k,int v)
14 {
15 int i,j,l;P tmp;
16 for(i=j=0;i<k;i++)
17 {
18 if(i>j)tmp=x[i],x[i]=x[j],x[j]=tmp;
19 for(l=k>>1;(j^=l)<l;l>>=1);
20 }
21 for(i=2;i<=k;i<<=1)
22 for(j=0;j<k;j+=i)
23 for(l=0;l<i>>1;l++)
24 {
25 tmp=x[j+l+(i>>1)]*w[v][k/i*l];
26 x[j+l+(i>>1)]=x[j+l]-tmp;
27 x[j+l]=x[j+l]+tmp;
28 }
29 }
30 int main(){
31 scanf("%d",&n);
32 for(i=0;i<n;i++)scanf("%d%d",&a[i],&b[i]);
33 for(k=1;k<n<<1;k<<=1);
34 for(i=0;i<=k;i++)w[1][k-i]=w[0][i]=(P){cos(PI*2*i/k),sin(PI*2*i/k)};
35 for(i=0;i<k;i++)x[i]=(P){a[i],0};FFT(x,k,0);
36 for(i=0;i<k;i++)y[i]=(P){b[i],0};FFT(y,k,0);
37 for(i=0;i<k;i++)x[i]=x[i]*y[i];FFT(x,k,1);
38
39 for(i=0;i<2*n-1;i++)c[i]=(ll)(x[i].x/k+0.5);
40 for(i=0;i<2*n-1;i++)printf("%lld\n",c[i]);
41 }
注意几点:
1. 理论上有c=8,实际算了下,c大概在80左右,还是NTT,FFT就更高了。
2. NTT中注意乘爆的地方,一定要加1LL*,否则呵呵
3. FFT其实是可以撑过1048575的,只要你的PI精度足够高并且被乘数<32767。亲自测试,不服来辩。
说什么FFT精度炸翔的,应该是这样子的:
const double PI=3.14159265359;
不炸翔才怪。调了一个上午,发现跟std比对后,第i个数的误差正比于sin(2*PI*i/N),然后就在那边调常数,过了大点小点又Wa。其实= =不想say。
4. 这个模板的效率还是蛮高的,蝶形变换的时间比普通的省了不少。
可以找这题练习:Fast Number Theoretic Transform