BZOJ2629 : binomial

时间：2024-01-07 14:50:14

根据Lucas定理，等价于在$P$进制下每一位分别求组合数最后乘积模$P$。

因为答案为$0$的并不好算，所以可以考虑用$n+1$减去其它所有的答案。

那么每一位的组合数都不能是$0$，那么这就保证了$k$的每一位都不大于$n$，所以无需考虑$k\leq n$这个限制。

求出模$P$下每个数的指标，考虑数位DP，设$f[i][j]$表示考虑了最后$i$位，组合数的指标之和模$\varphi(P)$为$j$的方案数。

那么转移就是卷积的形式，FFT加速即可。

时间复杂度$O(P\log P\log n)$。

#include<cstdio>

#include<cstring>

#include<algorithm>

#include<cmath>

using namespace std;

const int P=51061,F=P-1,O=29,N=131100;

const double pi=acos(-1.0);

char s[P];

int flag,n=131072,len,i,j,po[P],ind[P],fac[P],inv[P],a[P],ans[P],pos[N],f[N],g[N];

struct comp{

  double r,i;comp(double _r=0,double _i=0){r=_r,i=_i;}

  comp operator+(const comp&x){return comp(r+x.r,i+x.i);}

  comp operator-(const comp&x){return comp(r-x.r,i-x.i);}

  comp operator*(const comp&x){return comp(r*x.r-i*x.i,r*x.i+i*x.r);}

  comp conj(){return comp(r,-i);}

}A[N],B[N];

inline int C(int n,int m){return 1LL*fac[n]*inv[m]*inv[n-m]%P;}

void FFT(comp*a,int t){

  for(int i=1;i<n;i++)if(i<pos[i])swap(a[i],a[pos[i]]);

  for(int d=0;(1<<d)<n;d++){

    int m=1<<d,m2=m<<1;double o=pi*2/m2*t;comp _w(cos(o),sin(o));

    for(int i=0;i<n;i+=m2){

      comp w(1,0);

      for(int j=0;j<m;j++){

        comp&A=a[i+j+m],&B=a[i+j],t=w*A;

        A=B-t;B=B+t;w=w*_w;

      }

    }

  }

  if(t==-1)for(int i=0;i<n;i++)a[i].r/=n;

}

void mul(int k){

  if(!flag){

    flag=1;

    for(i=0;i<=k;i++)f[ind[C(k,i)]]++;

    for(i=0;i<F;i++)f[i]%=O;

    return;

  }

  for(i=0;i<F;i++)g[i]=0;

  for(i=0;i<=k;i++)g[ind[C(k,i)]]++;

  for(i=0;i<F;i++)g[i]%=O;

  for(i=0;i<F;i++)A[i]=comp(f[i],g[i]);

  for(i=F;i<n;i++)A[i]=comp(0,0);

  FFT(A,1);

  for(i=0;i<n;i++){

    j=(n-i)&(n-1);

    B[i]=(A[i]*A[i]-(A[j]*A[j]).conj())*comp(0,-0.25);

  }

  FFT(B,-1);

  for(i=0;i<F;i++)f[i]=0;

  for(i=0;i<n;i++)f[i%F]=(f[i%F]+int(B[i].r+0.5))%O;

}

int main(){

  for(po[0]=i=1;i<P-1;i++)po[i]=po[i-1]*2%P;

  for(i=0;i<P-1;i++)ind[po[i]]=i;

  for(fac[0]=fac[1]=inv[0]=inv[1]=1,i=2;i<P;i++)inv[i]=1LL*(P-inv[P%i])*(P/i)%P;

  for(i=2;i<P;i++)fac[i]=1LL*fac[i-1]*i%P,inv[i]=1LL*inv[i-1]*inv[i]%P;

  scanf("%s",s+1);

  len=strlen(s+1);

  for(i=1;i<=len;i++)a[i]=s[len-i+1]-'0',ans[0]=(ans[0]*10+s[i]-'0')%O;

  ans[0]++;

  ans[0]%=O;

  j=__builtin_ctz(n)-1;

  for(i=0;i<n;i++)pos[i]=pos[i>>1]>>1|((i&1)<<j);

  while(len){

    for(a[0]=0,i=len;i;i--)a[i-1]+=a[i]%P*10,a[i]/=P;

    while(len&&!a[len])len--;

    mul(a[0]/10);

  }

  for(i=0;i<F;i++)ans[po[i]]=f[i];

  for(i=1;i<P;i++)ans[0]=(ans[0]+O-ans[i])%O;

  for(i=0;i<P;i++)if(ans[i]<=9)putchar(ans[i]+'0');else putchar(ans[i]-10+'A');

  return 0;

}

　　

相关文章