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

来源:互联网 发布:淘宝保证金计划打不开 编辑:程序博客网 时间:2024/06/05 06:58

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);}
原创粉丝点击