拓展卢卡斯定理 / exlucas

恶心东西爬、、、

我们要求解一个 \(\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;
}