擴充盧卡斯定理 / exlucas

ChElsYqwq發表於2024-04-09

噁心東西爬、、、

我們要求解一個 \(\binom{n}{m}\mod M\)，\(M\) 是不太大的正整數，\(n,m\) 是可能比較大的正整數。

首先我們分解 \(M=\prod_{i=1}^kp_i^{x_i}\)，我們對於每一個 \(i\in[1,k]\) 求出 \(\binom{n}{m}\mod p_i^{x_i}\)，然後就會組成一個方程組，\(Ans\equiv\binom{n}{m}\pmod p_i^{x_i}\)，因為模數顯然互質，所以此時的 \(Ans\) 可以透過 CRT 求解。

問題是怎麼求這個 \(\binom{n}{m}\mod p^k\) 呢，我們知道有 \(a\) 在 \(\mod p\) 意義下有逆元的充要條件是 \(a\bot p\)，這個地方就顯得不一定有逆元了，不能直接求解，我們考慮變形：

\[\binom{n}{m}\equiv\frac{n!}{m!(n-m)!}\equiv\frac{\frac{n!}{p^x}}{\frac{m!}{p^y}\frac{(n-m)!}{p^z}}p^{x-y-z}\pmod {p^k} \]

然後求解 \(\frac{n!}{p^x}\pmod {p^k}\)，這個 \(x\) 儘量往大了取，先挑出 \(p\) 的倍數直接做。

\[\frac{n!}{p^x}\equiv p^{⌊\frac{n}{p}⌋}(⌊\frac{n}{p}⌋!)\prod_{i 不是 p 的倍數，i=1}^{n}i\pmod {p^k} \]

後面那一段其實是有一個長度為 \(p^k\) 的迴圈節的，手玩可得，好像是一個威爾遜定理的推論，所以我們可以拆開寫 >w<

\[\frac{n!}{p^x}\equiv p^{⌊\frac{n}{p}⌋}(⌊\frac{n}{p}⌋!)(\prod_{i 不是 p 的倍數，i=1}^{p^k} i)^{⌊\frac{n}{p^k}⌋}(\prod_{i 不是 p 的倍數,i=p^k⌊\frac{n}{p^k}⌋}^{n}i)\pmod {p^k} \]

寫成遞迴的函式形式，那麼 \(f(n)\equiv\frac{n!}{p^x}\pmod p\)，有下柿：

\[f(n)=f(⌊\frac{n}{p}⌋)(\prod_{i 不是 p 的倍數，i=1}^{p^k} i)^{⌊\frac{n}{p^k}⌋}(\prod_{i 不是 p 的倍數,i=p^k⌊\frac{n}{p^k}⌋}^{n}i)\pmod {p^k} \]

直接遞迴暴力做就行了謝謝喵，放回去原式，發現 \(\frac{\frac{n!}{p^x}}{\frac{m!}{p^y}\frac{(n-m)!}{p^z}}p^{x-y-z}\) 前面的分數部分可以做了，那麼還差一個 \(p^{x-y-z}\) 的指數，這個東西顯然是 \(f(n)\) 中省去的 \(p^{⌊\frac{n}{p}⌋}\) 的指數部分，可以直接累加這一部分。

然後終於他媽的做完了 /lh

#include<bits/stdc++.h>
#define int long long
#define up(i,l,r) for(int i=l; i<=r; ++i)
#define dn(i,r,l) for(int i=r; i>=l; --i)

using namespace std;

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

int inv(int a,int p) {
	int x, y; exgcd(a,p,x,y);
    return (x%p+p)%p;
}

int ksm(int a,int b,int p) {
	int res=1%p;
	for( ; b; b>>=1) {
		if(b&1) res=res*a%p;
		a=a*a%p;
	}
	return res;
}

int fac(int n,int pi,int pk) {
	if(!n) return 1;
	int res=1;
	up(i,2,pk) if(i%pi) res=res*i%pk;
	res=ksm(res,n/pk,pk);
	up(i,2,n%pk) if(i%pi) res=res*i%pk;
	return res*fac(n/pi,pi,pk)%pk;
}

int C(int n,int m,int pi,int pk) {
	int above=fac(n,pi,pk), l=fac(m,pi,pk), r=fac(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;
	return above*inv(l,pk)%pk*inv(r,pk)%pk*ksm(pi,k,pk)%pk;
}

int exlucas(int n,int m,int p) {
	int ans=0, x=p, t=sqrt(x);
	up(i,2,t) if(x%i==0) {
		int mul=1;
		while(x%i==0) x/=i, mul*=i;
		ans=(ans+C(n,m,i,mul)*(p/mul)%p*inv(p/mul,mul)%p)%p;
	}
	if(x>1) ans=(ans+C(n,m,x,x)*(p/x)%p*inv(p/x,x)%p)%p;
	return ans; 
}

signed main() {
	ios::sync_with_stdio(0);
	cin.tie(0);
	int n, m, p;
	cin >> n >> m >> p;
	cout << exlucas(n,m,p); 
	return 0;
}

相關文章