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;

}

巴特西

2023.2.26【模板】扩展Lucas定理

2023.2.26【模板】扩展Lucas定理

题目概述

算法流程

最新文章

热门文章