2023.2.26【模板】扩展Lucas定理

题目概述

求

由于

对于外面的\(p^{x - y - z}\)，只要求出\(x、y、z\)的值就可以计算了

观察以上函数可知，每次在\(f(n)\)这一层就会去掉\(\lfloor \frac np \rfloor\)个p因子

定义\(g(n)\)为\(n!\)中p因子的个数，则：

\[g(n) = g(\lfloor \frac np \rfloor) + \lfloor \frac np \rfloor\]

此结论对于其他题目也同样有效

所以原始式子就转化成了

\[\frac {f(n)}{f(m)f(n - m)} * p^{g(n) - g(m) - g(n - m)} \ mod \ p^{c_i}\]

因为去掉了p因子，所以\(f(m)\)和\(f(n - m)\)与\(p^{c_i}\)互质，可以求逆元

因为\(p^{c_i}\)不是质数，不能直接用费马小定理计算，所以我们采用\(exgcd\)求逆元

最后进行CRT合并答案

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
ll res[101],d[101],zs[101],tot = 0,M[101];
inline ll g(ll n,ll p)
{
	if(n == 0) return 0;
	return g(n / p,p) + n / p;
}
inline ll ksm(ll base,ll pts,ll mod)
{
	ll ret = 1;
	for(;pts > 0;pts >>= 1,base = base * base % mod)
		if(pts & 1)
			ret = ret * base % mod;
	return ret;
}
inline ll F(ll n,ll p,ll k)
{
	if(n == 0) return 1;
	ll P = ksm(p,k,1e18 + 1);
	ll mul = 1;
	for(ll i = 1;i <= P;i++)
		if(i % p)
			mul = mul * i % P;
	mul = ksm(mul,n / P,P);
	for(ll i = P * (n / P);i <= n;i++)
		if(i % p)
			mul = mul * (i % P) % P;
	return F(n / p,p,k) * mul % P;
}
inline void exgcd(ll a,ll b,ll &x,ll &y)
{
	if(b == 0)
	{
		x = 1;
    	y = 0;
   		return;
  	}
	ll tmp;
	exgcd(b,a % b,x,y);
	tmp = y;
	y = x - (a / b) * y;
	x = tmp;
}
inline ll exlucas(ll n,ll m,ll p)
{
	ll tmp = p;
	for(ll i = 2;i <= sqrt(p);i++)
	{
		if(tmp % i == 0)
		{
			++tot;
			d[tot] = i;
			while(tmp % i == 0)
			{
				tmp /= i;
				++zs[tot];
			}
		}
	}
	if(tmp != 1)
	{
		++tot;
		d[tot] = tmp;
		zs[tot] = 1; 
	}
	for(int i = 1;i <= tot;i++) 
	{
		ll P = ksm(d[i],zs[i],1e18 + 1);
 		ll inv1,inv2,yy;
 		exgcd(F(m,d[i],zs[i]),P,inv1,yy);
 		exgcd(F(n - m,d[i],zs[i]),P,inv2,yy);
		inv1 = (inv1 % P + P) % P;
  		inv2 = (inv2 % P + P) % P;
		res[i] = F(n,d[i],zs[i]) * inv1 % P * inv2 % P * ksm(d[i],g(n,d[i]) - g(m,d[i]) - g(n - m,d[i]),P) % P;
		M[i] = P;
	}
	ll ans = 0;
	for(int i = 1;i <= tot;i++)
	{
	    ll inv,yy;
	    exgcd(p / M[i],M[i],inv,yy);
	    inv = (inv % M[i] + M[i]) % M[i];
		ans = (ans + res[i] * (p / M[i]) % p * inv % p) % p;
	}
	return ans;
}
int main()
{
	ll n,m,p;
	cin>>n>>m>>p;
	cout<<exlucas(n,m,p);
	return 0;
}