卢卡斯定理

当 \(p\) 为素数时

\[\binom{n}{m} \bmod p=\binom{\lfloor n/p \rfloor}{\lfloor m/p \rfloor}\binom{n \bmod p}{m \bmod p} \]

背会就好
luogu 模板

#include <iostream>
#include <cstdio>
using namespace std;
typedef long long ll;
int n, m, p, jie[200005], ni[200005], T;
int lucas(int a, int b, int p){
	if(a<b)	return 0;
	else if(a<p)	return (ll)jie[a]*ni[b]*ni[a-b]%p;
	else	return (ll)lucas(a/p,b/p,p)*lucas(a%p,b%p,p)%p;
}
int main(){
	cin>>T;
	while(T--){
		cin>>n>>m>>p;
		jie[0] = jie[1] = ni[0] = ni[1] = 1;
		for(int i=2; i<=n+m; i++)
			jie[i] = ((ll)jie[i-1] * i) % p;
		for(int i=2; i<=n+m; i++)
			ni[i] = (ll)(p-p/i) * ni[p%i] % p;
		for(int i=2; i<=n+m; i++)
			ni[i] = (ll)ni[i-1] * ni[i] % p;
		cout<<lucas(n+m, m, p)<<endl;
	}
	return 0;
}

当 \(p\) 没有限制时

参考zyf2000。

将 \(p\) 分解成 \(p_1^{c_1}p_2^{c_2}\cdots p_k^{c_k}\)。
根据每个 \(p_i^{c_i}\) 求出答案然后用中国剩余定理合并。

然而 \(p_i^{c_i}\) 并不是素数。

先来研究研究求 \(n! \bmod p_i^{c_i}\)。

比方说， \(19! \bmod 3^2\)。

\(19! \equiv abc \pmod p_i^{c_i}\)

其中，\(a\) 是 \(3^6\)，\(b\) 是 \(1 \times 2 \times 3 \times 4 \times 5 \times 6\)，\(c\) 是 \(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\)。

因此，求解 \(n! \bmod p_i^{c_i}\) 可以分为三个部分：

• \(p_i^{\lfloor n/p_i \rfloor}\)
• \(\lfloor n/p_i \rfloor !\)
• 其余的数

我们发现， 第三部分在模意义下是以 \(p_i^{c_i}\) 为长度循环的（这里的长度不是数字个数，而是权值）。即 \(1 \times 2 \times 4 \times 5 \times 7 \times 8 \equiv 10 \times 11 \times 13 \times 14 \times 16 \times 17 \pmod {p_i^{c_i}}\)。因此求出一个周期的然后快速幂，最后剩下几个暴力算就行了。

我们还发现，这样 \(p_i^{c_i}\) 和 \(n!\) 不一定互素。那么，我们就先解决掉 \(n!\) 中的 \(p_i\) 因子，算阶乘的时候就不要算第一部分了，最后再乘上。

\(n!\) 中 \(p\) 的个数显然为 \(\lfloor n/p \rfloor + \lfloor n/p^2 \rfloor + \cdots\)

cf 2015ICL,Finals,Div.1#J Ceizenpok’s formula

#include <iostream>
#include <cstdio>
using namespace std;
typedef long long ll;
ll ans;
ll n, m, mod, x, y;
void exgcd(ll a, ll b){
	if(!b){
		x = 1;
		y = 0;
		return ;
	}
	exgcd(b, a%b);
	ll z=x;
	x = y;
	y = z - a / b * y;
}
ll inv(ll a, ll p){
	if(!a)	return 0;
	exgcd(a, p);
	ll re=(x%p+p)%p;
	if(!re)	re += p;
	return re;
}
ll ksm(ll a, ll b, ll p){
	ll re=1;
	while(b){
		if(b&1)	re = (re * a) % p;
		a = (a * a) % p;
		b >>= 1;
	}
	return re;
}
ll fac(ll x, ll pi, ll pk){
	if(!x)	return 1;
	ll re=1;
	if(x/pk){
		for(ll i=2; i<=pk; i++)
			if(i%pi)
				re = re * i % pk;
		re = ksm(re, x/pk, pk);
	}
	for(ll i=2; i<=x%pk; i++)
		if(i%pi)
			re = re * i % pk;
	return re*fac(x/pi,pi,pk)%pk;
}
ll C(ll n, ll m, ll mod, ll pi, ll pk){
	if(n<m)	return 0;
	ll a=fac(n,pi,pk), b=fac(m,pi,pk), c=fac(n-m,pi,pk);
	ll zhi=0;
	for(ll i=n; i; i/=pi)	zhi += i / pi;
	for(ll i=m; i; i/=pi)	zhi -= i / pi;
	for(ll i=n-m; i; i/=pi)	zhi -= i / pi;
	return a*inv(b,pk)%pk*inv(c,pk)%pk*ksm(pi,zhi,pk)%pk;
}
int main(){
	cin>>n>>m>>mod;
	ll tmp=mod;
	for(ll i=2; i<=mod; i++)
		if(tmp%i==0){
			ll pk=1;
			while(tmp%i==0)	pk *= i, tmp /= i;
			ans = (ans + C(n,m,mod,i,pk)*(mod/pk)%mod*inv(mod/pk,pk)%mod) % mod;
		}
	cout<<ans<<endl;
	return 0;
}