[学习笔记]扩展Lucas定理

我们首先来回顾一下Lucas定理：
\(\binom{n}{m} = \binom{\lfloor\frac{n}{p}\rfloor}{\lfloor\frac{m}{p}\rfloor}\binom{n\ mod\ p}{m\ mod\ p}(mod\ p)(p\ is\ prime)\)

其给出了在膜数为质数时，组合数取膜的优秀性质。
其使用的地方为当 \(n,m \geq p\)时。

ll C(ll x,ll y){
	if(y > x)return 0;
	return s[x] * inv[y] % p * inv[x - y] % p;
}

inline ll Lucas(ll x,ll y){
	if(y == 0)
	return 1;
	return Lucas(x / p,y / p) * C(x % p,y % p) % p;
}

inline void solve(){
	s[0] = 1;
	scanf("%lld%lld%lld",&n,&m,&p);
	for(int i = 1;i <= p - 1;++i)
	s[i] = s[i - 1] * i % p;
	inv[p - 1] = pow(s[p - 1],p - 2);
	for(int i = p - 2;i >= 0;--i)
	inv[i] = inv[i + 1] * (i + 1) % p;
	std::cout<<Lucas(n + m,n)<<std::endl;
}

那么我们自然思考当膜数并非质数的时候，我们有扩展Lucas来应对他。

我们首先把 \(p\) 质数分解：
\(p = q_1^{\alpha_1}q_2^{\alpha_2}q_3^{\alpha_3}....q_r^{\alpha_r}\)

Step1:

那么我们就有若干个同余方程
\(\left\{ \begin{aligned} a1\equiv\binom{n}{m}\pmod{q_1^{\alpha_1}}\\ a2\equiv\binom{n}{m}\pmod{q_2^{\alpha_2}}\\ a3\equiv\binom{n}{m}\pmod{q_3^{\alpha_3}}\\ .....\\ a4\equiv\binom{n}{m}\pmod{q_r^{\alpha_r}} \end{aligned} \right.\)

那么依据中国剩余定理我们就能求出 \(\binom{n}{m}\)。

Step2:

那么我们的问题转而求

\(\binom{n}{m}\bmod q^k(q\ is\ prime)\)。

那么根据定义：
我们有 \(\frac{n!}{m!(n - m)!}\bmod q^k\)

但是我们发现不一定有逆元。
所以我们将原式转换为:

\(\frac{\frac{n!}{q^x}}{\frac{m!}{q^y}\frac{(n - m)!}\ {q^z}}q^{x-y-z}\bmod q^k\)

其中\(x,y,z\)表示在\(n!\)里的\(p\)的因子个数，其他类似。

那么我们转为求 \(\frac{n!}{q^x}\bmod q^k\)

对\(n!\)变形，则有：
\(n! = p^{\lfloor \frac{n}{p} \rfloor}(\lfloor\frac{n}{p} \rfloor)(\prod\limits_{i = 1,!(i\equiv 0\pmod{p})}{i})\)
显然后面这项拥有循环节。

则有:\(n! = p^{\lfloor \frac{n}{p} \rfloor}(\lfloor\frac{n}{p} \rfloor)(\prod\limits_{i = 1,!(i\equiv 0\pmod{p})}i)^{\lfloor\frac{n}{p^k}\rfloor}(\prod\limits_{i = p^k{\lfloor\frac{n}{p^k}\rfloor},!(i\equiv 0\pmod{p})}i)\)

定义：
\(f(n) = \frac{n!}{p^x}\)

所以\(f(n) = f(\lfloor\frac{n}{p}\rfloor))(\prod\limits_{i = 1,!(i\equiv 0\pmod{p})}i)^{\lfloor\frac{n}{p^k}\rfloor}(\prod\limits_{i = p^k{\lfloor\frac{n}{p^k}\rfloor},!(i\equiv 0\pmod{p})}i)\)

设 \(g(n)\)为\(n!\)的\(p\)的因子个数：
则有 \(g(n) = \lfloor\frac{n}{p}\rfloor + g(\lfloor\frac{n}{p}\rfloor)\)

所以答案为 \(\frac{f(n)}{f(m)f(n -m)}p^{g(n) - g(m) - g(n - m)}\bmod p^k\)
逆元扩展\(exgcd\)就行了。

#include<iostream>
#include<cstdio>
#define ll long long 

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

ll gcd(ll a,ll b){
	if(b == 0)return a;
	return gcd(b,a % b);
}

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

inline ll lcm(ll a,ll b){
	return a / gcd(a,b) * b;
}

inline ll fastpow(ll a,ll b,ll p){
	ll ans = 1;
	a %= p;
	while(b){
		if(b & 1)ans = (ans * a) % p;
		b >>= 1;
		a = (a * a) % p;
	}
	return ans;
}

inline ll read(){
	ll ans = 0;
	char a = getchar();
	while(!(a <= '9' && a >= '0'))a = getchar();
	while(a <= '9' && a >= '0')ans = (ans << 3) + (ans << 1) + (a - '0'),a = getchar();
	return ans;
}

inline ll f(ll n,ll p,ll pk){
	if(n == 0)return 1;
	ll rou = 1;
	ll res = 1;
	for(ll i = 1;i <= pk;++i)
	if(i % p)rou = rou * i % pk;
	rou = fastpow(rou,n / pk,pk);
	for(ll i = pk * (n / pk);i <= n;++i)
	if(i % p)res = res * (i % pk) % pk;
	return f(n / p,p,pk) * rou % pk * res % pk;
}

inline ll g(ll n,ll p){
	if(n < p)return 0;
	return g(n / p,p) + (n / p);
}

inline ll c_pk(ll n,ll m,ll p,ll pk){
	ll fn = f(n,p,pk),fm = inv(f(m,p,pk),pk),fnm = inv(f(n - m,p,pk),pk);
	ll mi = fastpow(p,g(n,p) - g(m,p) - g(n - m,p),pk);
	return fn % pk * fm % pk * fnm % pk * mi % pk;
}

ll A[1001],B[1001];

// x = B(mod A)

inline ll exlucas(ll n,ll m,ll p){
	ll P = p,tot = 0;
	for(ll i = 2;i * i <= P;++i){
		if(!(p % i)){
			ll pk = 1;
			while(!(p % i))
			pk *= i,p /= i;
			A[++tot] = pk;
			B[tot] = c_pk(n,m,i,pk);
		}
	}
	if(p != 1)
	A[++tot] = p,B[tot] = c_pk(n,m,p,p);
	ll ans = 0;
	for(ll i = 1;i <= tot;++i){
		ll M = P / A[i],t = inv(M,A[i]);
		ans = (ans + B[i] * M % P * t % P) % P;
	}
	return ans;
}

int main(){
	ll n = read(),m = read(),p = read();
	std::cout<<exlucas(n,m,p)<<std::endl;
}