本来不打算入坑。离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