本来不打算入坑。离NOI也没几天了。OD说这个考的很多,让我看看,最起码基本的做做。
没办法硬着头皮看吧。当初愣是让NEYC的GGN讲了三遍FFT还是一头雾水。就当个半黑箱用吧。
这几篇写的很好:
用来求卷积
FFT模板:
cd a[N],b[N]; int rev[N]; void getrev(int bit){ for(int i=0;i<(1<<bit);i++){ rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1)); } } void fft(cd *a,int n,int dft){ pos(i,0,n){ if(i<rev[i]) swap(a[i],a[rev[i]]); } for(int step=1;step<n;step<<=1){ cd wn=exp(cd(0,dft*pi/step)); for(int j=0;j<n;j+=step<<1){ cd wnk(1.0,0.0); pos(k,j,j+step-1){ cd x=a[k],y=wnk*a[k+step]; a[k]=x+y;a[k+step]=x-y; wnk*=wn; } } } if(dft==-1) pos(i,0,n-1) a[i]/=n; }
多项式相乘就是卷积的一种
[BZOJ 2194] 快速傅立叶之二 裸题
#include<iostream> #include<cstdio> #include<cstring> #include<cmath> #include<complex> using namespace std; #define pos(i,a,b) for(int i=(a);i<=(b);i++) #define pos2(i,a,b) for(int i=(a);i>=(b);i--) #define N 1001000 const double pi=3.1415926535897932; typedef complex<double> cd; cd a[N],b[N]; int rev[N]; int len; int a1[N],b1[N]; void getrev(int bit){ for(int i=0;i<(1<<bit);i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1)); } void fft(cd *a,int n,int dft){ pos(i,0,n) if(i<rev[i]) swap(a[i],a[rev[i]]); for(int step=1;step<n;step<<=1){ cd wn=exp(cd(0,dft*pi/step)); for(int j=0;j<n;j+=step<<1){ cd wnk(1.0,0.0); pos(k,j,j+step-1){ cd x=a[k],y=wnk*a[k+step]; a[k]=x+y;a[k+step]=x-y; wnk*=wn; } } } if(dft==-1) pos(i,0,n-1) a[i]/=n; } int c[N]; int main(){ scanf("%d",&len); int bit=1,s=2; for(bit=1;(1<<bit)<len*2-1;bit++) s<<=1; getrev(bit); pos(i,0,len-1){ scanf("%d%d",&a1[len-1-i],&b1[i]); } pos(i,0,len-1) a[i]=(double)a1[i],b[i]=(double)b1[i]; fft(a,s,1);fft(b,s,1); pos(i,0,s-1) a[i]*=b[i]; fft(a,s,-1); pos2(i,len-1,0){ printf("%d\n",(int)(a[i].real()+0.5)); } return 0; }
NTT 快速数论变换
用来解决精度问题 在模意义下
不错的文章快速数论变换(NTT)
(先贴上一份不知道对不对的代码)
void getrev(int bit){ for(int i=0;i<(1<<bit);i++){ rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1)); } } int qpow(int aa,int k){ int res=1; while(k){ if(k&1) res=(res*aa)%p; aa=(aa*aa)%p; k>>=1; } return res; } int getni(int x){ return qpow(x,p-2); } void getg(){ //求p的原根 对p-1质因数分解 //若恒有 g^(p-1/pi) mod p !=1 则g为原根 } void ntt(int *a,int n,int dft){ pos(i,0,n){ if(i<rev[i]) swap(a[i],a[rev[i]]); } for(int step=1;step<n;step<<=1){ int wn=qpow(g,(p-1)/step); if(dft==-1) wn=getni(wn); for(int j=0;j<n;j+=step<<1){ int wnk=1; pos(k,j,j+step-1){ int x=a[k],y=wnk*a[k+step]%p; a[k]=x+y;a[k+step]=x-y; if(a[k]>=p) a[k]-=p; if(a[k+step]<0) a[k+step]+=p; wnk=(wnk*wn)%p; } } } int nin=getni(n); if(dft==-1) pos(i,0,n-1) a[i]=a[i]*nin%p; }
原文:https://www.cnblogs.com/Hallmeow/p/9281340.html