codeforces 100633J. Ceizenpok’s formula （扩展Lucas定理+中国剩余定理）

最新推荐文章于 2019-09-04 23:36:18 发布

clover_hxy

最新推荐文章于 2019-09-04 23:36:18 发布

阅读量589

点赞数

分类专栏：数论组合数取模

本文链接：https://blog.youkuaiyun.com/clover_hxy/article/details/55106859

版权

数论同时被 2 个专栏收录

158 篇文章

订阅专栏

组合数取模

17 篇文章

订阅专栏

Dr. Ceizenp'ok from planet i1c5l became famous across the whole Universe thanks to his recent discovery — the Ceizenpok’s formula. This formula has only three arguments: n, k and m, and its value is a number of k-combinations of a set of n modulo m.

While the whole Universe is trying to guess what the formula is useful for, we need to automate its calculation.

Input

Single line contains three integers n, k, m, separated with spaces (1 ≤ n ≤ 1018, 0 ≤ k ≤ n, 2 ≤ m ≤ 1 000 000).

Output

Write the formula value for given arguments n, k, m.

        Examples 
      
          input 
        
2 1 3

          output 
        
2

          input 
        
4 2 5

          output 
        
1

题解：扩展Lucas定理+中国剩余定理

快速乘要少用，死慢啊！！

各种地方炸掉，我也是醉了。。。

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
#define N 20
#define LL long long
using namespace std;
LL n,m,M,a[N];
LL p[1000003],mod[1000003],cnt;
void get_prime(LL k)
{
	LL x=k;
	for (LL i=2;i*i<=x;i++) 
	 if (x%i==0) {
	    ++cnt; p[cnt]=i; mod[cnt]=1;
	    while (x%i==0) {
	    	x/=i;
	    	mod[cnt]*=i;
		}
	 }
	if (x>1) {
		p[++cnt]=x;
		mod[cnt]=x;
	}
}
LL mul(LL a,LL b,LL mod)
{
	LL res=0;
	while (b) {
		if (b&1) res=(res+a)%mod;
		a=(a+a)%mod;
		b>>=1;
	}
	return res;
}
LL quickpow(LL num,LL x,LL p)
{
	LL base=num%p; LL ans=1;
	while (x) {
		if (x&1) ans=(ans*base)%p;
		x>>=1;
		base=(base*base)%p;
	}
	return ans;
}
void exgcd(LL a,LL b,LL &x,LL &y)
{
	if (b==0) {
		x=1; y=0; return;
	}
	exgcd(b,a%b,x,y);
	LL t=y; 
	y=x-(a/b)*y;
	x=t;
}
LL inverse(LL a,LL b)
{
	LL x,y;
	exgcd(a,b,x,y);
	x=(x%b+b)%b;
	if (!x) x+=b; 
	return x;
}
LL solve(int k1,LL n)
{
	if (n==0) return 1;
	LL ans=1;
	for (LL i=2;i<mod[k1];i++) 
	 if (i%p[k1]) ans=ans*i%mod[k1];
	ans=quickpow(ans,n/mod[k1],mod[k1]);
	for (LL i=2;i<=n%mod[k1];i++)
	 if (i%p[k1]) ans=ans*i%mod[k1];
	return ans*solve(k1,n/p[k1])%mod[k1];
}
LL calc(int k,LL n,LL m)
{
	if (m>n) return 0;
    LL a=solve(k,n),b=solve(k,m),c=solve(k,n-m);
    LL k1=0,ans;
    for (LL i=n;i;i/=p[k]) k1+=i/p[k];
    for (LL i=m;i;i/=p[k]) k1-=i/p[k];
    for (LL i=n-m;i;i/=p[k]) k1-=i/p[k];
    ans=a*inverse(b,mod[k])%mod[k]*inverse(c,mod[k])%mod[k]*quickpow(p[k],k1,mod[k])%mod[k];
    return ans;
}
LL china()
{
	LL x,y,ans=0;
    for (int i=1;i<=cnt;i++){
    	LL r=M/mod[i];
    	ans=(ans+r*inverse(r,mod[i])%M*a[i]%M)%M;
	}
	return ans%M;	
}
int main()
{
//	freopen("a.in","r",stdin);
//	freopen("my.out","w",stdout);
	cin>>n>>m>>M;
	get_prime(M);
	//for (int i=1;i<=cnt;i++) cout<<p[i]<<" "<<mod[i]<<endl;
	for (int i=1;i<=cnt;i++)
	 a[i]=calc(i,n,m);
	printf("%I64d\n",china());
}