扩展Lucas定理

我们首先来回顾一下Lucas定理：
(inom{n}{m} = inom{lfloorfrac{n}{p} floor}{lfloorfrac{m}{p} floor}inom{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{ egin{aligned} a1equivinom{n}{m}pmod{q_1^{alpha_1}}\ a2equivinom{n}{m}pmod{q_2^{alpha_2}}\ a3equivinom{n}{m}pmod{q_3^{alpha_3}}\ .....\ a4equivinom{n}{m}pmod{q_r^{alpha_r}} end{aligned} ight.)

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

Step2:

那么我们的问题转而求

(inom{n}{m}mod q^k(q is prime))。

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

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

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

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

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

对(n!)变形，则有：
(n! = p^{lfloor frac{n}{p} floor}(lfloorfrac{n}{p} floor)(prodlimits_{i = 1,!(iequiv 0pmod{p})}{i}))
显然后面这项拥有循环节。

则有:(n! = p^{lfloor frac{n}{p} floor}(lfloorfrac{n}{p} floor)(prodlimits_{i = 1,!(iequiv 0pmod{p})}i)^{lfloorfrac{n}{p^k} floor}(prodlimits_{i = p^k{lfloorfrac{n}{p^k} floor},!(iequiv 0pmod{p})}i))

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

所以(f(n) = f(lfloorfrac{n}{p} floor))(prodlimits_{i = 1,!(iequiv 0pmod{p})}i)^{lfloorfrac{n}{p^k} floor}(prodlimits_{i = p^k{lfloorfrac{n}{p^k} floor},!(iequiv 0pmod{p})}i))

设 (g(n))为(n!)的(p)的因子个数：
则有 (g(n) = lfloorfrac{n}{p} floor + g(lfloorfrac{n}{p} floor))

所以答案为 (frac{f(n)}{f(m)f(n -m)}p^{g(n) - g(m) - g(n - m)}mod 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;
}