【学习笔记】exLucas

exLucas

前置知识：EXCRT（或 CRT）

题面

给定整数 \(n,m,p\) 的值，求出 \(\binom{n}{m} \bmod p\) 的值，但不保证 \(p\) 为质数，\(n,m \leq 10^{18}\)。

做法

此时用 Lucas 定理显然是不行的（说实话 exLucas 的做法跟 Lucas 定理几乎没什么关系）。

因此我们需要将其进行一些转化。

我们令 \(p=\prod\limits_{i=1}^{k} p_i^{\alpha_i}(p_i\in prime)\)。则有 \(\binom{n}{m} \bmod p\) 等价于

\[\begin{cases} \binom{n}{m} \equiv x_1 \pmod{p_1^{\alpha_1}}\\ \binom{n}{m} \equiv x_2 \pmod{p_2^{\alpha_2}}\\ \ \ \ \ \ \ \dots \\ \binom{n}{m} \equiv x_k \pmod{p_k^{\alpha_k}}\\ \end{cases} \]

因此我们只需求出所有的 \(x_i\) 即可。

对于单独的一个式子 \(\binom{n}{m} \equiv x \pmod{p^{\alpha}}\)，需要将其求解。

若直接求出 \(\binom{n}{m} \bmod p^\alpha\) 显然是会 TLE 的，但直接求 \(x! \bmod p^\alpha\) 与它的逆元会发现 \(x!\) 与 \(p^\alpha\) 不一定互质。我们需将其转化一下。

设 \(f_i\) 为 \(i!\) 中含 \(p\) 因子的个数，则有 \(f_i = \lfloor { \frac{i}{p} } \rfloor + f_{ \lfloor { \frac{i}{p} } \rfloor }\)。

设 \(k_0\) 为 \(p^k \mid \binom{n}{m}\) 中 \(k\) 的最大值，有 \(k_0=f_n-f_m-f_{n-m}\)。

则 \(\binom{n}{m}\) 为 \(p^{k_0}\) 乘上一个数。

我们令 \(fact_i=\frac{i!}{p^{f_i}}\)，可以发现 \(fact_i=1 \times 2 \times \dots \times (p-1) \times (p+1) \times \dots \times (2p-1) \times(2p+1) \times \dots\)

因为 \(p\) 为质数且 \(fact_i\) 不含 \(p\) 的因子，所以 \(fact_i \perp p^\alpha\)

再设

\[w_t=\prod\limits_{i=1}^t \begin{cases} 1 \ ,\ \ p \mid i\\ i \ \ ,\ \ p \nmid i \end{cases} \]

则有 \(fact_i=w_i \times fact_{\lfloor { \frac{i}{p} } \rfloor}\)。单个求解可以通过递归实现。

因为 \(n\) 很大，用数组难以存下，因此我们可以先求出 \(w_{0 \sim p^k-1}\)，则有 \(w_i=(w_{p^k-1})^{\lfloor { \frac{i}{p^k} } \rfloor} \times w_{n \ \bmod \ p^k}\)。因此可以用 exLucas 求解的条件是 p 不大才行。

所以 \(\binom{n}{m} \bmod p^k = p^{k_0}\times fact_n \times fact_m^{-1} \times fact_{n-m}^{-1}\)。

最后通过 EXCRT（或 CRT）求解即可。

Code

#include<bits/stdc++.h>
#define IOS cin.tie(0),cout.tie(0),ios::sync_with_stdio(0)
#define mod 998244353
#define ll long long
#define db double
#define pb push_back
#define MS(x,y) memset(x,y,sizeof x)
using namespace std;
const int N=1e6+5,M=1e5+5;
const ll INF=1ll<<60;
ll n,m,cnt,tP,p[N],k[N],w[N];
ll exgcd(ll a,ll b,ll &x,ll &y){
	if(!b){x=1,y=0;return a;}
	else{ll tp=exgcd(b,a%b,y,x);y-=a/b*x;return tp;}
}
namespace exCRT{
	ll A=1,B=0;
	void add(ll a,ll b){
		ll x=0,y=0;
		ll gd=exgcd(A,a,x,y);
		x=(B-b)/gd*x;
		B=B-A*x;
		A=a/gd*A;
		B=(B%A+A)%A;
	}
}
ll qpow(ll a,ll b,ll P){//快速幂 
	ll ans=1;
	for(;b;b>>=1,a=a*a%P) if(b&1) ans=ans*a%P;
	return ans;
}
ll getf(ll xx,ll P){//求解 f 
	if(xx<P) return 0;
	return (xx/P)+getf(xx/P,P);
}
ll getw(ll xx,ll pk){//求解 w 
	return qpow(w[pk-1],xx/pk,pk)*w[xx%pk]%pk;
}
ll fact(ll xx,ll pk,ll P){//求解 fact 
	if(!xx || xx==1) return 1;
	return getw(xx,pk)*fact(xx/P,pk,P)%pk;
}
ll inv(ll xx,ll pk){//逆元 
	ll x=0,y=0;
	exgcd(xx,pk,x,y);
	return (x%pk+pk)%pk;
}
void init(){cnt=0;MS(p,0);MS(k,0);}
ll exLucas(ll n,ll m,ll tP){
	init();
	for(int i=2;i*i<=tP;i++){//分解质因数 
		if(!(tP%i)){
			p[++cnt]=i;
			while(!(tP%i)) k[cnt]++,tP/=i;
		}
	}
	if(tP>1) p[++cnt]=tP,k[cnt]=1;
	ll fn,fm,fnm,kmax,pk,b;
	for(int i=1;i<=cnt;i++){
		MS(w,0);
		fn=getf(n,p[i]),fm=getf(m,p[i]),fnm=getf(n-m,p[i]);
		kmax=fn-fm-fnm;//求k0 
		pk=qpow(p[i],k[i],LONG_LONG_MAX);//p^k 
		w[0]=1;
		for(int j=1;j<pk;j++) w[j]=w[j-1]*((ll)j%p[i]?j:1)%pk;//求 w 
		ll ifm=inv(fact(m,pk,p[i]),pk),ifnm=inv(fact(n-m,pk,p[i]),pk);
		ll pkmax=qpow(p[i],kmax,pk);
		b=pkmax*fact(n,pk,p[i])%pk*ifm%pk*ifnm%pk;//(n,m) mod p^k
		exCRT::add(pk,b);
	}
	return exCRT::B;
}
int main(){
	IOS;cin>>n>>m>>tP;
	cout<<exLucas(n,m,tP)<<'\n';
	return 0;
}

\[\large{\mathbb{END}} \]