BZOJ2629 : binomial

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