快速阶乘与（扩展）卢卡斯定理
(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} imes(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
",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
",(B+A)%A);
return 0;
}