首页 > 代码库 > FFT小结

FFT小结

2024-07-13 01:36:18   220人阅读

先上模板

 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其实是可以撑过1048576的，只要你的π精度足够高。

说什么FFT精度炸翔的，应该是这样子的：

const double PI=3.14159265359;

不炸翔才怪。调了一个上午，发现跟std比对后，第i个数的误差正比于sin(2*PI*i/N)，然后就在那边调常数，过了大点小点又Wa。其实= =不想say。

4. 这个模板的效率还是蛮高的，蝶形变换的时间比普通的省了不少。

 

可以找这题练习：Number Theory Transform


联系
我们