大组合数问题

组合数的公式：$$C_n^r=\frac{n\ast (n-1)\ast ...\ast (n-r+1)}{r!}=\frac{n!}{r!(n-r)!}$$观察上述式子，读者应该发现了，我们主要还是采用第一个算式，因为乘法操作相对较少，速度比较快，也相对来说不容易溢出。

思路一

暴力做乘法，并且使用组合数本身的性质进行优化，使用long long，保证在乘法过程中可以连乘更多的数。另外我们也可以优先做除法，再做乘法。当然这种优化方式能带来的支持运算的范围增大是相当小的。
A. $C_n^0=C_n^n=1$
B. $C_n^r=C_n^{n-r}$

#include<bits/stdc++.h>
using namespace std;
int main()
{
	int n,r;
	long long mon=1,son=1; 
	cin>>n>>r;
	if(r==0||r==n) cout<<1;
	if(r>n-r) r=n-r;
	for(int i=1;i<=r;i++,--n)
	{
		son*=i;
		mon*=n;
	}
	cout<<mon/son;
	return 0;
} 

思路二

其实我们就是要想办法让数在连乘过程中不要因为达到上界发生溢出而导致错误，所以我们可以考虑在做乘法的过程中同时除。这样会导致当前没法整除，而在计算机中，整数中的整除是下取整的，所以用整数直接做除法是不正确的。解决办法有三个：

1. 我们将除数做缓存，每次可以整除的时候，拿出来除，缩小被除数，并从缓存表中删去这个除数
2. 直接用浮点数存，虽然浮点数会有精度损失，但是损失的量是十分小的，而组合数的结果一定是整数，因此我们将浮点结果上取整就好啦~~ ( 但是double类型的存储范围要比long long广，不过需要同时考虑到浮点数本身带来的精度损失，因此我们只当作一种思维拓展，而不把他们视为正规解法 )。
3. .我们如果对上面的公式两边取对数，就可以得到$$\ln (C_n^r)=\ln (n!)-\ln (r!)-\ln ((n-r)!)=\sum _{i=1}^n\ln (i)-\sum _{i=1}^r\ln (i)-\sum _{i=1}^{n-r}\ln (i)$$使用对数的性质可以帮助我们压缩存储范围，同时改乘除法为加减法，但问题是计算机本身的存储浮点数带来的精度损失问题。

思路三

套大整数的模板。但是大整数采用数组实现，所以速度是比较慢的，另外，有时会碰到取模的要求，在大整数类型中做除法和取模操作是十分不方便的。

思路四

采用杨辉三角递推式。对于组合数，有这么一个性质：$$C_n^r=C_{n-1}^r+C_{n-1}^{r-1}$$我们以取物问题为例，理解上述公式：$C_n^r$ 表示从n个物体中拿出r个，那么如果我们针对第n个物体考虑，可以将这个原问题划分为两个小的子问题：

1. 我们将第n个物体当作我们拿出的第r个物体，此时我们需要在前n-1个物体中拿出r-1个，即$C_{n-1}^r$；
2. 我们保证不拿第n个物体，此时我们需要在前n-1个物体中拿出r个，即$C_{n-1}^{r-1}$

#include<bits/stdc++.h>
using namespace std;
int C[1010][1010];
int main()
{
	int n,r,i,j;
	cin>>n>>r;
	C[1][0]=C[1][1]=1;
	for(i=2;i<=n;++i)
	{
		C[i][0]=C[i][i]=1;
		for(j=1;j<n;j++) C[i][j]=C[i-1][j]+C[i-1][j-1];
	}
	cout<<C[n][r];
	return 0;
} 

相信读者也发现了，递推式本身过于消耗空间，因此递推法能解决的范围也很有限。当然细心的读者应该发现了，这边的递推数组是一个下三角矩阵，因此可以采用坐标转换的方式化为一维存储，可以省去一半的空间。
注意，这个方法需要在打表计算的同时就进行取模操作$(a+b)\mathrm{m}\mathrm{o}\mathrm{d}p=(a\mathrm{m}\mathrm{o}\mathrm{d}p+b\mathrm{m}\mathrm{o}\mathrm{d}p)$，因为在不取模进行累加时，数据的增长是异常迅速的，因此很快就会数据溢出报错。

