exLucas

参考博客

exLucas：

求 \(C_n^m\bmod d\) ( \(d\) 不一定为质数)

1.

将 \(d\) 质因数分解为 \(d=p_1^{c_1}\times p_2^{c_2}\times\cdots\times p_k^{c_k}\)

\(\forall i,j\in[1,k]\) , \(p_i^{c_i}\) 与 \(p_j^{c_j}\) 互质，所以可以构造出如下同余方程：

\[\begin{cases}a_1\equiv C_n^m\ (\bmod \ p_1^{c_1})\\a_2\equiv C_n^m\ (\bmod \ p_2^{c_2})\\\cdots\\a^k\equiv C_n^m\ (\bmod \ p_k^{c_k})\end{cases} \]

在求出所有的 \(a_i\) 后，就可以利用中国剩余定理求解出 \(C_n^m\) 。

所以问题转化为求出 \(C_n^m\bmod p^\alpha\) 的值。

2.

但是同余方程 \(ax\equiv 1(\bmod \ p)\) (乘法逆元) 有解的充要条件为 \(\gcd(a,p)=1\) 。

然而题目显然无法保证有解，所以无法直接求 \({\rm inv}_{m!}\) 和 \({\rm inv}_{(n-m)!}\)

所以可将原式转化为：

\[\frac{\frac{n!}{p^x}}{\frac{m!}{p^y}\times \frac{(n-m)!}{p^z}}\times p^{x-y-z}\bmod p^\alpha \]

其中 \(x\) 为 \(n!\) 中包含 \(p\) 的个数，\(y,z\) 同理。

3.

若对于一个 \(n\) 可以求出 \(\frac{n!}{p^x}\) ，那么就可以求逆元了，所以现在要求 \(\frac{n!}{p^x}\bmod p^k\) 。

对 \(n!\) 进行变形可得：

\[\begin{array}{l}n!=1\times 2\times \cdots\times n\\ =(p\times 2p\times \cdots)\times\prod_{i=1,i\not\equiv 0(\bmod p)}^n i\\=p^{\lfloor\frac{n}{p}\rfloor}\times (\lfloor\frac{n}{p}\rfloor)!\times\prod_{i=1,i\not\equiv 0(\bmod p)}^n i\\=p^{\lfloor\frac{n}{p}\rfloor}\times (\lfloor\frac{n}{p}\rfloor)!\times(\prod_{i=1,i\not\equiv 0(\bmod p)}^{p^k} i)^{\lfloor\frac{n}{p^k}\rfloor}\times(\prod_{i=p^k\lfloor\frac{n}{p^k}\rfloor,i\not\equiv 0(\bmod p)}^ni)\end{array} \]

设 \(f(n)=\frac{n!}{p^x}\ (f(0)=1)\)

\[\therefore f(n)=f(\lfloor\frac{n}{p}\rfloor)\times(\prod_{i=1,i\not\equiv 0(\bmod p)}^{p^k} i)^{\lfloor\frac{n}{p^k}\rfloor}\times(\prod_{i=p^k\lfloor\frac{n}{p^k}\rfloor,i\not\equiv 0(\bmod p)}^ni)\\ \]

\[\therefore 原式=\frac{f(n)}{f(m)f(n-m)}p^{x-y-z}\bmod p^k \]

这样就可以用 \(exgcd\) 求逆元了。

最后利用中国剩余定理得到答案。单次时间复杂度 \(O(p\log p)\)

代码

#include<bits/stdc++.h>
#define ll long long
using namespace std;

inline ll read(){
	ll s=0,k=1;
	char c=getchar();
	while(c>'9'||c<'0'){
		if(c=='-')k=-1;
		c=getchar();
	}
	while(c>='0'&&c<='9'){
		s=(s<<3)+(s<<1)+(c^48);
		c=getchar();
	}
	return s*k;
}

ll n,m,p; 

ll ksm(ll a,ll b,ll mod){
	ll t=1;
	for(;b;b>>=1,a=a*a%mod)
		if(b&1) t=t*a%mod;
	return t;
}

void exgcd(ll a,ll b,ll &x,ll &y){
	if(!b){
		x=1,y=0;
		return ;
	}
	exgcd(b,a%b,x,y);
	ll t=x;
	x=y;
	y=t-a/b*y;
}

ll inv(ll n,ll mod){
	ll x,y;
	exgcd(n,mod,x,y);
	return (x%mod+mod)%mod; 
}

ll CRT(int n,ll a[],ll P[]){
	ll ans=0;
	for(int i=1;i<=n;i++) ans=(ans+a[i]*(p/P[i])%p*inv(p/P[i],P[i])%p)%p;
	return ans;
}

ll f(ll n,ll pi,ll pk){
	if(!n) return 1;
	ll t=1;
	for(ll i=2;i<=pk;i++) if(i%pi) (t*=i)%=pk;
	t=ksm(t,n/pk,pk);
	for(ll i=n/pk*pk+1;i<=n;i++) if(i%pi) (t*=(i%pk))%=pk;
	return t*f(n/pi,pi,pk)%pk;
}

ll C(ll n,ll m,ll pi,ll pk){
	int k=0;
	for(ll i=n;i;i/=pi) k+=i/pi;
	for(ll i=m;i;i/=pi) k-=i/pi;
	for(ll i=n-m;i;i/=pi) k-=i/pi;
	return f(n,pi,pk)*inv(f(m,pi,pk),pk)%pk*inv(f(n-m,pi,pk),pk)%pk*ksm(pi,k,pk)%pk;
}

ll exLucas(ll n,ll m,ll p){
	int cnt=0;
	ll a[20],P[20],tmp=p;
	for(int i=2;i*i<=p;i++)
		if(tmp%i==0){
			P[++cnt]=1;
			while(tmp%i==0) P[cnt]*=i,tmp/=i;
			a[cnt]=C(n,m,i,P[cnt]);
		} 
	if(tmp>1) P[++cnt]=tmp,a[cnt]=C(n,m,tmp,tmp);
	return CRT(cnt,a,P);
}

int main(){
	n=read();m=read();p=read();
	printf("%lld\n",exLucas(n,m,p));
	return 0;
}