拓展卢卡斯定理 / exlucas

恶心东西爬、、、

我们要求解一个 $(\binom{n}{m}) \mod M$ ， $M$ 是不太大的正整数， $n, m$ 是可能比较大的正整数。

首先我们分解 $M = \prod_{i = 1}^{k} p_{i}^{x_{i}}$ ，我们对于每一个 $i \in [1, k]$ 求出 $(\binom{n}{m}) \mod p_{i}^{x_{i}}$ ，然后就会组成一个方程组， $A n s \equiv (\binom{n}{m}) (\mod p)_{i}^{x_{i}}$ ，因为模数显然互质，所以此时的 $A n s$ 可以通过 CRT 求解。

问题是怎么求这个 $(\binom{n}{m}) \mod p^{k}$ 呢，我们知道有 $a$ 在 $\mod p$ 意义下有逆元的充要条件是 $a ⊥ 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} (\mod p^{k})

然后求解 $\frac{n!}{p^{x}} (\mod p^{k})$ ，这个 $x$ 尽量往大了取，先挑出 $p$ 的倍数直接做。

\frac{n!}{p^{x}} \equiv p^{⌊ \frac{n}{p} ⌋} (⌊ \frac{n}{p} ⌋!) \prod_{i 不 是 p 的 倍 数 ， i = 1}^{n} i (\mod 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) (\mod p^{k})

写成递归的函数形式，那么 $f (n) \equiv \frac{n!}{p^{x}} (\mod 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) (\mod 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;
}