当 (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;
}