2023.2.26【模板】扩展Lucas定理
题目概述
求\(\binom {n}{m} mod\) \(p\) 的值，不保证\(p\)为质数
算法流程
(扩展和普通算法毫无关系)
由于\(p\)不是质数，我们考虑[SDOI2010]古代猪文 - 洛谷中的处理方法：将\(p\)质因数分解得：
\[p = {p_1}^{c_1}{p_2}^{c_2}{p_3}^{c_3}....{p_k}^{c_k} \]所以我们考虑计算$\binom nm mod $ \({p_i}^{c_i}\)的值，再用CRT合并即可
展开上式：
\[\frac {n!}{m!(n - m)!} mod\ {p_i}^{c_i} \]我们发现由于\(m!(n - m)!\)中可能含有因数p，我们无法求出\(m!(n - m)!\)模\({p_i}^{c_i}\)意义下的逆元，所以我们考虑除去三个数中所有的p因子，假设\(p^x | n!\)且\(p^{x+1} \nmid n!\)，即x是\(n!\)中p因子的个数（对于\(m!\)和\((n - m)!\)同理）
\[\frac {\frac{n!}{p^x}}{\frac{m!}{p^y}\frac{(n - m)!}{p^z}}p^{x - y - z}\ mod\ {p_i}^{c_i} \]由于\(\frac{n!}{p^x}、\frac{m!}{p^y}、\frac{(n - m)!}{p^z}\)三式同构，我们考虑计算其中一个式子（以下用\(p\)替换\(p_i\)）
\[\frac {n!}{p^x}\ mod \ {p}^{c_i} \]展开为
\[\frac {1*2*3*...*n}{p^x}\ mod \ {p}^{c_i} \]提出p的倍数
\[\frac {(p * 2p * 3p * .. * \lfloor {\frac np} \rfloor p) * \Pi_{i = 1;i \not\equiv 0}^{n}}{p^x}\ mod\ {p}^{c_i} \]即
\[\frac {\lfloor {\frac np} \rfloor! * p^{\lfloor \frac np \rfloor} * \Pi_{i = 1;i \not\equiv 0}^{n}}{p^x}\ mod\ {p}^{c_i} \]如果暴力计算\(\Pi_{i = 1;i \not\equiv 0}^{n}\)复杂度过高，不难发现其有一个循环节，即每过p个数就会少乘上第p个数，又因为\({p_i}^{c_i}+ r \equiv r\ mod\ {p_i}^{c_i}\)，所以我们以\({p_i}^{c_i}\)作为这个循环节
\[\frac {\lfloor {\frac np} \rfloor! * p^{\lfloor \frac np \rfloor} * {[\Pi_{i = 1;i \not\equiv 0}^{p^{c_i}}]}^{\lfloor\frac {n}{p^{c_i}}\rfloor}\Pi_{i = {p^{c_i}}\lfloor\frac{n}{p^{c_i}}\rfloor;i \not\equiv 0}^{n}}{p^x}\ mod\ {p}^{c_i} \]对于\(\Pi_{i = 1;i \not\equiv 0}^{p^{c_i}}\)和\(\Pi_{i = {p^{c_i}}\lfloor\frac{n}{p^{c_i}}\rfloor;i \not\equiv 0}^{n}\)，暴力计算即可
不管\(x\)取何值，最终p因子都会消除，所以计算时去掉\(p^{\lfloor \frac np \rfloor}\)
因为\(\lfloor \frac np \rfloor!\)中可能含有p因子，所以我们将其进行递归：
设\(f(n) = \frac {n!}{p^x}\ mod \ {p}^{c_i}\)，则：
\[f(n) = {f(\lfloor {\frac np} \rfloor) * {[\Pi_{i = 1;i \not\equiv 0}^{p^{c_i}}]}^{\lfloor\frac {n}{p^{c_i}}\rfloor}\Pi_{i = {p^{c_i}}\lfloor\frac{n}{p^{c_i}}\rfloor;i \not\equiv 0}^{n}}\ mod\ {p}^{c_i} \]根据此式递推即可，时间复杂度为\(O(log_pn)\)，不会证明qwq
对于外面的\(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;
}