省选专练之数学BZOJ3738 Ontak2013 Kapitał

和礼物差不多

但是题目中有一句话:

去除零

于是在阶乘的时候要搞事情

对逆元做阶乘

然后扩展Lucas水过

#include<bits/stdc++.h>
using namespace std;
typedef int INT;
#define int long long 
int n,m,k;
int mod;
int f[11]={};
int fac2[513]={};
int fac5[1953130]={};
//int Quick_Power(int x,int k,int p){
//	int ret=1;
//	while(k){
//		if(k%2==1)ret=ret*x%p;
//		k/=2;
//		x=x*x%p;
//	}
//	return ret;
//}
//void exgcd(int a,int b, int &x,int &y){
//	if(b==0)x=1,y=0;
//	else exgcd(b,a%b,y,x),y-=(a/b)*x;
//}
//int Inv(int a,int p){
//	int x,y;
//	exgcd(a,p,x,y);
//	return (x+p)%p;
//}
//int Fac(int n,int p,int np){
//	if(n==0)return 1;
//	int ret=1;
//	for(int i=2;i<=np;i++)if(i%p)ret=ret*i%np;
//	ret=Quick_Power(ret,n/np,np);
//	int r=n%np;
//	for(int i=2;i<=r;i++)if(i%p)ret=ret*i%np;
//	return ret*(Fac(n/p,p,np))%np;
//}
//int C(int n,int m,int p,int np){
	cout<<n<<" "<<m<<" "<<p<<" "<<np<<'\n';
//	int x=Fac(n,p,np);
//	int y=Fac(m,p,np);
//	int z=Fac(n-m,p,np);
//	int cnt=0;
//	for(int i=n;i;i/=p)cnt+=i/p;
//	for(int i=m;i;i/=p)cnt-=i/p;
//	for(int i=n-m;i;i/=p)cnt-=i/p;
//	int ret=x*Inv(y,np)%np*Inv(z,np)%np*Quick_Power(p,cnt,np)%np;
//	return ret*(mod/np)%mod*Inv(mod/np,np)%mod;
//}
//int Lucas(int n,int m){
//	int x=mod;
//	int ret=0;  
//	for(int i=2;i<=x;i++){
//		if(x%i==0){
//			int pr=1;
//			while(x%i==0)x/=i,pr*=i;
//			ret=(ret+C(n,m,i,pr)%mod)%mod;
//		}
//	}
//	return ret;
//}
int Pow(int a,int k,int p){
	int ret=1;
	while(k){
		if(k%2==1){
			ret=ret*a%p;
		}
		k/=2;
		a=a*a%p;
	}
	return ret;
}
//void exgcd(int a,int b,int &d,int &x,int &y){
//	if(b==0)d=a,x=1,y=0;
//	else exgcd(b,a%b,d,y,x),y-=(a/b)%x;
//}
void exgcd(int a,int b,int &d,int &x,int &y){
    if(b==0) d=a,x=1,y=0;
    else exgcd(b,a%b,d,y,x),y-=(a/b)*x;
}
int Inv(int a,int n){
	int d,x,y;
	exgcd(a,n,d,x,y);
//	return (x+n)%n;
	return (d==1)?(x+n)%n:-1;
}
int Fac(int n,int p,int pr){
	if(n==0)return 1;
	int ret=1;
	if(p==2)ret=fac2[pr];
	else ret=fac5[pr];
//	for(int i=2;i<=pr;i++)if(i%p)ret=ret*i%pr;
	ret=Pow(ret,n/pr,pr);
	int r=n%pr;
	if(p==2)ret=ret*fac2[r]%pr;
	else ret=ret*fac5[r]%pr;
//	for(int i=2;i<=r;i++)if(i%p)ret=ret*i%pr;
	return ret*Fac(n/p,p,pr)%pr;
}
int get(int n,int m,int p)
{
    int cnt=0;
    for(int i=n;i;i/=p)cnt+=(i/p);
    for(int i=m;i;i/=p)cnt-=(i/p);
    for(int i=n-m;i;i/=p)cnt-=(i/p);
    return cnt;
}
int C(int n,int m,int p,int pr){
	if(n<m)return 0;
	int a=Fac(n,p,pr);
	int b=Fac(m,p,pr);
	int c=Fac(n-m,p,pr);
	int c2,c5;
	c2=get(n,m,2);
	c5=get(n,m,5);
	int cnt;
	if(p==2)cnt=c2;
	else cnt=c5;
	int ret;
    if(p==2)
    {
        if(c5>=c2)ret=a*Inv(b,pr)%pr*Inv(c,pr)%pr*Pow(Inv(5,pr),c2,pr)%pr;
        else ret=a*Inv(b,pr)%pr*Inv(c,pr)%pr*Pow(Inv(5,pr),c5,pr)%pr*Pow(2,c2-c5,pr)%pr;
    }
    else
    {
        if(c2>=c5)ret=a*Inv(b,pr)%pr*Inv(c,pr)%pr*Pow(Inv(2,pr),c5,pr)%pr;
        else ret=a*Inv(b,pr)%pr*Inv(c,pr)%pr*Pow(Inv(2,pr),c2,pr)%pr*Pow(5,c5-c2,pr)%pr;
    }
	return ret*(mod/pr)%mod*Inv(mod/pr,pr)%mod; 
}
int Lucas(int n,int m){
	int x=mod;
	int ret=0;
	for(int i=2;i<=x;i++){
		if(x%i==0){
			int pr=1;
			while(x%i==0)pr*=i,x/=i;
			ret=(ret+C(n,m,i,pr)%mod)%mod;
		}
	}
	return ret;
}
INT main(){
	f[0]=1;
	for(int i=1;i<=10;i++)f[i]=f[i-1]*10;
    fac2[0]=1;
	for(int i=1;i<=512;i++)
   		fac2[i]=fac2[i-1]*((i&1)?i:1)%512;
    fac5[0]=1;
    for(int i=1;i<=1953125;i++)
    	fac5[i]=fac5[i-1]*((i%5)?i:1)%1953125;
	cin>>n>>m>>k;
	mod=f[k];
	int ans=Lucas(n+m,n);
	int w=0;
	int tmp=ans;
	while(tmp)tmp/=10,w++;
    for(int i=1;i<=k-w;i++)printf("0");
    printf("%lld",ans);
}

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值