[BZOJ]3738: [Ontak2013]Kapitał 扩展lucas

Description
给出三个数字N,M,K。求C(N+M,N)去掉所有末尾的0后对10^K取模的结果.

1<=N,M<=10^15,1<=k<=9

题解：

扩展lucas直接上，去掉末尾的0就消去同样个数的2和5就行了，然后阶乘要预处理。

代码：

#include<bits/stdc++.h>
using namespace std;
#define LL long long
#define pa pair<int,int>
const LL M=1000000000LL;
LL fac2[520],fac5[1953130];
LL read()
{
    LL x=0,f=1;char ch=getchar();
    while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}
    while(ch>='0'&&ch<='9')x=(x<<3)+(x<<1)+(ch^48),ch=getchar();
    return x*f;
}
LL Pow(LL x,LL y,LL mod)
{
    if(y==0)return 1;
    if(y==1)return x;
    LL t=Pow(x,y>>1,mod),ans=t*t%mod;
    if(y&1)ans=ans*x%mod;
    return ans;
}
LL exgcd(LL a,LL b,LL &x,LL &y)
{
    if(!a)
    {
        x=0,y=1;
        return b;
    }
    LL tx,ty,d=exgcd(b%a,a,tx,ty);
    y=tx;
    x=ty-(b/a)*tx;
    return d;
}
LL inv(LL a,LL mod)
{
    if(!a)return 1;
    LL x,y,d=exgcd(a,mod,x,y);
    return (x%mod+mod)%mod;
}
LL fac(LL n,LL p,LL pk)
{
    if(!n)return 1;
    LL re;
    if(p==2)re=fac2[pk];
    else re=fac5[pk];
    re=Pow(re,n/pk,pk);
    if(p==2)re=re*fac2[n%pk]%pk;
    else re=re*fac5[n%pk]%pk;
    return re*fac(n/p,p,pk)%pk;
}
LL get(LL n,LL m,LL p)
{
    LL cnt=0;
    for(LL i=n;i;i/=p)cnt+=(i/p);
    for(LL i=m;i;i/=p)cnt-=(i/p);
    for(LL i=n-m;i;i/=p)cnt-=(i/p);
    return cnt;
}
LL C(LL n,LL m,LL p,LL pk)
{
    if(n<m)return 0;
    LL a=fac(n,p,pk),b=fac(m,p,pk),c=fac(n-m,p,pk),cnt;
    LL c2=get(n,m,2),c5=get(n,m,5),re=1;
    if(p==2)cnt=c2;else cnt=c5;
    if(p==2)
    {
        if(c5>=c2)return a*inv(b,pk)%pk*inv(c,pk)%pk*Pow(inv(5,pk),c2,pk)%pk;
        else return a*inv(b,pk)%pk*inv(c,pk)%pk*Pow(inv(5,pk),c5,pk)%pk*Pow(2,c2-c5,pk)%pk;
    }
    else
    {
        if(c2>=c5)return a*inv(b,pk)%pk*inv(c,pk)%pk*Pow(inv(2,pk),c5,pk)%pk;
        else return a*inv(b,pk)%pk*inv(c,pk)%pk*Pow(inv(2,pk),c2,pk)%pk*Pow(5,c5-c2,pk)%pk;
    }
}
LL lucas(LL n,LL m,LL mod)
{
    LL M=mod,re=0;
    for(LL i=2;i<=5;i++)
    if(mod%i==0)
    {
        LL t=1;
        while(mod%i==0)mod/=i,t*=i;
        re=(re+C(n,m,i,t)*inv(M/t,t)%M*(M/t)%M)%M;
    }
    if(mod>1)re=(re+C(n,m,mod,mod)*inv(M/mod,mod)%M*(M/mod)%M)%M;
    return re;
}
LL n,m,k;
int main()
{
    fac2[0]=fac5[0]=1;
    for(LL i=1;i<=512;i++)
    if(i%2)fac2[i]=fac2[i-1]*i%512LL;
    else fac2[i]=fac2[i-1];
    for(LL i=1;i<=1953125;i++)
    if(i%5)fac5[i]=fac5[i-1]*i%1953125LL;
    else fac5[i]=fac5[i-1];
    m=read(),n=read(),k=read();
    n+=m;
    LL f=1;for(int i=1;i<=k;i++)f*=10;
    LL t=lucas(n,m,M)%f;
    int w=0;LL tmp=t;while(tmp)tmp/=10,w++;
    for(int i=1;i<=k-w;i++)printf("0");
    printf("%lld",t);
}