思路五

用数论的知识，所有的数都可以使用素因数进行分解：即对任意一个数n，存在一个素数序列$P_i$，使得$$n=\prod _{i=1}^{k}P{i}^{Qi}$$我们使用质因数分解的办法，将除数和被除数分解成小的素数因子后，再统计计算结果。因而我们需要一个素数表，我们可以使用欧拉筛快速获得这张表。而也正是因为这样，我们的算法也有了局限性，由于内存的限制，我们没办法获得很大的素数，也就导致了当需要被分解的数中存在大素数因子的时候，我们就束手无测。（或者当素数表使用完之后，继续向后枚举素数，这样会降低算法性能），不过一般情况下，至少处理百万级的因数时，是完全可以的。

#include<bits/stdc++.h>
using namespace std;
const int MaxN = 1000010;
bool isPrime[MaxN];
int Primes[MaxN],cnt,Factor[MaxN];
void EulrSift()//欧拉筛
{
	int i,j;
	for(i=2;i<=MaxN;i++)
		if(!isPrime[i]){
			Primes[cnt++]=i;
			for(j=0;j<cnt&&Primes[j]*i<MaxN;++j)
			{
				isPrime[Primes[j]*i]=true; 
				if(i%Primes[j]) break;
			}
		}
}
void PrimeResolve(int num,int delta)//质因数分解
{
	for(int i=0;num!=1&&i<cnt;i++)
		while(num%Primes[i]==0&&num!=1)
		{
			num/=Primes[i];
			Factor[i]+=delta;
		}
}
int quickMul(int a,int b)//快速幂
{
	int r = 1;
	while(b)
	{
		if(b&1) r*=a;
		b>>=1;
		a*=a;
	}
	return r;
}
int main()
{
	int n,r,m,i,ans=1;
	cin>>n>>r;
	m=n;
	EulrSift();
	for(i=1;i<=r;i++,n--)
	{
		PrimeResolve(n,1);
		PrimeResolve(i,-1);
	}
	for(i=0;i<cnt&&Primes[i]<=m;i++)
		if(Factor[i]!=0 )
			ans*=quickMul(Primes[i],Factor[i]);
	cout<<ans;
	return 0;
}

思路六

从这里开始，我们会接触到求大组合数+取模的问题。观察组合数的运算公式，我们发现运算过程中有除法，但是在取模条件下，我们可以把除法使用逆元的方法转化为乘法来做。费马定理的要求过于苛刻，我们可以采用扩展欧几里得来求。欧拉定理比费马小定理的通用性更强，但是有时候模数过大，可能导致欧拉定理无法使用，另外，对于欧拉定理而言，一般我们不使用他的关键原因在于，本身欧拉定理告诉我们，如果$\gcd (A,P)=1$，那么一定会有：$$A^{\phi (P)}=1(\mathrm{m}\mathrm{o}\mathrm{d}P)$$但是即便我们使用欧拉筛取求解这个$\phi (P)$，也需要$O(\sqrt{P})$的复杂度，并不值得。但是我们还是复习一遍欧拉函数吧(233)。
我们针对问题：$C_n^r\mathrm{m}\mathrm{o}\mathrm{d}p$，给出下面两种解决办法：
扩展欧几里得版：

#include<bits/stdc++.h>
using namespace std;
int extgcd(int A,int B,int &x,int &y)
{
	if(B==0)
	{
		x=1;
		y=0;
		return A; 
	}
	int gcd=extgcd(B,A%B,x,y);
	int x0=y;
	int y0=x-A/B*y;
	x=x0;
	y=y0;
	return gcd;
}
//A*x=1(mod P)=>A*x+P*y=1
int getInv(int A,int p)
{
	int x,y,gcd=extgcd(A,p,x,y);
	if(gcd!=1) return -1;//扩展欧几里得无法求出逆元 
	p/=gcd;
	return (x%p+p)%p;//取最小正整数解 
}
int main()
{
	int i,n,r,p,ans=1,inv;
	cin>>n>>r>>p;
	for(i=1;i<=r;i++,n--)
	{
		ans=ans*n%p;
		if(ans%r==0) ans/=r;//如果可以除掉,就不要求逆元 
		else if(r>=2)//1无需求逆元
		{
			inv=getInv(r,p);
			if(inv==-1)
			{
				cout<<"No solution"<<endl;
				return 0;
			}
			ans=ans*inv%p;
		}
	}
	cout<<ans;
	return 0;
}

