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 }
NTT
 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 }
FFT

 

注意几点:

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

posted @ 2014-07-20 11:47  n+e  阅读(497)  评论(0编辑  收藏  举报