# [知识点]在街角的咖啡店邂逅多项式

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 快速数论变换

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<<1));
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;
}


#include<iostream>
#include<cstring>
#include<cstdio>
using namespace std;
#define pos(i,aa,bb) for(int i=(aa);i<=(bb);++i)
#define N 300000
#define int long long
const int mod=1004535809;
int g=3;
int l,len;
int jc[N],jc_inv[N];
int C[N],G[N],G_inv[N];
int qpow(int a,int k){
int res=1;
while(k){
if(k&1){
res=(res*a)%mod;
}
a=(a*a)%mod;
k>>=1;
}
return res%mod;
}
void getjc(){
jc[0]=1;jc_inv[0]=qpow(jc[0],mod-2);
pos(i,1,l){
jc[i]=(jc[i-1]*i)%mod;
jc_inv[i]=qpow(jc[i],mod-2);
}
}
int rev[N],tmp[N];
void getrev(int bit){
for(int i=0;i<(1<<bit);i++)
rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
}
void init(int n){
int bit;len=2;
for(bit=1;(1<<bit)<n*2-1;bit++) len<<=1;
getrev(bit);
}
void ntt(int *a,int dft){
pos(i,0,len-1){
if(i<rev[i]) swap(a[i],a[rev[i]]);
}
for(int step=1;step<len;step<<=1){
int wnk=qpow(g,(mod-1)/(step<<1));
if(dft==-1) wnk=qpow(wnk,mod-2);
for(int j=0;j<len;j+=step<<1){
int wn=1;
pos(k,j,j+step-1){
int x=a[k],y=(wn*a[k+step])%mod;
a[k]=x+y;a[k+step]=x-y;
if(a[k]>=mod) a[k]-=mod;
if(a[k+step]<0) a[k+step]+=mod;
wn=(wn*wnk)%mod;
}
}
}
if(dft==-1){
int ni=qpow(len,mod-2);
pos(i,0,len-1) a[i]=a[i]*ni%mod;
}
}
void get_inv(int *a,int *b,int n){
if(n==1){
b[0]=qpow(a[0],mod-2);
return;
}
get_inv(a,b,(n+1)>>1);
init(n);
pos(i,0,n-1) tmp[i]=a[i];
pos(i,n,len-1) tmp[i]=0;
ntt(b,1);ntt(tmp,1);
pos(i,0,len-1) b[i]=(b[i]*(2-b[i]*tmp[i]%mod+mod))%mod;
ntt(b,-1);
pos(i,n,len-1) b[i]=0;
}
signed main(){
scanf("%lld",&l);
getjc();
pos(i,1,l) C[i]=qpow(2,(i*(i-1)/2)%(mod-1))*jc_inv[i-1]%mod;
//C[0]=1;
pos(i,0,l) G[i]=qpow(2,(i*(i-1)/2)%(mod-1))*jc_inv[i]%mod;
get_inv(G,G_inv,l+1);
init(l+1);
ntt(C,1);ntt(G_inv,1);
pos(i,0,len-1) C[i]=(C[i]*G_inv[i])%mod;
ntt(C,-1);
printf("%lld",C[l]*jc[l-1]%mod);
return 0;
}


posted @ 2018-07-08 21:27 Hallmeow 阅读(...) 评论(...) 编辑 收藏