欧拉函数版：

#include<bits/stdc++.h>
using namespace std;
int cnt,phi[100010],Primes[100010];
int Phi(int n)//质因数分解 结合 欧拉函数性质
{ 
	int res=n;
	for(int i=2;i*i<=n;i++)
	{ 
		if(n%i==0) res=res/i*(i-1);
		while( n % i == 0 ) n/=i;
	}
	if(n>1)//除干净了因子后 n = 1 或者是最后一个质因数的幂次 
	{
		res=res/n*(n-1);
	}
	return res;
}
int quickMul(int A,int B,int P)
{
	int r=1;
	while(B)
	{
		if(B&1) r=r*A%P;
		B>>=1;
		A=A*A%P;
	}
	return r;
}
int gcd(int A,int B)
{
	while(B^=A^=B^=A%=B);
	return A;
}
int main()
{
	int n,r,p,i,ans=1,phi,inv;//son用来存储那些没有逆元的因子 
	cin>>n>>r>>p;
	phi=Phi(p);
	for(i=1;i<=r;i++,n--)
	{
		ans=ans*n ;
		if(gcd(i,p)!= 1)
		{
			cout<<"No solution"<<endl; 
			return 0;
		}
		else
		{
			inv=quickMul(i,phi,p);
			ans=ans*inv%p;
		}
		ans%=p;
	}
	cout<<ans;
	return 0;
} 

有些读者可能会想到，本身求逆元是需要时间的，能否在能整除时做除，不能做整除的时候求逆元改为乘法，同时每一步保持取模，答案是否定的，我们不能在乘法、除法、逆元分别取模后再合并。
以$\frac{12}{3}(\mathrm{m}\mathrm{o}\mathrm{d}9)=4$为例，如果我们做完分子的取模后，再除掉分母后再取模，那么答案为：$$\frac{12(\mathrm{m}\mathrm{o}\mathrm{d}9)}{3}(\mathrm{m}\mathrm{o}\mathrm{d}9)=\frac{3}{3}(\mathrm{m}\mathrm{o}\mathrm{d}9)=1$$原因在于，取模本身不满足和除法得交换律：

1. 因此，我们只能选择，在可以整除的时候，整除，但是一旦除法参与了，那么取模运算必须放到最后进行。
2. 全部改用逆元做乘法。

我们使用这些算法的目的就是为了避免乘法溢出，因此方案1肯定是不切实际的。

思路7：卢卡斯定理

卢卡斯定理提供了一种在模数为素数的情况下，由于n和r远大于p，所以在除法过程中求逆元的时候很容易遇到与p不满足互质的情况，例如2p，3p，从而由于不满足互质导致无法求解逆元，但是卢卡斯定理提供了一种这种情形下还能求解大组合数的方案。
定理内容：对于非负整数n,r,p(其中p为素数)，那么有如下的同余方程：$$C_n^r=\prod _{i=0}^k{C}_{n_i}^{r_i}(\mathrm{m}\mathrm{o}\mathrm{d}p)$$其中，$n_i$,$r_i$需要满足如下等式：$$n_k{p}^k+n_{k-1}{p}^{k-1}+...+n_{1}p+n_0$$$$r_k{p}^k+r_{k-1}{p}^{k-1}+...+r_{1}p+r_0$$不难发现，这就是分别将n,r做p进制下的分解罢了。
在实际的算法实现上，我们会这么使用这个定理，与上述是等价的：$$Luca{s}_n^r=C_{n\mathrm{m}\mathrm{o}\mathrm{d}p}^{r\mathrm{m}\mathrm{o}\mathrm{d}p}Luca{s}_{n/p}^{r/p}$$注意上述公式中的除法，是整除，其实上述公式完全等价于卢卡斯定理原本的公式。
需要注意的是，在这个定理中，如果发生了$n<r$的情况，那么规定$C_n^r=0$。

