卢卡斯定理

当 (p) 为素数时

[inom{n}{m} mod p=inom{lfloor n/p floor}{lfloor m/p floor}inom{n mod p}{m mod 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! mod p_i^{c_i})。

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

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

其中，(a) 是 (3^6)，(b) 是 (1 imes 2 imes 3 imes 4 imes 5 imes 6)，(c) 是 (1 imes 2 imes 4 imes 5 imes 7 imes 8 imes 10 imes 11 imes 13 imes 14 imes 16 imes 17 imes 19)。

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

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

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

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

(n!) 中 (p) 的个数显然为 (lfloor n/p floor + lfloor n/p^2 floor + 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;
}