Bzoj2629 binomial

  1 #include<iostream>
  2 #include<algorithm>
  3 #include<cstring>
  4 #include<cstdio>
  5 #include<cmath>
  6 #define LL long long
  7 using namespace std;
  8 const double pi=acos(-1.0);
  9 const int mod=29;
 10 const int P=51061;
 11 const int ED=P-1;
 12 const int mxn=150005;
 13 int read(){
 14     int x=0,f=1;char ch=getchar();
 15     while(ch<‘0‘ || ch>‘9‘){if(ch==‘-‘)f=-1;ch=getchar();}
 16     while(ch>=‘0‘ && ch<=‘9‘){x=x*10+ch-‘0‘;ch=getchar();}
 17     return x*f;
 18 }
 19 int fac[mxn],inv[mxn];
 20 void init(){
 21     fac[0]=fac[1]=1;inv[0]=inv[1]=1;
 22     for(int i=2;i<mxn;i++){
 23         fac[i]=(LL)fac[i-1]*i%P;
 24         inv[i]=((-P/i*(LL)inv[P%i]%P)+P)%P;
 25     }
 26     for(int i=2;i<mxn;i++)
 27         inv[i]=(LL)inv[i-1]*inv[i]%P;
 28     return;
 29 }
 30 //
 31 struct com{
 32     double x,y;
 33     com(){}
 34     com(double _x,double _y):x(_x),y(_y){}
 35     com operator + (const com &b){return com(x+b.x,y+b.y);}
 36     com operator - (const com &b){return com(x-b.x,y-b.y);}
 37     com operator * (const com &b){return com(x*b.x-y*b.y,x*b.y+y*b.x);}
 38     com operator / (const double v){return com(x/v,y/v);}
 39 }A[mxn],B[mxn];
 40 int N,len,rev[mxn];
 41 void FFT(com *a,int flag){
 42     for(int i=0;i<N;i++)if(i<rev[i])swap(a[i],a[rev[i]]);
 43     for(int i=1;i<N;i<<=1){
 44         int p=i<<1;
 45         com wn=com(cos(pi/i),flag*sin(pi/i));
 46         for(int j=0;j<N;j+=p){
 47             com w=com(1,0);
 48             for(int k=0;k<i;k++,w=w*wn){
 49                 com x=a[j+k],y=w*a[j+k+i];
 50                 a[j+k]=x+y;
 51                 a[j+k+i]=x-y;
 52             }
 53         }
 54     }
 55     if(flag==-1)
 56         for(int i=0;i<N;i++) a[i].x/=N;
 57     return;
 58 }
 59 int C(int n,int m){
 60     if(n<m)return 0;
 61     return (LL)fac[n]*inv[m]%P*inv[n-m]%P;
 62 }
 63 int pw[mxn],ind[mxn];
 64 int f[mxn];
 65 bool flag=0;
 66 void solve(int x){
 67 //    printf("x:%d\n",x);
 68     memset(B,0,sizeof B);
 69     for(int i=0;i<=x;i++)
 70         B[ind[C(x,i)]].x++;
 71     for(int i=0;i<ED;i++){
 72         A[i].x=f[i]%mod;A[i].y=0;
 73         B[i].x=((int)B[i].x)%mod;
 74     }
 75     for(int i=ED;i<N;i++)A[i].x=A[i].y=0;
 76     FFT(A,1);FFT(B,1);
 77     for(int i=0;i<N;i++)A[i]=A[i]*B[i];
 78     FFT(A,-1);
 79     memset(f,0,sizeof f);
 80     for(int i=0;i<N;i++){
 81         (f[i%ED]+=(int)(A[i].x+0.5)%mod)%=mod;
 82     }
 83     return;
 84 }
 85 char s[mxn];
 86 int num[mxn];
 87 int ans[mxn];
 88 int main(){
 89     int i,j;
 90     init();
 91     pw[0]=1;
 92     for(i=1;i<ED;i++)pw[i]=pw[i-1]*2%P;
 93     for(i=0;i<ED;i++)ind[pw[i]]=i;
 94     //
 95     int m=P<<1;
 96     for(N=1,len=0;N<=m;N<<=1)len++;
 97     for(i=0;i<N;i++)rev[i]=(rev[i>>1]>>1)|((i&1)<<(len-1));
 98     scanf("%s",s+1);
 99     int len=strlen(s+1);
100     for(i=1;i<=len;i++){
101         num[i]=s[len-i+1]-‘0‘;
102         ans[0]=(ans[0]*10+(s[i]-‘0‘))%mod;
103     }
104     ans[0]=(ans[0]+1)%mod;
105     f[0]=1;
106     while(len){
107         num[0]=0;
108         for(i=len;i;i--)num[i-1]+=num[i]%P*10,num[i]/=P;
109         while(len && !num[len])len--;
110         solve(num[0]/10);
111     }
112     for(i=0;i<ED;i++)ans[pw[i]]=f[i]%mod;
113     for(i=1;i<P;i++){
114         ans[0]=(ans[0]-ans[i])%mod;
115     }
116     ans[0]=(ans[0]+mod)%mod;
117     for(i=0;i<P;i++){
118         if(ans[i]>9)putchar(‘A‘+ans[i]-10);
119             else putchar(‘0‘+ans[i]);
120     }
121     return 0;
122 }