公式证明：求解$C_p^x(\mathrm{m}\mathrm{o}\mathrm{d}p)$，设$x<p$，$x\ne 0$$$C_p^x=\frac{p!}{x!(p-x)!}=\frac{p\ast (p-1)!}{x\ast (x-1)!\ast (p-x)!}=\frac{p}{x}\ast C_{p-1}^{x-1}$$由于我们已经设$x<p$且p为素数，因此我们知道x肯定有关于p的逆元，所以上述式子可以转换为：$$C_p^x(\mathrm{m}\mathrm{o}\mathrm{d}p)=\frac{p}{x}\ast C_{p-1}^{x-1}(\mathrm{m}\mathrm{o}\mathrm{d}p)=p\ast inv(x\mid p)\ast C_{n-1}^{r-1}(\mathrm{m}\mathrm{o}\mathrm{d}p)$$$inv(x\mid p)$表示x在模p下的逆元。
此时等式右端是p的倍数，因此$C_p^x=0(\mathrm{m}\mathrm{o}\mathrm{d}p)$。
根据二项式定理：$$(1+x{)}^n=\sum _{i=0}^n{C}_n^i\ast x^i$$我们知道：$$(1+x{)}^p=\sum _{i=0}^p{C}_p^i\ast x^i$$根据我们上面所证，这个上述式子的展开式中，只有$i=0$和$i=p$的两项模p不为0。故我们得到：$$(1+x{)}^p=1+x^p(\mathrm{m}\mathrm{o}\mathrm{d}p)$$我们先预设：$$\{\begin{array}{l}n/p=n_k\\ n\mathrm{m}\mathrm{o}\mathrm{d}p=re{s}_k\end{array}$$我们对$(1+x{)}^n$做p进制分解：
$$\begin{array}{rl}(1+x{)}^n& =(1+x{)}^{n_k\ast p+re{s}_k}\\ & =(1+x{)}^{n_k\ast p}+(1+x{)}^{re{s}_k}\\ & =[(1+x{)}^p{]}^{n_k}\ast (1+x{)}^{re{s}_k}......................(1)\\ & =(1+x^p{)}^{n_k}\ast (1+x{)}^{re{s}_k}........................(2)\\ & =\sum _{i=0}^{n_k}{C}_{n_k}^i\ast (x^p{)}^i\ast \sum _{j=0}^{re{s}_k}\ast x^j.......................(3)\\ & =\sum _{i=0}^{n_k}\sum _{j=0}^{re{s}_k}{C}_{n_k}^i\ast C_{re{s}_k}^j\ast x^{ip+j}(\mathrm{m}\mathrm{o}\mathrm{d}p)....(4)\end{array}$$(1)->(2)：在%P的条件下，我们用了上述推出的结论$(1+x{)}^p=1+x^p(\mathrm{m}\mathrm{o}\mathrm{d}p)$
(2)->(3)：使用两次二项式展开定理
(3)->(4)：多项式的乘法就是分别乘法后求和
注意(2)式中我们用了上面推导的结论：$(1+x{)}^p=1+x^p(\mathrm{m}\mathrm{o}\mathrm{d}p)$
所以我们得到了： $$\begin{array}{rl}(1+x{)}^p& =\sum _{r=0}^n{C}_n^r\ast x^r\\ & =\sum _{i=0}^{n_k}{C}_{n_k}^i{x}^{ip}\ast \sum _{j=0}^{re{s}_k}{C}_{re{s}_k}^j{x}^j\\ & =\sum _{i=0}^{n_k}\sum _{j=0}^{re{s}_k}{C}_{n_k}^i{C}_{re{s}_k}^j{x}^{ip+j}(\mathrm{m}\mathrm{o}\mathrm{d}p),(x<p,re{s}_k<p)\end{array}$$由于上述式是恒等式，因此比较两端$x^k$的系数可知：$$C_n^r{x}^r=C_{n_k}^i{C}_{re{s}_k}^j{x}^{ip+j}(\mathrm{m}\mathrm{o}\mathrm{d}p),(x<p,re{s}_k<p)$$即：$$C_n^r=C_{n_k}^i{C}_{re{s}_k}^j(\mathrm{m}\mathrm{o}\mathrm{d}p),(x<p,re{s}_k<p)$$$$n=p{n}_k+re{s}_k$$$$r=ip+j$$则$$C_n^r=C_{n_k}^i{C}_{re{s}_k}^j=C_{n/p}^{r/p}{C}_{n\mathrm{m}\mathrm{o}\mathrm{d}p}^{r\mathrm{m}\mathrm{o}\mathrm{d}p}$$证毕。由于$r\mathrm{m}\mathrm{o}\mathrm{d}p<p$，$n\mathrm{m}\mathrm{o}\mathrm{d}p<p$，因此对于$C_{n\mathrm{m}\mathrm{o}\mathrm{d}p}^{r\mathrm{m}\mathrm{o}\mathrm{d}p}$我们就可以比较容易地计算（因为p是质数，所以对与$x<p$,一定有$\gcd (x,p)=1$，互质条件下可以求逆元）
在卢卡斯定理的加持下，我们可以保证在如下式子中：$$C_n^r=\prod _{i=0}^k{C}_{n_i}^{r_i}(\mathrm{m}\mathrm{o}\mathrm{d}p)$$一定会有$n_i,r_i<p$(按照p进制数理解，每一位上的权值不可能大于它所表示的数的进制值)。
而又由于p为素数，因此我们又可以用逆元啦，完美！

#include<bits/stdc++.h>
using namespace std;
int quickMul(int A,int B,int P)
{
	int res=1;
	while(B)
	{
		if(B&1) res=res*A%P;
		B>>=1;
		A=A*A%P;
	}
	return res;
}
int inv(int n,int p)// 费马小定理求解逆元  n^(p-1)=1(mod p) 
{
	return quickMul(n,p-2,p);
}
int C(int n,int r,int p)
{
	int res=1;
	for(int i=1;i<=r;i++,n--) res=res*n%p*inv(i,p)%p;
	return res;
}
int Lucas(int n,int r,int p)
{
	return r==0?1:Lucas(n/p,r/p,p)*C(n%p,r%p,p)%p; 
}
int main()
{
	int n,r,p;
	cin>>n>>r>>p;
	if(n<r||n==0)
	{
		cout<<-1;// 视具体题目要求而定，毕竟数学上规定 0!=1 
		return 0;
	}
	if(p==1)
	{
		cout<<0; 
		return 0;
	}
	cout<<Lucas(n,r,p);
	return 0;
}

来看一道模板题吧~

洛谷 P3807 【模板】卢卡斯定理

题目背景

这是一道模板题。

题目描述

给定整数 $n,m,p$ 的值，求出 $C_{n+m}^n\mathrm{m}\mathrm{o}\mathrm{d}p$ 的值。
输入数据保证 $p$为质数。
注：$C$ 表示组合数。

输入格式

本题有多组数据。
第一行一个整数 $T$，表示数据组数。
对于每组数据:
一行，三个整数 $n$, $m$, $p$。

输出格式

对于每组数据，输出一行，一个整数，表示所求的值。

输入输出样例

输入#1

2
1 2 5
2 1 5

输出#1

3
3

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
ll quickMul(ll A,ll B,ll P)
{
	ll res = 1;
	while(B)
	{
		if(B&1) res=res*A%P;
		B>>=1;
		A=A*A%P;
	}
	return res;
}
 
ll inv(ll n,ll p)// 费马小定理求解逆元  n^(p-1)=1(mod p) 
{
	return quickMul(n,p-2,p);
}
ll C(ll n,ll r,ll p)
{
	ll res=1;
	for(int i=1;i<=r;++i,--n) res=res*n%p*inv(i,p)%p;
	return res;
}
ll Lucas(ll n,ll r,ll p)
{
	return r==0?1:Lucas(n/p,r/p,p)*C(n%p,r%p,p)%p; 
}
int main()
{
	ll T,n,m,p;
	cin>>T;
	while(T--)
	{
		cin>>n>>m>>p;
		n+=m;
		if(p==1) cout<<0<<endl;
		else cout<<Lucas(n,m,p)<<endl;
	}
	return 0;
}

奇技淫巧：

在使用卢卡斯定理的时候，我们还可以预处理阶乘，这样可以加快速度，在预处理阶乘时，我们采用公式二：$$C_n^r=\frac{n!}{r!(n-r)!}$$事实上我们可以得到：$$\begin{array}{rl}{C}_n^r(\mathrm{m}\mathrm{o}\mathrm{d}p)& =\frac{n!}{r!(n-r)!}(\mathrm{m}\mathrm{o}\mathrm{d}p)\\ & =n!(\mathrm{m}\mathrm{o}\mathrm{d}p)\ast inv(r!\mid p)(\mathrm{m}\mathrm{o}\mathrm{d}p)\ast inv(n-r!(\mathrm{m}\mathrm{o}\mathrm{d}p)\mid p)(\mathrm{m}\mathrm{o}\mathrm{d}p)\\ & =n!(\mathrm{m}\mathrm{o}\mathrm{d}p)\ast inv(r!(\mathrm{m}\mathrm{o}\mathrm{d}p)\mid p)(\mathrm{m}\mathrm{o}\mathrm{d}p)\ast inv(n-r!(\mathrm{m}\mathrm{o}\mathrm{d}p)\mid p)(\mathrm{m}\mathrm{o}\mathrm{d}p)\end{array}$$即我们需要证明$\frac{x!}{y!}(\mathrm{m}\mathrm{o}\mathrm{d}p)=x!\ast inv(y!(\mathrm{m}\mathrm{o}\mathrm{d}p)\mid p)(\mathrm{m}\mathrm{o}\mathrm{d}p)$
设$z=y!(\mathrm{m}\mathrm{o}\mathrm{d}p)$，即$y=z+kp$，那么有$z<p$，且p为素数，故z必有逆元，则有：$$z\ast inv(z\mid p)=1(\mathrm{m}\mathrm{o}\mathrm{d}p)$$那么$$\begin{array}{rl}y!\ast inv(z)& =(kp+z)inv(z\mid p)\\ & =z\ast inv(z)+kp\ast inv(z\mid p)\\ & =z\ast inv(z)(\mathrm{m}\mathrm{o}\mathrm{d}p)\\ & =1\end{array}$$证毕。

#include<bits/stdc++.h>
using namespace std;
const int MaxN=1e5+10;
typedef long long ll;
ll Factorial[MaxN<<1];
ll quickMul(ll A,ll B,ll P)
{
	ll res=1;
	while(B)
	{
		if(B&1) res=res*A%P;
		B>>=1;
		A=A*A%P;
	}
	return res;
}
ll inv(ll n,ll p)
{
	return quickMul(n,p-2,p);
}
ll C(ll n,ll r,ll p)// 这里需要注意 n < r 时 是非法状态 返回 0 
{
	return n<r?0:Factorial[n]*inv(Factorial[r],p)%p*inv(Factorial[n-r],p)%p; 
}
ll Lucas(ll n,ll r,ll p)
{
	return r==0?1:Lucas(n/p,r/p,p)*C(n%p,r%p,p)%p; 
}
int main()
{
	ll T,n,m,p,i;
	cin>>T;
	Factorial[1]=1;
	while(T--)
	{
		cin>>n>>m>>p;
		n+=m;
		memset(Factorial,0,sizeof(Factorial));
		Factorial[1]=1;
		if(p==1) cout<<0<<endl;
		else
		{
			m=m>(n>>1)?n-m:m;
			for(i=2;i<=n;i++) Factorial[i]=Factorial[i-1]*i%p;
			cout<<Lucas(n,m,p)<<endl;
		}
	}
	return 0;
}

————————————————
版权声明：本文为CSDN博主「TMZT」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
原文链接：https://blog.csdn.net/qq_39304630/article/details/77417017
有更改

思路八：扩展Lucas定理

在Lucas定理中，模数只能是质数。但很多情况下，模数不一定是质数。这个时候就需要用到Lucas定理的扩展版——扩展Lucas定理。
扩展卢卡斯定理用于求如下式子（其中p不一定是质数）：$$C_m^n\mathrm{m}\mathrm{o}\mathrm{d}p$$
我们将这个问题由总体到局部地分为三个层次解决。

层次一：原问题

首先对p进行质因数分解：$$p=\prod _i{p}_i^{k_i}$$显然 $p_i^{k_i}$ 是两两互质的，所以如果分别求出 $C_m^n\mathrm{m}\mathrm{o}\mathrm{d}{p}_i^{k_i}$，就可以构造出若干个形如 $C_n^m=a_i\mathrm{m}\mathrm{o}\mathrm{d}{p}_i^{k_i}$的方程，然后用中国剩余定理即可求解。

层次二：组合数模质数幂

现在的问题就转化成了求如下式子（其中p是质数）：$$C_n^m\mathrm{m}\mathrm{o}\mathrm{d}{p}^k$$脑补一下组合数公式 $C_n^m=\frac{n!}{m!(n-m)!}$，发现由于 $m!$ 和 $(n-m)!$ 可能包含质因子p，所以不能直接求他们对于 $p^k$ 的逆元。此时我们可以将 $n!$ 、$m!$、$(n-m)!$ 中的质因子 $p$ 全部提出来，最后再乘回去即可。即变为下式（$k1$ 为 $n!$ 中质因子 $p$ 的次数，$k2$、$k3$同理）：$$\frac{\frac{n!}{p^{k1}}}{\frac{m!}{p^{k2}}\times \frac{(n-m)!}{p^{k3}}}\times p^{p1-p2-p3}$$$\frac{m!}{p^{k2}}$ 和 $\frac{(n-m)!}{p^{k3}}$ 和 $p^k$ 是互质的，可以直接求逆元。

层次三：阶乘除去质因子后模质数幂

现在看看如何计算形如下式的式子。$$\frac{n!}{p^a}\mathrm{m}\mathrm{o}\mathrm{d}{p}^k$$先考虑如何计算 $n!\mathrm{m}\mathrm{o}\mathrm{d}{p}^k$。
举个例子：n=22，p=3，k=2
把这个写出来：
$22!=1\times 2\times 3\times 4\times 5\times 6\times 7\times 8\times 9\times 10\times 11\times 12\times 13\times 14\times 15\times 16\times 17\times 18\times 19\times 20\times 21\times 22$
把其中所有p（也就是3）的倍数提取出来，得到：
$22!=3^7\times (1\times 2\times 3\times 4\times 5\times 6\times 7)\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\times 20\times 22)$
可以看出上式分为三个部分：第一个部分是3的幂，次数是小于等于22的3的倍数的个数，即 $\lfloor \frac{n}{p}\rfloor$ 。
第二个部分是一个阶乘7!，即 $\lfloor \frac{n}{p}\rfloor !$，可以递归解决。
第三个部分是n!中与p互质的部分的乘积，这一部分具有如下性质：
$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}{p}^k)$
在模 $3^2$ 的意义下10和1同余，11和2同余……写成下式就比较显然。（t是任意正整数）$$\prod _{i,(i,p)=1}^{p^k}i\equiv \prod _{i,(i,p)=1}^{p^k}(i+t{p}^k)(\mathrm{m}\mathrm{o}\mathrm{d}{p}^k)$$${\prod }_{i,(i,p)=1}^{p^k}i$ 一共循环了 $\lfloor \frac{n}{p^k}\rfloor$ 次，暴力求出 ${\prod }_{i,(i,p)=1}^{p^k}i$ ，然后用快速幂求它的 ${\prod }_{i,(i,p)=1}^{p^k}i$ 次幂。
最后还要乘上19×20×22（即 ${\prod }_{i,(i,p)=1}^{n\mathrm{m}\mathrm{o}\mathrm{d}{p}^k}$），显然这一段的长度一定小于 $p^k$ ，暴力乘上去即可。
如上三部分的乘积就是 $n!$。最终要求的是 $\frac{n!}{p^a}\mathrm{m}\mathrm{o}\mathrm{d}{p}^k$ ，分母全部由上述第一部分和第二部分贡献（第三部分和p互质）。而递归计算第二部分的时候已经除去了第二部分中的因子p，所以最终的答案就是上述第二部分递归返回的结果和第三部分的乘积（与第一部分无关）。
结合代码方便理解：

