思路
扩展Lucas和Lucas定理其实没什么关系
我们要求的是这样的一个问题
[left(egin{matrix}n\mend{matrix}
ight) mod P
]
p不一定是素数
所以需要CRT合并
问题转化为
[xequiv left(egin{matrix}n\mend{matrix}
ight) (mod p_1^{k_1}) \
xequiv left(egin{matrix}n\mend{matrix}
ight) (mod p_2^{k_2})\
dots\
xequiv left(egin{matrix}n\mend{matrix}
ight) (mod p_t^{k_t})
]
然后因为(p_1^{k_1},p_2^{k_2},dots,p_t^{k_t})互质，所以直接CRT
现在要求的是$ left(egin{matrix}nmend{matrix}
ight) (mod p_i^{k_i})(
由组合数的公式可知，要求的是)n! (modp_i^{k_i} )$
为了避免没有逆元，要先把阶乘中(p_i)全部消去，最后再乘回来（jc求质因数）
然后可以发现
有一部分是可以递归处理的（就是(n/pi)!）
有一部分在模pk意义下是有循环节的，枚举pk的长度，计算即可，出现了n/pk次
还有一部分剩下的，长度不会超过pk，暴力计算即可
然后就没了
代码
#include <cstdio>
#include <cstring>
#include <algorithm>
#define int long long
using namespace std;
int pow(int a,int b,int MOD){
int ans=1;
while(b){
if(b&1)
ans=(ans*a)%MOD;
a=(a*a)%MOD;
b>>=1;
}
return ans;
}
int exgcd(int a,int b,int &x,int &y){
if(b==0){
x=1;
y=0;
return a;
}
int req=exgcd(b,a%b,x,y);
int t=x;
x=y;
y=t-a/b*y;
return req;
}
int inv(int a,int p){
if(!a)
return 0;
int x,y;
exgcd(a,p,x,y);
x=((x%p+p)%p);
if(!x)
x+=p;
return x;
}
int mul(int n,int pi,int pk){//get n!/pi^a%p^k
if(!n)
return 1;
int ans=1;
for(int i=2;i<=pk;i++)
if(i%pi)
ans=(ans*i)%pk;
ans=pow(ans,n/pk,pk);
for(int i=2;i<=n%pk;i++)
if(i%pi)
ans=(ans*i)%pk;
return ans*mul(n/pi,pi,pk)%pk;
}
int C(int n,int m,int Mod,int pi,int pk){
if(m>n)
return 0;
int jcn=mul(n,pi,pk),jcm=mul(m,pi,pk),jcnm=mul(n-m,pi,pk),k=0;
for(int i=n;i;i/=pi)
k+=i/pi;
for(int i=m;i;i/=pi)
k-=i/pi;
for(int i=n-m;i;i/=pi)
k-=i/pi;
int ans=jcn*inv(jcm,pk)%pk*inv(jcnm,pk)%pk*pow(pi,k,pk)%pk;
return ans*(Mod/pk)%Mod*inv(Mod/pk,pk)%Mod;
}
int exLucas(int n,int m,int Mod){
int ans=0;
for(int i=2,t=Mod;i<=Mod;i++){
if(!(t%i)){
int midpk=1;
while(!(t%i)){
midpk*=i;
t/=i;
}
ans=(ans+C(n,m,Mod,i,midpk))%Mod;
}
}
return ans;
}
int n,m,MOD;
signed main(){
scanf("%lld %lld %lld",&n,&m,&MOD);
printf("%lld
",exLucas(n,m,MOD));
return 0;
}