bzoj3491: PA2007 Subsets

#include<cstdio>
#include<cmath>
#include<algorithm>
typedef unsigned long long u64;
typedef long double ld;
const ld pi=std::acos(-1);
struct cplx{
    ld a,b;
    cplx(ld x=0,ld y=0):a(x),b(y){}
    cplx operator+(cplx x){return cplx(a+x.a,b+x.b);}
    cplx operator-(cplx x){return cplx(a-x.a,b-x.b);}
    cplx operator*(cplx x){return cplx(a*x.a-b*x.b,a*x.b+b*x.a);}
    cplx conj(){return cplx(a,-b);}
}A[2111],B[2111],E[2][13][1111];
u64 n,pwx[77],ts[77];
int m,k,x,ppx=0,C[77][77],N,K;
int ans[1111],tmp[1111],tmp2[1111],rev[2111];
void dft(cplx*a,int t){
    for(int i=0;i<N;++i)if(i<rev[i])std::swap(a[i],a[rev[i]]);
    for(int i=1,z=0;i<N;i<<=1,++z){
        cplx*e=E[t][z];
        for(int j=0;j<N;j+=i<<1){
            cplx*b=a+j,*c=b+i;
            for(int k=0;k<i;++k){
                cplx x=b[k],y=c[k]*e[k];
                b[k]=x+y;
                c[k]=x-y;
            }
        }
    }
    if(t)for(int i=0;i<N;++i)a[i].a/=N;
}
void mul(int*a,int*b){
    for(int i=0;i<=k;++i)B[i]=cplx(a[i],b[i]);
    for(int i=k+1;i<N;++i)B[i]=0;
    dft(B,0);
    B[N]=B[0];
    for(int i=0;i<N;++i){
        cplx x=B[i],y=B[N-i].conj();
        A[i]=cplx(0,0.25)*(y+x)*(y-x);
    }
    dft(A,1);
    for(int i=0;i<=k;++i)a[i]=((long long)(A[i].a+0.49))%m;
}
int main(){
    scanf("%llu%d%d%d",&n,&m,&k,&x);
    C[0][0]=1;
    for(int i=0;i<70;++i){
        for(int j=0;j<=i;++j){
            (C[i+1][j]+=C[i][j])%=m;
            (C[i+1][j+1]+=C[i][j])%=m;
        }
    }
    for(pwx[ppx++]=1;(pwx[ppx]=pwx[ppx-1]*x)<=n;++ppx);
    for(int i=ppx-1;i>=0;--i){
        ts[i]=n/pwx[i];
        ts[i]-=ts[i]/x;
    }
    for(N=2,K=0;N<k*2+5;N<<=1,++K);
    for(int i=1;i<N;++i)rev[i]=rev[i>>1]>>1|(i&1)<<K;
    for(int i=1,z=0;i<N;i<<=1,++z){
        for(int j=0;j<i;++j){
            E[0][z][j]=cplx(cos(j*pi/i),sin(j*pi/i));
            E[1][z][j]=E[0][z][j].conj();
        }
    }
    ans[0]=1;
    for(int i=0;i<ppx;++i){
        ts[i]-=ts[i+1];
        for(int j=0;j<=k;++j)tmp[j]=(i+2-j>=0?C[i+2-j][j]:0),tmp2[j]=0;
        tmp2[0]=1;
        for(u64 t=ts[i];t;t>>=1,mul(tmp,tmp))if(t&1)mul(tmp2,tmp);
        mul(ans,tmp2);
    }
    printf("%d\n",ans[k]);
    return 0;
}