ll fac(const ll n,const ll p,const ll pk)
{
	if (!n)
		return 1;
	ll ans=1;
	for(int i=1;i<pk;i++)
		if(i%p)
			ans=ans*i%pk;
	ans=power(ans,n/pk,pk);
	for(int i=1;i<=n%pk;i++)
		if(i%p)
			ans=ans*i%pk;
	return ans*fac(n/p,p,pk)%pk;
}

层次二：组合数模质数幂

回到这个式子$$\frac{\frac{n!}{p^{k1}}}{\frac{m!}{p^{k2}}\times \frac{(n-m)!}{p^{k3}}}\times p^{p1-p2-p3}$$可以很容易地把它转换成代码（注意i要开long long）：

ll C(const ll n,const ll m,const ll p,const ll pk)
{
	if(n<m) return 0;
	ll f1=fac(n,p,pk),f2=fac(m,p,pk),f3=fac(n-m,p,pk),cnt=0;
	for(ll i=n;i;i/=p)
		cnt+=i/p;
	for(ll i=m;i;i/=p)
		cnt-=i/p;
	for (ll i=n-m;i;i/=p)
		cnt-=i/p;
	return f1*inv(f2,pk)%pk*inv(f3,pk)%pk*power(p,cnt,pk)%pk;
}

