快速阶乘与（扩展）卢卡斯定理
\(p\)为质数时
考虑 \(n!~mod~p\) 的性质
当\(n>>p\)时，不妨将\(n!\)中的因子\(p\)提出来
\(n!\) 可以写成 \(a*p^e\) ， \(a\)与\(p\)互质
如何求解\(a\)和\(e\)？
显然，\(e=n/p+n/p^2+n/p^3+……\)
因为\(1\)~\(n\)有\(n/p\)个\(p\)的倍数，贡献为\(1\)，\(n/p^2\)个\(p^2\)的倍数，贡献为\(2\)……
事实上，可以每次先将\(1\)~\(n\)中\(p\)的倍数的因子提出来，\(p,2p,3p...(n/p)p\)就变成了\(1,2,3...n/p\)，\(e+=n/p\),而\(1\)~\(n/p\)这些数的\(a\)和\(e\)又可以递归求解
再考虑\(a\)的求法：
将\(1\)~\(n\)中\(p\)的倍数拿走，递归处理求得\(1\)~ \(n/p\)的\(a\)，再乘上剩下的数\(1,2,3...p-1,p+1,p+2,...2p-1,2p+1,2p+2,...\)就行了
而剩下的数在对\(p\)取模意义下是从\(1\)~\(p-1\)的循环，于是可以预处理\(n< p\)时的\(n!\)
剩下的数的积即为: \(((p-1)!)^{n/p}\times(n~mod~p)!\)
根据威尔逊定理，\((p-1)!\equiv 1 ~(mod~p)\)，可以省去快速幂
递归求解即可
#include<iostream>
#include<cstring>
#include<cstdio>
#define int long long
using namespace std;
const int MAXN=100010;
int fact[MAXN];
inline int qpow(int x,int k,int p){
int s=1%p;
while(k){
if(k&1) s=s*x%p;
k>>=1;
x=x*x%p;
}
return s;
}
inline int mod_fact(int n,int &e,int p){
e=0;
if(n<p) return fact[n];
int res=mod_fact(n/p,e,p);
e+=n/p;
if(n/p%2!=0) return res*(p-fact[n%p])%p; //((p-1)!)^(n/p)=(-1)^(n/p)
return res*fact[n%p]%p;
}
inline int C(int n,int m,int p){
if(m>n) return 0;
int e1,e2,e3;
int a1=mod_fact(n,e1,p),a2=mod_fact(m,e2,p),a3=mod_fact(n-m,e3,p);
if(e1>e2+e3) return 0; //p^(e1-e2-e3)是p的倍数,return 0
return a1*qpow(a2*a3%p,p-2,p)%p;
}
int n,m,p;
signed main()
{
int T;
scanf("%lld",&T);
while(T--){
scanf("%lld%lld%lld",&n,&m,&p);
fact[0]=1;
for(int i=1;i<=p;++i) //预处理阶乘
fact[i]=fact[i-1]*i%p;
printf("%lld\n",C(n+m,m,p));
}
return 0;
}\(p\)为任意数时
考虑将\(p\)分解为\(p_1^{a_1}p_2^{a_2}p_3^{a_3}...p_k^{a_k}\)
分别求解\(C(n,m)~mod ~ p_i^{a_i}\)
然后用\(CRT\)合并就可以了
求解\(C(n,m)~mod~p_i^{a_i}\)时，仍然递归求解\(1\)~\(n/p_i\)，但是\(a\)是对\(p_i^{a_i}\)取模，而非对\(p_i\)取模
非\(p_i\)的倍数的数乘积计算与上面也有不同，首先\(p_i^{a_i}\)是合数，逆元只能用\(exgcd\)求了
因为是对\(p_i^{a_i}\)取模,循环是\(p_i^{a_i}-a_i\)，而且不能用威尔逊定理，需要先从\(1\)~\(p_i^{a_i}\)暴力乘，再用快速幂直接计算，剩下的再暴力乘进去即可，此时预处理\(fact[i]\)就不必要了
代码:
#include<iostream>
#include<cstring>
#include<cstdio>
#include<cmath>
#define int long long
using namespace std;
const int MAXN=1000010;
inline int qpow(int x,int k,int p){
int s=1;
while(k){
if(k&1) s=s*x%p;
k>>=1;
x=x*x%p;
}
return s;
}
inline int exgcd(int a,int b,int &x,int &y){
if(!b){
x=1; y=0;
return a;
}
int t=exgcd(b,a%b,y,x);
y-=a/b*x;
return t;
}
inline int mod_fact(int n,int &e,int p,int MOD){
if(!n) return 1;
int res=1;
for(int i=2;i<=MOD;++i) //暴力累乘一个循环
if(i%p) res=res*i%MOD;
res=qpow(res,n/MOD,MOD); //共有n/MOD个循环
e+=n/p;
for(int i=2;i<=n%MOD;++i)
if(i%p) res=res*i%MOD; //乘上循环外面的数
return res*mod_fact(n/p,e,p,MOD)%MOD; //递归处理n/p个p的倍数
}
inline int C(int n,int m,int p,int MOD){
if(m>n) return 0;
int e1=0,e2=0,e3=0;
int a1=mod_fact(n,e1,p,MOD),a2=mod_fact(m,e2,p,MOD),a3=mod_fact(n-m,e3,p,MOD);
if(e1-e2-e3>0){
int temp=pow(p,e1-e2-e3);
if(temp>=MOD) return 0; //p^(e1-e2-e3)是p^ai的倍数
a1=a1*temp%MOD;
}
int x,y; exgcd(a2*a3%MOD,MOD,x,y);
return a1*x%MOD;
}
int n,m,p,a[110],b[110],cnt;
signed main()
{
scanf("%lld%lld%lld",&n,&m,&p);
int k=sqrt(p),x=p;
for(int i=2;i<=k;++i)
if(x%i==0){
int t=1;
while(x%i==0)
x/=i,t*=i;
a[++cnt]=t;
b[cnt]=C(n,m,i,t);
}
if(x!=1){
a[++cnt]=x;
b[cnt]=C(n,m,x,x);
}
int A=a[1],B=b[1];
for(int i=2;i<=cnt;++i){
int k1,k2;
int g=exgcd(A,a[i],k1,k2);
int t=A;
A=A*a[i]/g;
B=(t*k1%A*(b[i]-B)/g%A+B)%A;
}
printf("%lld\n",(B+A)%A);
return 0;
}