HDU 3944 DP?------Lucas 大组合数取余

// Lucas定理求解C(n,k)%p, p是素数
// C(n,k)%p = C(n%p,k%p)*C(n/p,k/p)%p
// C(n,k)%p = n!/(k!*(n-k)!)%p = n!*((k!*(n-k)!)^(p-2))%p
// pow_mod(a,b,p)求a^b%p的值, 快速幂取模算法
// com_mod(n,k,p)求C(n,k)%p的值
// lucas(n,k,p)使用lucas定理求C(n,k)%p的值


#include<stdio.h>
#include<string.h>
#include<algorithm>
using namespace std;
const int maxn=10005;
int tot;
int prime[1250];
bool flag[maxn];
int fac[maxn][1250];
int Mul(int a,int b,int p)
{
	int ret=1;
	while(b)
	{
		if(b&1)ret=ret*a%p;
		a=a*a%p;
		b>>=1;
	}
	return ret;
}
void Init()
{
	int i,j;
	memset(flag,0,sizeof(flag));
	tot=0;
	for(i=2;i<maxn;i++)
		if(!flag[i])
		{
			prime[tot++]=i;
			for(j=i*i;j<maxn;j+=i)
				flag[j]=1;
		}
}
void Gao()
{
	Init();
	int i,j;
	for(i=0;i<1250;i++)
		fac[0][i]=1;	
	for(i=0;i<tot;i++)
	{
		for(j=1;j<maxn;j++)
		{
			fac[j][i]=fac[j-1][i]*j%prime[i];
		}
	}
}
int Lucas(int a,int b,int k,int p)
{
	int ret=1;
	while(a && b)
	{
		int ta=a%p;
		int tb=b%p;
		if(tb>ta) return 0;
		ret=ret*fac[ta][k]%p;
		ret=ret*Mul(fac[tb][k]*fac[ta-tb][k]%p,p-2,p)%p;
		a/=p;
		b/=p;
	}
	return ret;
}

int Bin(int x)
{
	int l=0,r=tot-1;
	while(l<=r)
	{
		int mid=(l+r)>>1;
		if(prime[mid]<x) l=mid+1;
		else if(prime[mid]>x) r=mid-1;
		else return mid;
	}
}

int main()
{
	int cas=1;
	int i,j,k;
	int a,b,p;
	Gao();
	while(scanf("%d%d%d",&a,&b,&p)!=EOF)
	{
		if(b>a-b)b=a-b;
		k=Bin(p);
		//k=(int)(lower_bound(prime,prime+tot,p)-prime);
		printf("Case #%d: %d\n",cas++,(Lucas(a+1,b,k,p)+a-b)%p);
	}
	return 0;
}