LG4720 【模板】扩展卢卡斯定理

扩展卢卡斯定理

求 (C_n^m mod{p})，其中 (C) 为组合数。

(1≤m≤n≤10^{18},2≤p≤1000000) ，不保证 (p) 是质数。

Fading的题解

设

[p=p_1^{alpha_1}p_2^{alpha_2}cdots p_k^{alpha_k} ]

求出

[left{egin{align*} C_n^m & mod & {p_1^{alpha_1}} \ C_n^m & mod & {p_2^{alpha_2}}\ vdots\ C_n^m & mod & {p_k^{alpha_k}} end{align*} ight. ]

最后用中国剩余定理合并即可。

假设我们现在要求

[C_n^mmod{P^k}\ =frac {n!}{m!(n-m)!}mod{P^k} ]

我们无法求得(m!,(n-m)!)的逆元，理由很简单，不一定有解。

怎么办呢？发现(a)对(p)有逆元的充要条件为(gcd(a,p)=1)。所以我们换一个式子：

[frac {frac {n!}{P^{x}}}{frac {m!}{P^{y}}frac {(n-m)!}{P^{z}}}P^{x-y-z}mod{P^k} ]

其中(x)为(n!)中(P)因子的个数（包含多少(P)因子）。其他都是同理的。

所以如果我们对于一个(n)，可以求出(frac {n!}{P^x})，那不就完事了吗？（这样我们就可以求逆元了）

问题等价于求(frac {n!}{P^x}mod{P^k})。我们对(n!)进行变形：

[n!=1cdot 2cdot 3cdots n=(Pcdot 2Pcdot 3Pcdots)(1cdot 2cdots) ]

左边是(1sim n)中所有(leq n)的(P)的倍数，右边是其它剩余的。

因为(1sim n)中有(lfloor frac nP floor)个(P)的倍数，

[ herefore P^{lfloor frac nP floor}(1 cdot 2 cdot 3cdots)(1cdot 2cdots) =P^{lfloor frac nP floor}(lfloor frac nP floor)!prod_{i=1,i otequiv 0pmod {P}}^ni ]

显然后面这个(mod P)是有循环节的。

[=P^{lfloor frac nP floor}(lfloor frac nP floor)!(prod_{i=1,i otequiv 0pmod {P}}^{P^k}i)^{lfloor frac n{P^k} floor}(prod_{i=P^klfloor frac n{P^k} floor,i otequiv 0pmod {P}}^ni) ]

其中((prod_{i=1,i otequiv 0pmod {P}}^{P^k}i)^{lfloor frac n{P^k} floor})是循环节(1sim P)中所有无(P)因子数的乘积，((prod_{i=P^klfloor frac n{P^k} floor,i otequiv 0pmod {P}}^ni))是余项的乘积。

比如

[22!equiv(1cdot 2cdot 4cdot 5cdot 7cdot 8)(10cdot 11cdot 13cdot 14cdot 16cdot 17)(19cdot 20cdot 22)(3cdot 6cdot 9cdot 12cdot 15cdot 18cdot 21)pmod {3^2}\ equiv(1cdot 2cdot 4cdot 5cdot 7cdot 8)^2(19cdot 20cdot 22)3^7(1cdot 2cdot 3cdot 4cdot 5cdot 6cdot 7)pmod {3^2}\ equiv3^77!(1cdot 2cdot 4cdot 5cdot 7cdot 8)^2(19cdot 20cdot 22)pmod {3^2} ]

正好是一样的，$lfloor frac {22}{3^2} floor=2 $。理解了吧...

发现前面这个(P^{lfloor frac nP floor})是要除掉的，然而((lfloor frac nP floor)!)里可能还包含(P)。

所以，我们定义函数：

[f(n)=frac {n!}{P^x}\ =f(lfloor frac nP floor)(prod_{i=1,i otequiv 0pmod {P}}^{P^k}i)^{lfloor frac n{P^k} floor}(prod_{i=P^klfloor frac n{P^k} floor,i otequiv 0pmod {P}}^ni)) ]

