(扩展)Lucas定理——求大组合数取模

Lucas定理

若$P$为质数
$$C_n^r\equiv C_{\lfloor n/P\rfloor }^{\lfloor r/P\rfloor }\times C_{nmodP}^{rmodP}(modP)$$
要用生成函数证明（不会，所以不证了）

实现直接递归：

void exgcd(long long a,long long b,long long &x,long long &y)
{
	if(b==0)
		x=1,y=0;
	else
	{
		long long x1,y1;
		exgcd(b,a%b,x1,y1);
		x=y1;
		y=x1-(a/b)*y1;
	}
}
long long inv(long long x,long long p)
{
	long long i,t;
	exgcd(x,p,i,t);
	i=(i%p+p)%p;
	return i;
}


long long C(long long n,long long r,long long MOD)
{
	if(n<r)
		return 0;
	if(r==0)
		return 1;
	long long ans=1LL;
	ans=C(n/MOD,r/MOD,MOD);
	n%=MOD,r%=MOD;
	long long a=1,b=1;
	for(int i=n;i>n-r;i--)
		a=(a*i)%MOD;
	for(int i=1;i<=r;i++)
		b=(b*i)%MOD;
	ans=(ans*a%MOD*inv(b,MOD))%MOD;
	return ans;
}

扩展Lucas定理

可以解决模数$P$不是质数的组合数取模

将$P$质因数分解：$P=p_1^{k_1}{p}_2^{k_2}...p_t^{k_t}$
如果我们求出一下方程中的$a_1,a_2...a_t$，就可以利用中国剩余定理求出$C_n^r$
$$\{\begin{array}{c}{C}_n^r\equiv a_1(mod{p}_1^{k_1})\\ C_n^r\equiv a_2(mod{p}_2^{k_2})\\ C_n^r\equiv a_3(mod{p}_3^{k_3})\\ ...\\ C_n^r\equiv a_t(mod{p}_t^{k_t})\end{array}$$

现在来求$C_n^r\equiv a(mod{p}^k)$
$$C_n^r\equiv \frac{n!}{(n-r)!r!}(mod{p}^k)$$
令$a\equiv n!(mod{p}^k)$，$b\equiv r!(mod{p}^k)$，$c\equiv (n-r)!(mod{p}^k)$

将组合数分成三个阶乘计算：$C_n^r\equiv a{b}^{-1}{c}^{-1}(mod{p}^k)$
（$b^{-1}$表示逆元）
但此时如果$b$与$p^k$不互质，就无法求出逆元，所以需要把$a$,$b$,$c$中的质因子$p$全部提出来，求了逆元之后再乘回去。

举个例子：计算 $19!(mod{3}^2)$
$$\begin{array}{rl}19!& =1\times 2\times 3\times 4\times ...\times 19\\ & =(1\times 2\times 4\times 5\times 7\times 8\times 10\times 11\times 13\times 14\times 16\times 17\times 19)\times 3^6\times (1\times 2\times 3\times 4\times 5\times 6)\end{array}$$
即把阶乘中的$p$倍数提出来，后面$6!$继续递归求解；
第一个括号里的式子有周期性的：
$$1\times 2\times 4\times 5\times 7\times 8\equiv 10\times 11\times 13\times 14\times 16\times 17(mod{3}^2)$$
所以暴力求出一个区间的值$t$，结果乘$t^{\lfloor \frac{19}{9}\rfloor }$

细节见代码：

long long pow_mod(long long a,long long b,long long p)
{
	long long res=1;
	while(b)
	{
		if(b&1LL)
			res=(res*a)%p;
		b>>=1LL;
		a=(a*a)%p;
	}
	return res;
}
void exgcd(long long a,long long b,long long &x,long long &y)
{
	if(b==0)
	{
		x=1,y=0;
		return;
	}
	exgcd(b,a%b,y,x);
	y-=(a/b)*x;
}
long long inv(long long a,long long p)//逆元必须用扩展欧几里得，因为p不为质数
{
	long long i,t;
	exgcd(a,p,i,t);
	i=(i%p+p)%p;
	return i;
}

long long Fac(long long N,long long p,long long pk)
{
	if(N==0)
		return 1;
	long long ans=1;
	if(N>=pk)//周期性
	{
		for(long long i=1;i<pk;i++)
			if(i%p)
				ans=(ans*i)%pk;
		ans=pow_mod(ans,N/pk,pk);
	}
	for(long long i=1;i<=N%pk;i++)//除了周期性的，剩下的几个数
		if(i%p)
			ans=(ans*i)%pk;
	ans=(ans*Fac(N/p,p,pk))%pk;//递归求解
	return ans;
}
long long C(long long N,long long R,long long p,long long pk)
{
	if(R==0)
		return 1;
	long long a=Fac(N,p,pk),b=Fac(N-R,p,pk),c=Fac(R,p,pk),k=0,ans;//求出几个阶乘，不包含质因子p
	//下面单独计算质因子p的个数
	for(long long i=N;i;i/=p)
		k+=i/p;
	for(long long i=N-R;i;i/=p)
		k-=i/p;
	for(long long i=R;i;i/=p)
		k-=i/p;
	ans=a*inv(b,pk)%pk*inv(c,pk)%pk*pow_mod(p,k,pk)%pk;
	return ans;
}
long long C(long long N,long long R,long long M)
{
	if(N<R)
		return 0;
	long long ans=0,m=M,temp;
	for(long long p=2,pk;p*p<=m;p++)
		if(m%p==0)//枚举模数M的因子p
		{
			pk=1;
			while(m%p==0)
				m/=p,pk*=p;
			temp=C(N,R,p,pk);//求C(N,R)mod p^k
			ans=(ans+temp*(M/pk)%M*inv(M/pk,pk)%M)%M;//中国剩余定理
		}
	if(m!=1)
	{
		temp=C(N,R,m,m);
		ans=(ans+temp*(M/m)%M*inv(M/m,m)%M)%M;
	}
	return ans;
}