层次一：原问题

完整代码（题目：洛谷4720【模板】扩展卢卡斯）：

#include<bits/stdc++.h>
using namespace std;
const int N = 1e6;
typedef long long ll;
ll n, m, p;
ll power(ll a,ll b,ll p=0x7fffffffffffffff)
{
	ll ans=1;
	while(b)
	{
		if(b&1) ans=ans*a%p;
		a=a*a%p;
		b>>=1;
	}
	return ans;
}
ll fac(ll n,ll p,ll pk)
{
	if(!n) return 1;
	ll ans=1;
	for(int i=1;i<pk;i++) if(i%p) ans=ans*i%pk;
	ans=power(ans,n/pk,pk);
	for(int i=1;i<=n%pk;i++) if(i%p) ans=ans*i%pk;
	return ans*fac(n/p,p,pk)%pk;
}
ll exgcd(ll a,ll b,ll &x,ll &y)
{
	if(!b)
	{
		x=1;
		y=0;
		return a;
	}
	ll xx,yy,g=exgcd(b,a%b,xx,yy);
	x=yy;
	y=xx-a/b*yy;
	return g;
}
ll inv(const ll a, const ll p)
{
	ll x,y;
	exgcd(a,p,x,y);
	return (x%p+p)%p;
}
ll C(ll n,ll m,ll p,ll pk)
{
	if(n<m) return 0;
	ll f1=fac(n,p,pk),f2=fac(m,p,pk),f3=fac(n-m,p,pk),cnt=0;
	for(ll i=n;i;i/=p) cnt+=i/p;
	for(ll i=m;i;i/=p) cnt-=i/p;
	for(ll i=n-m;i;i/=p) cnt-=i/p;
	return f1*inv(f2,pk)%pk*inv(f3,pk)%pk*power(p,cnt,pk)%pk;
}
ll a[N], c[N],cnt;
ll crt()
{
	ll M=1,ans=0;
	for(int i=0;i<cnt;i++) M*=c[i];
	for(int i=0;i<cnt;i++) ans=(ans+a[i]*(M/c[i])%M*inv(M/c[i],c[i])%M)%M;
	return ans;
}
ll exlucas(ll n,ll m,ll p)
{
	ll tmp=sqrt(p);
	for(int i=2;p>1&&i<=tmp;i++)
	{
		ll tmp=1;
		while(p%i==0)
		{
			p/=i;
			tmp*=i;
		}
		if(tmp>1)
		{
			a[cnt]=C(n,m,i,tmp);
			c[cnt++]=tmp;
		}
	}
	if(p>1)
	{
		a[cnt]=C(n,m,p,p);
		c[cnt++]=p;
	}
	return crt();
}
int main()
{
	cin>>n>>m>>p;
	cout<<exlucas(n,m,p);
	return 0;
}

————————————————
版权声明：本文为CSDN博主「Inspector_Javert」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
原文链接：https://blog.csdn.net/hqddm1253679098/article/details/82897638
有改动