FFT小结

时间:2021-05-24 20:45:01

先上模板

 #include<cstdio>
#include<cmath>
const int P=(<<)*+;
typedef long long ll;
ll power(ll t,int k,int mod){ll f=;for(;k;k>>=,t=t*t%mod)if(k&)f=f*t%mod;return f;}
int a[<<],b[<<],n,m,k,w[][<<],f;
void FFT(int x[],int k,int v)
{
int i,j,l,tmp;
for(i=j=;i<k;i++)
{
if(i>j)tmp=x[i],x[i]=x[j],x[j]=tmp;
for(l=k>>;(j^=l)<l;l>>=);
}
for(i=;i<=k;i<<=)
for(j=;j<k;j+=i)
for(l=;l<i>>;l++)
{
tmp=1LL*x[j+l+(i>>)]*w[v][k/i*l]%P;
x[j+l+(i>>)]=(1LL*x[j+l]-tmp+P)%P;
x[j+l]=(1LL*x[j+l]+tmp)%P;
}
}
int main(){
scanf("%d",&n);
for(i=;i<n;i++)scanf("%d%d",&a[i],&b[i]);
for(k=;k<n<<;k<<=);
w[][]=w[][k]=;j=power(,(P-)/k,P);
for(i=;i<k;i++)w[][i]=1LL*w[][i-]*j%P;
for(i=;i<=k;i++)w[][i]=w[][k-i];
FFT(a,k,);FFT(b,k,);
for(i=;i<k;i++)a[i]=1LL*a[i]*b[i]%P;
FFT(a,k,);j=power(k,P-,P);
for(i=;i<*n-;i++)printf("%d\n",1LL*a[i]*j%P);
}

NTT

 #include<cstdio>
#include<cmath>
typedef long long ll;
typedef double ld;
const ld PI=*asin();
struct P{ld x,y;};
P operator+(const P&a,const P&b){return (P){a.x+b.x,a.y+b.y};}
P operator-(const P&a,const P&b){return (P){a.x-b.x,a.y-b.y};}
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};} int a[<<],b[<<],n,m,k,f;ll c[<<];
P w[][<<],x[<<],y[<<];
void FFT(P*x,int k,int v)
{
int i,j,l;P tmp;
for(i=j=;i<k;i++)
{
if(i>j)tmp=x[i],x[i]=x[j],x[j]=tmp;
for(l=k>>;(j^=l)<l;l>>=);
}
for(i=;i<=k;i<<=)
for(j=;j<k;j+=i)
for(l=;l<i>>;l++)
{
tmp=x[j+l+(i>>)]*w[v][k/i*l];
x[j+l+(i>>)]=x[j+l]-tmp;
x[j+l]=x[j+l]+tmp;
}
}
int main(){
scanf("%d",&n);
for(i=;i<n;i++)scanf("%d%d",&a[i],&b[i]);
for(k=;k<n<<;k<<=);
for(i=;i<=k;i++)w[][k-i]=w[][i]=(P){cos(PI**i/k),sin(PI**i/k)};
for(i=;i<k;i++)x[i]=(P){a[i],};FFT(x,k,);
for(i=;i<k;i++)y[i]=(P){b[i],};FFT(y,k,);
for(i=;i<k;i++)x[i]=x[i]*y[i];FFT(x,k,); for(i=;i<*n-;i++)c[i]=(ll)(x[i].x/k+0.5);
for(i=;i<*n-;i++)printf("%lld\n",c[i]);
}

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