bzoj 2629: binomial （FFT+DP+Lucas定理+短除法）

题目描述

传送门

题目描述：对于给定的n和p，求对于所有的$0<=i<p$,满足$C(n,k)\% \space p =i$的k的个数
　　注：$C(n,k)={n!\over{k!*(n-k)!}}$ $n<=p^{10},p=51061$

题解

这道题因为n很大所以不能直接计算组合数。
但是我们可以利用Lucas定理进行变形，$C(n,k)=C(n\%p,k\%p)*C(n/p,k/p)$ 不断的除p，就可以得到一串C相乘的形式。分解成一串$n\%p$其实就相当于是n的p进制数。因为$C(n,k)$在k>n的时候是等于0的，所以我们对于每一位来说只需要计算到p进制下这一位的值即可。这样组合出来的p进制数k一定是<=n的，对于剩余的情况都是余数为0的情况，最后用总数-$\sum_{i=1}^{p-1} ans[i]$即可得到$ans[0]$
另$f[i][j]$表示计算到p进制下的第i位，余数为j的方案数。
对于每个i，预处理$b[j]$ 表示$C(ni,mi)$ 中余数为j的个数
$f[i+1][j*k\%p]+=f[i][j]*b[k]$
如果我们求出p的原根，那么就可以把j*k变成num[j]+num[k]的形式，那么上面的式子就是一个卷积的形式，就可以用FFT来优化

代码

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
#define N 150000
#define LL long long
#define pi acos(-1)
using namespace std;
char ch[40]={'0','1','2','3','4','5','6','7','8','9','A','B','C','D','E','F','G','H','I','J','K','L','M','N','O','P','Q','R','S'};
struct data{
    double x,y;
    data(double X=0,double Y=0) {
        x=X,y=Y;
    }
}b[N],c[N];
data operator -(data a,data b) {
    return data(a.x-b.x,a.y-b.y);
}
data operator +(data a,data b){
    return data(a.x+b.x,a.y+b.y);
}
data operator *(data a,data b){
    return data(a.x*b.x-a.y*b.y,a.y*b.x+a.x*b.y);
}
int n,m,up[200],p,t[2000],k,l,a[2000],md[N],ans[N],f[20][N],R[N],g,num[N],L,mp[N];
LL jc[N],inv[N];
char s[2000];
void calc()
{
    for (l=0;k;) {
        for (int i=k-1;i>0;--i) {
            t[i-1]+=t[i]%p*10;
            t[i]/=p;
        }
        a[++l]=t[0]%p;
        t[0]/=p;
        for (;k>0&&!t[k-1];--k);
    }
} 
int quickpow(LL num,int x,LL p)
{
    LL ans=1; LL base=num%p;
    while (x) {
        if (x&1) ans=ans*base%p;
        x>>=1;
        base=base*base%p;
    }
    return ans;
}
LL C(int n,int m)
{
    return (LL)jc[n]*inv[m]%p*inv[n-m]%p;
}
int get_g(int x)
{
    if (x==2) return 1;
    for (int i=2;i;i++) {
        bool pd=1;
        for (int j=2;j*j<x;j++) 
         if (quickpow(i,(x-1)/j,x)==1) {
            pd=false;
            break;
         } 
        if (pd) return i;
    }
}
void FFT(data a[N],int opt)
{
    for (int i=0;i<n;i++)   
    if (i<R[i]) swap(a[R[i]],a[i]);  
    for (int i=1;i<n;i<<=1) {  
        double s=(double)pi/i;
        data wn=data(cos(s),opt*sin(s));  
        for (int p=i<<1,j=0;j<n;j+=p) {  
            data w=data(1,0);  
            for (int k=0;k<i;k++,w=w*wn) {  
                data x=a[j+k],y=w*a[j+k+i];  
                a[j+k]=x+y; a[j+k+i]=x-y;  
            }  
        }  
    }  
}
int main()
{
    freopen("tyrrell.in","r",stdin);
    freopen("tyrrell.out","w",stdout);
    scanf("%s",s);
    k=strlen(s);
    for (int i=k-1;i>=0;i--) t[k-i-1]=s[i]-'0';
    scanf("%d",&p);
    int cnt=0; 
    for (int i=k-1;i>=0;i--) {
        cnt=cnt*10+t[i];
        cnt%=29;
    }
    cnt=(cnt+1)%29;
    calc();
    g=get_g(p); 
    for (int i=0,j=1;i<p-1;i++,j=(j*g)%p) 
     num[j]=i,mp[i]=j;
    m=2*p;
    for (n=1;n<=m;n<<=1) L++;
    for (int i=0;i<n;i++) R[i]=(R[i>>1]>>1)|((i&1)<<(L-1));
    jc[0]=1;
    for (int i=1;i<=p;i++) jc[i]=(jc[i-1]*i)%p;
    for (int i=0;i<=p;i++) inv[i]=quickpow(jc[i],p-2,p)%p;
    f[0][0]=1;
    for (int i=1;i<=l;i++) {
        memset(md,0,sizeof(md));
        for (int j=0;j<=a[i];j++)
         md[num[C(a[i],j)%p]]++;
        for (int j=0;j<p;j++) md[j]%=29;
        for (int j=0;j<p;j++) b[j].x=md[j],b[j].y=0;
        for (int j=p;j<=n;j++) b[j].x=0,b[j].y=0;
        FFT(b,1);
        for (int j=0;j<p;j++) c[j].x=f[i-1][j],c[j].y=0;
        for (int j=p;j<=n;j++) c[j].x=0,c[j].y=0;
        FFT(c,1);
        for (int j=0;j<=n;j++) b[j]=b[j]*c[j];
        FFT(b,-1);
        for (int j=0;j<n;j++) f[i][j%(p-1)]+=(int)(b[j].x/n+0.5); 
        for (int j=0;j<p;j++) f[i][j]%=29;
    }
    int sum=0;
    for (int i=0;i<p;i++) ans[mp[i]]=f[l][i]%29,sum+=ans[mp[i]];
    sum%=29;
    ans[0]=(cnt-sum+29)%29;
    for (int i=0;i<p;i++) putchar(ch[ans[i]%29]);
    printf("\n");
}