这样就好了。时间复杂度是(O(log_{P}n))。边界(f(0)=1)。

看看原式

[frac {frac {n!}{P^{x}}}{frac {m!}{P^{y}}frac {(n-m)!}{P^{z}}}P^{x-y-z}mod{P^k}\ =frac {f(n)}{f(m)f(n-m)}P^{x-y-z}mod{P^k} ]

由于(f(m))一定与(P^k)互质，所以我们可以直接用(exgcd)求逆元啦。

最后我们还有一个问题，(P^{x-y-z})咋求?其实很好求。

比如说，我们要求(f(n)=frac {n!}{P^x})中的(x)。

设(g(n)=x)(上述的(x))，再看看阶乘的式子

[n!=P^{lfloor frac nP floor}(lfloor frac nP floor)!(prod_{i=1,i otequiv 0pmod {P}}^{P^k}i)^{lfloor frac n{P^k} floor}(prod_{i=P^klfloor frac n{P^k} floor,i otequiv 0pmod {P}}^ni) ]

很显然我们要的就是(P^{lfloor frac nP floor})，而且((lfloor frac nP floor)!)可能还有(P)因子，所以我们可以得到以下递推式

[g(n)=lfloor frac nP floor+g(lfloor frac nP floor) ]

时间复杂度是(O(log_{P}n))。边界(g(n)=0(n<P))

所以答案就是

[frac {frac {n!}{P^{x}}}{frac {m!}{P^{y}}frac {(n-m)!}{P^{z}}}P^{g(x)-g(y)-g(z)}mod{P^k} ]

最后用中国剩余定理合并答案即可得到(C_n^m)。总时间复杂度(O(P log P))。

#include<bits/stdc++.h>
#define co const
#define il inline
template<class T> T read(){
	T x=0,w=1;char c=getchar();
	for(;!isdigit(c);c=getchar())if(c=='-') w=-w;
	for(;isdigit(c);c=getchar()) x=x*10+c-'0';
	return x*w;
}
template<class T> il T read(T&x){
	return x=read<T>();
}
using namespace std;
#define int long long

#define gcd __gcd
int pow(int a,int b,int mod){
	a=(a%mod+mod)%mod;
	int ans=1;
	for(;b;b>>=1,a=a*a%mod)
		if(b&1) ans=ans*a%mod;
	return ans;
}
void exgcd(int a,int b,int&x,int&y){
	if(!b) x=1,y=0;
	else exgcd(b,a%b,y,x),y-=a/b*x;
}
int inv(int a,int mod){
	if(!a) return 0;
	if(gcd(a,mod)!=1) return -1;
	int x,y;
	exgcd(a,mod,x,y);
	return (x%mod+mod)%mod;
}
int fac(int n,int pi,int pk){ //n! % pi^ki pk=pi^ki
	if(n<=1) 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*fac(n/pi,pi,pk)%pk;
}
int binom(int n,int m,int pi,int pk){
	int up=fac(n,pi,pk),dl=fac(m,pi,pk),dr=fac(n-m,pi,pk);
	int ex=0;
	for(int i=n;i;i/=pi) ex+=i/pi;
	for(int i=m;i;i/=pi) ex-=i/pi;
	for(int i=n-m;i;i/=pi) ex-=i/pi;
	return up*inv(dl,pk)%pk*inv(dr,pk)%pk*pow(pi,ex,pk)%pk;
}
signed main(){
	int n=read<int>(),m=read<int>(),P=read<int>();
	int ans=0,p=P;
	for(int i=2;i<=p;++i)if(p%i==0){
		int pi=i,pk=1;
		while(p%i==0) p/=i,pk*=i;
		ans=(ans+binom(n,m,pi,pk)*(P/pk)%P*inv(P/pk,pk)%P)%P;
	}
	printf("%lld
",ans);
	return 0;
}