扩展卢卡斯定理

引入

  有卢卡斯定理，$\left(\begin{array}{c}n\\ m\end{array}\right)\equiv \left(\begin{array}{c}n\mathrm{m}\mathrm{o}\mathrm{d}P\\ m\mathrm{m}\mathrm{o}\mathrm{d}P\end{array}\right)\times \left(\begin{array}{c}\lfloor \frac{n}{P}\rfloor \\ \lfloor \frac{m}{P}\rfloor \end{array}\right)(\mathrm{m}\mathrm{o}\mathrm{d}P)$。卢卡斯定理的适用条件是：$P$ 为质数。对于 $P$ 是一般正整数的情况，可以使用扩展卢卡斯定理计算 $C_n^m(\mathrm{m}\mathrm{o}\mathrm{d}P)$。

算法分析

原问题的转化

  要求 $C_n^m\mathrm{m}\mathrm{o}\mathrm{d}P$，首先对 $P$ 进行质因数分解，即 $P=p_1^{k_1}{p}_2^{k_2}\cdots p_{cnt}^{k_{cnt}}$，记 $a_i=C_n^m\mathrm{m}\mathrm{o}\mathrm{d}{p}_i^{k_i}$，可以获得一个一元线性同余方程组。
$$(S):\{\begin{array}{c}{C}_n^m\equiv a_1(\mathrm{m}\mathrm{o}\mathrm{d}{p}_1^{k_1})\\ C_n^m\equiv a_2(\mathrm{m}\mathrm{o}\mathrm{d}{p}_2^{k_2})\\ ⋮\\ C_n^m\equiv a_{cnt}(\mathrm{m}\mathrm{o}\mathrm{d}{p}_{cnt}^{k_{cnt}})\end{array}$$
  由于 $p_1,p_2,\cdots ,p_{cnt}$ 是互不相同的质数，因此 $p_1^{k_1},p_2^{k_2},\cdots ,p_{cnt}^{k_{cnt}}$ 两两互质，由中国剩余定理，其解集为 $A=\left\{x\in \mathbb{Z}\mid kP+\sum _{i=1}^{cnt}{a}_i{t}_i{P}_i,k\in \mathbb{Z}\right\}$；其中，$P_i=\frac{P}{p_i^{k_i}}$，$P_i{t}_i\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}P)$。尽管求出的 $C_n^m$ 并不是准确值，仅仅确定了 $C_n^m\in A$，但是可以注意到，此集合中的元素在模 $P$ 意义下的值都为 $\left(\sum _{i=1}^{cnt}{a}_i{t}_i{P}_i\right)\mathrm{m}\mathrm{o}\mathrm{d}P$；根据逻辑三段论，$C_n^m\equiv \sum _{i=1}^{cnt}{a}_i{t}_i{P}_i(\mathrm{m}\mathrm{o}\mathrm{d}P)$。所以，一元线性同余方程组 $\{\begin{array}{c}x\equiv a_1(\mathrm{m}\mathrm{o}\mathrm{d}{p}_1^{k_1})\\ x\equiv a_2(\mathrm{m}\mathrm{o}\mathrm{d}{p}_2^{k_2})\\ ⋮\\ x\equiv a_{cnt}(\mathrm{m}\mathrm{o}\mathrm{d}{p}_{cnt}^{k_{cnt}})\end{array}$的最小正整数解即为 $C_n^m\mathrm{m}\mathrm{o}\mathrm{d}P$。
  问题在于，要如何先求出 $a_i=C_n^m\mathrm{m}\mathrm{o}\mathrm{d}{p}_i^{k_i}$。

计算 $C_n^m\mathrm{m}\mathrm{o}\mathrm{d}{p}^k$

  下面讨论如何计算 $C_n^m\mathrm{m}\mathrm{o}\mathrm{d}{p}^k$，其中 $p$ 为质数。
  首先由组合数的公式 $C_n^m=\frac{n!}{m!(n-m)!}$，分别计算出 $n!,m!,(n-m)!$ 在模 $p^k$ 意义下的值，再求 $m!,(n-m)!$ 模 $p^k$ 意义下的逆元，就可以得到答案。但是 $m!$ 和 $(n-m)!$ 可能包含质因子 $p$，根据裴蜀定理，此时它们在模 $p^k$ 意义下的逆元不存在。不妨将$n!,m!,(n-m)!$ 中的质因子 $p$ 除尽，对 $C_n^m$ 做如下变换。
$${\text{C}}_n^m=\frac{\frac{n!}{p^x}}{\frac{m!}{p^y}\times \frac{(n-m)!}{p^z}}\times p^{x-y-z}$$
  上式中，$x$ 为 $n!$ 中质因子 $p$ 的个数，$y$ 为 $m!$ 中质因子 $p$ 的个数，$z$ 为 $(n-m)!$ 中质因子 $p$ 的个数。$\frac{n!}{p^x}$，$\frac{m!}{p^y}$，$\frac{(n-m)!}{p^z}$ 均与 $p^k$ 互质，都存在模 $p^k$ 意义下的逆元。
  以上，我们已经解决了阶乘求逆元的问题。只要求得 $\frac{n!}{p^x},\frac{m!}{p^y},\frac{(n-m)!}{p^z}$ 在模 $p^k$ 意义下的值，再求两次逆元，即可得到 $C_n^m\mathrm{m}\mathrm{o}\mathrm{d}{p}^k$。但是，问题在于：如何求形如 $\frac{n!}{p^x}\mathrm{m}\mathrm{o}\mathrm{d}{p}^k$ 的式子。

计算 $\frac{n!}{p^x}\mathrm{m}\mathrm{o}\mathrm{d}{p}^k$

  $n!=1\cdot 2\cdot 3\cdot \cdots \cdot n$=$(p\cdot 2p\cdot 3p\cdots )(1\cdot 2\cdot \cdots )$，左括号内为 $1\sim n$ 中 $p$ 的倍数的乘绩，右括号内反之。$1\sim n$ 中 $p$ 的倍数有 $\lfloor \frac{n}{p}\rfloor$ 个，故 $n!=p^{\lfloor \frac{n}{p}\rfloor }\lfloor \frac{n}{p}\rfloor !\prod _{i=1,p\nmid i}^{n}i$。
  计算 $\frac{n!}{p^x}$ 时，可以直接忽略 $p^{\lfloor \frac{n}{p}\rfloor }$ ，但是 $\lfloor \frac{n}{p}\rfloor !$ 还可能包含质因子 $p$，难以直接忽略其中的质因子；不妨定义函数 $fac(n,p,pk)=\frac{n!}{p^x}\mathrm{m}\mathrm{o}\mathrm{d}pk$，则 $fac(n,p,pk)=fac\left(\lfloor \frac{n}{p}\rfloor ,p,pk\right)\prod _{i=1,p\nmid i}^{n}i$，如此一来可以进行递归求解。
  $\left(\prod _{i=1,p\nmid i}^{n}i\right)\mathrm{m}\mathrm{o}\mathrm{d}{p}^k$ 是有循环节的，比如当 $p=3,t=2,n=19$ 时，有：$19!=3^6\times (1\times 2\times 3\times 4\times 5\times 6)\times (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)$，右括号中 $1\times 2\times 4\times 5\times 7\times 8\equiv 10\times 11\times 13\times 14\times 16\times 17(\mathrm{m}\mathrm{o}\mathrm{d}{3}^2)$，可以看到，$\left(\prod _{i=1,p\nmid i}^{n}i\right)\mathrm{m}\mathrm{o}\mathrm{d}{p}^k$ 以 $p^k$ 为一个周期。$\left(\prod _{i=1,p\nmid i}^{n}i\right)\mathrm{m}\mathrm{o}\mathrm{d}{p}^k$ 是以 $p^k$ 为周期的，每个循环节长度为 $p^k $，所以只需要计算最后不满足一个周期的数的乘积（例子中只需要计算 $19$）。显然，完整的周期数为 $\lfloor \frac{n}{p^k}\rfloor$， 每个周期模 $p^k$ 意义下的成绩为 $\prod _{i=1,p\nmid i}^{p^k}i\mathrm{m}\mathrm{o}\mathrm{d}{p}^k$ ；不满足一个周期的数的个数为 $n\mathrm{m}\mathrm{o}\mathrm{d}{p}^k$，模 $p^k$ 意义下的乘积为 $\prod _{i=1,p\nmid i}^{n\mathrm{m}\mathrm{o}\mathrm{d}{p}^k}i\mathrm{m}\mathrm{o}\mathrm{d}{p}^k$。

洛谷P4720 【模板】扩展卢卡斯

题目背景

  这是一道模板题。

题目描述

  求 $C_n^m\mathrm{m}\mathrm{o}\mathrm{d}P$，其中 $C$ 为组合数。

输入格式

  一行三个整数 $n,m,P$，含义由题所述。

输出格式

  一行一个整数，表示答案。

输入输出样例

输入 #1

5 3 3

输出 #1

1

输入 #2

666 233 123456

输出 #2

61728

说明/提示

  $1\le m\le n\le 10^{18}$，$2\le P\le 10^6$，不保证 $P$ 是质数。

代码

#include<iostream>
#define N 12
#define ll long long
using namespace std;
ll b[N],a[N];
int cnt;
//-------------------------------------------------------------------------------------------------------
//扩展欧几里得算法
//用于求逆元
void ex_gcd(ll a,ll b,ll &x,ll &y)
{
	if(!b) 
	{
		x=1;
		y=0;
		return;
	}
	ex_gcd(b,a%b,x,y);
	ll temp=x;
	x=y;
	y=temp-a/b*y;
}
//-------------------------------------------------------------------------------------------------------
//-------------------------------------------------------------------------------------------------------
//运用扩展欧几里得算法求逆元
//ax≡1(mod b)
ll inv(ll a,ll b)
{
	ll x,y;
	ex_gcd(a,b,x,y);
	//通解 x=x0+kb
	return (x%b+b)%b;
}
//-------------------------------------------------------------------------------------------------------
//-------------------------------------------------------------------------------------------------------
//快速幂并取模
ll qpow_mod(ll a,ll b,ll mod)
{
	ll res=1;
	while(b)
	{
		if(b&1) res=res*a%mod;
		a=a*a%mod;
		b>>=1;
	}
	return res%mod;
}
//-------------------------------------------------------------------------------------------------------
//-------------------------------------------------------------------------------------------------------
//阶乘除去质因子模质数幂
//pk表示p的k次幂
ll fac_mod(ll n,ll p,ll pk)
{
	if(!n) return 1;
	ll res=1;
	int i;
	//----------------------------------------
	//一个循环节内的乘积
	//周期为pk
	for(i=1;i<=pk;i++) if(i%p) res=res*i%pk;
	//n/pk个完整周期
	res=qpow_mod(res,n/pk,pk);
	//不满足一个周期的数的乘积
	for(i=1;i<=n%pk;i++) if(i%p) res=res*i%pk;
	//-----------------------------------------
	//-----------------------------------------
	//递归求 (n/p)!除去质因子p并对p的k次幂取模
	return res*fac_mod(n/p,p,pk)%pk;
	//-----------------------------------------
}
//-------------------------------------------------------------------------------------------------------
//-------------------------------------------------------------------------------------------------------
//组合数对p的k次幂取模
ll C_mod(ll n,ll m,ll p,ll pk)
{
	ll x,y,z;
	x=y=z=0;
	ll i;
	//------------------------------
	//阶乘除尽质因子
	//最后乘p的(x-y-z)次幂
	//n!中p质因子的数目
	for(i=p;i<=n;i*=p) x+=n/i;
	//m!中p质因子的数目
	for(i=p;i<=m;i*=p) y+=m/i;
	//(n-m)!中p质因子的数目
	for(i=p;i<=n-m;i*=p) z+=(n-m)/i;
	//-------------------------------
	//n!/m!/(n-m)!
	//求两次逆元
	//三次阶乘除以质因子模质数幂
	return fac_mod(n,p,pk)*inv(fac_mod(m,p,pk),pk)%pk*inv(fac_mod(n-m,p,pk),pk)%pk*qpow_mod(p,x-y-z,pk);
	//--------------------------------------------------------------------------------------------------
}
//-------------------------------------------------------------------------------------------------------
//中国剩余定理
ll CRT()
{
	int i;
	ll mul=1;
	ll res=0;
	for(i=1;i<=cnt;i++) mul*=a[i];
	for(i=1;i<=cnt;i++)
	{
		ll Mi=mul/a[i];
		ll x,y;
		res+=b[i]*inv(Mi,a[i])*Mi;
	}
	return res%mul;
}
//-------------------------------------------------------------------------------------------------------
//-------------------------------------------------------------------------------------------------------
//扩展卢卡斯定理
ll ex_Lucas(ll n,ll m,ll P)
{
	ll i;
	//-------------------------------------
	//分解质因数
	for(i=2;i*i<=P;i++)
	{
		if(P%i==0)
		{
			ll pk=1;
			while(P%i==0)
			{
				pk*=i;
				P/=i;
			}
			cnt++;//不同的质因子个数
			//----------------------------
			//质因子i
			//pk表示i的k次幂
			a[cnt]=pk;//模数
			b[cnt]=C_mod(n,m,i,pk);//余数
			//----------------------------
		}
	}
	if(P>1) 
	{
		cnt++;
		//质因子P
		//pk=P
		a[cnt]=P;
		b[cnt]=C_mod(n,m,P,P);
	}
	return CRT();
}
//-------------------------------------------------------------------------------------------------------
//-------------------------------------------------------------------------------------------------------
int main()
{
	ll n,m,P;
	cin>>n>>m>>P;
	cout<<ex_Lucas(n,m,P)<<endl;
	return 0;
}