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

本来不打算入坑。离NOI也没几天了。OD说这个考的很多,让我看看,最起码基本的做做。

没办法硬着头皮看吧。当初愣是让NEYC的GGN讲了三遍FFT还是一头雾水。就当个半黑箱用吧。

这几篇写的很好:

有关多项式的算法

从多项式乘法到快速傅里叶变换

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

  

多项式求逆

此文章非常棒 多项式求逆元

打了道[BZOJ 3456]城市规划 看题解感觉出现了新思路一样

#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 阅读(...) 评论(...) 编辑 收藏