首页 > 代码库 > BZOJ 4332 FFT+快速幂

BZOJ 4332 FFT+快速幂

思路:

最裸的方程:f[i][j]=Σf[i-1][j-k]*F[k]

诶呦 这不是卷积嘛 

f[i]就可以用f[i-1]卷F 求到

但是这样还是很慢 

设p[i] 为Σ f[j](1<=j<=i)

发现p可以倍增推

于是  就 倍增一下  就完了...

http://www.cnblogs.com/Skyminer/p/6561689.html

hz神犇的题解写得非常详细..

 

 

//By SiriusRen#include <cmath>#include <cstdio>#include <algorithm>using namespace std;const int N=32888;const double pi=acos(-1);int n,M,P,A,O,S,U,L,R[N],F[N],g[N],p[N],t[N],ans;struct cplxd{double x,y;cplxd(){}cplxd(double X,double Y){x=X,y=Y;}}ca[N],cb[N],cc[N];cplxd operator+(cplxd a,cplxd b){return cplxd(a.x+b.x,a.y+b.y);}cplxd operator-(cplxd a,cplxd b){return cplxd(a.x-b.x,a.y-b.y);}cplxd operator*(cplxd a,cplxd b){return cplxd(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}cplxd operator/(cplxd a,int b){return cplxd(a.x/b,a.y/b);}void FFT(cplxd *a,int f){    for(int i=0;i<n;i++)if(i<R[i])swap(a[i],a[R[i]]);    for(int i=1;i<n;i<<=1){        cplxd wn=cplxd(cos(pi/i),f*sin(pi/i));        for(int j=0;j<n;j+=(i<<1)){            cplxd w(1,0);            for(int k=0;k<i;k++,w=w*wn){                cplxd x=a[j+k],y=w*a[j+k+i];                a[j+k]=x+y,a[j+k+i]=x-y;            }        }    }    if(f==-1)for(int i=0;i<n;i++)a[i]=a[i]/n;}void Pw(int *a,int *b,int *c){    for(int i=0;i<n;i++)ca[i]=cplxd(a[i],0);    for(int i=0;i<n;i++)cb[i]=cplxd(b[i],0);    FFT(ca,1),FFT(cb,1);    for(int i=0;i<n;i++)cc[i]=ca[i]*cb[i];    FFT(cc,-1);    for(int i=0;i<=M;i++)c[i]=((int)(0.3+cc[i].x))%P;}void pow(int k){    if(k==1)return;    pow(k>>1);    Pw(p,g,t),Pw(g,g,g);    for(int i=0;i<=M;i++)(p[i]+=t[i])%=P;    if(k&1){        Pw(g,F,g);        for(int i=0;i<=M;i++)(p[i]+=g[i])%=P;    }}signed main(){    scanf("%d%d%d%d%d%d",&M,&P,&A,&O,&S,&U);    for(n=1;n<=M*2;n<<=1)L++;    for(int i=1;i<=n;i++)R[i]=(R[i>>1]>>1)|((i&1)<<(L-1));    for(int i=1;i<=M;i++)p[i]=g[i]=F[i]=(i*i*O+S*i+U)%P;    pow(A),printf("%d\n",p[M]);}

 

BZOJ 4332 FFT+快速幂