【luogu P4720】【模板】扩展卢卡斯定理/exLucas(数论)(CRT)

【模板】扩展卢卡斯定理/exLucas

题目链接:luogu P4720

题目大意

要你求一个组合数在模一个数的结果。
模数不一定是质数。

思路

我们考虑质因数分解模数。

考虑分成了 p=p1a1,p2a2,...

然后你考虑求出:
{Cnmmodp1a1Cnmmodp2a2

不难看出求出来的结果我们可以直接用 CRT(中国剩余定理)合并。

那接着就变成如何求 Cnmmodpk
首先我们用普通的式子:n!m!(nm)!modpk

显然是不能直接搞,因为我们不能求出下面的两个的逆元。
然后又逆元的条件是这个数和模数互质,所以我们可以换一下式子:
n!pxm!py(nm)!pzpxyzmodpk

然后 xn! 中包含 p 因子的个数。(y,z 同理)

那问题就是求 n!pxmodpk
考虑处理一下 n!
=(p2p3p...)(123...)
(左边是 1np 的倍数,右边就是不是倍数的)

然后提出左边的 p
=pnp(np)!i=1,i0(modp)ni

然后后面的显然是有循环节:
=pnp(np)!(i=1,i0(modp)pki)npk(i=pknpk,i0(modp)ni)

两个分别代表循环节和余项。

然后因为取模,pnp 会没掉,但是 (np)! 中可能还有 p
那就是一个递归下去的过程,设 F(n)=n!px

然后递推就是:
F(n)=F(np)(i=1,i0(modp)pki)npk(i=pknpk,i0(modp)ni)
(这个的复杂度是 logpn,边界是 f(0)=1

然后在原来的式子,就变成了 F(n)F(m)F(nm)pxyzmodpk
这下分母的两个就是一定跟 pk 互质,就可以用扩欧求啦。

然后你会发现还有一个 pxyz
那我们就是要求 F(n)=n!px 中的 x

想想上面的式子,我们要的就是被除掉的部分,就是 pnp 中的 np
每个 np 加起来就是结果,所以:
g(n)=np+g(np)
(复杂度也是 logpn,边界是 g(n)=0(n<p)

然后再带入到原来的式子就是:
F(n)F(m)F(nm)pG(n)G(m)G(nm)modpk

然后就好啦!!!

代码

#include<cstdio> #define ll long long using namespace std; ll n, m, p; ll a[1001], b[1001]; ll ksm(ll x, ll y, ll p) { ll re = 1; while (y) { if (y & 1) re = re * x % p; x = x * x % p; y >>= 1; } return re; } ll exgcd(ll a, ll b, ll &x, ll &y) { if (!b) { x = 1; y = 0; return a; } ll re = exgcd(b, a % b, y, x); y -= x * (a / b); return re; } ll inv(ll a, ll p) {//别用快速幂求,要用扩欧 ll x, y; exgcd(a, p, x, y); return (x % p + p) % p; } ll F(ll n, ll p, ll pk) { if (!n) return 1; ll rou = 1, el = 1; for (ll i = 1; i <= pk; i++) {//循环节部分 if (i % p) rou = rou * i % pk; } rou = ksm(rou, n / pk, pk); for (ll i = pk * (n / pk); i <= n; i++) {//余项部分 if (i % p) el = el * (i % pk) % pk; } return F(n / p, p, pk) * rou % pk * el % pk; } ll G(ll n, ll p) { if (n < p) return 0; return n / p + G(n / p, p); } ll work_C(ll n, ll m, ll p, ll pk) { ll fz = F(n, p, pk), fm1 = inv(F(m, p, pk), pk), fm2 = inv(F(n - m, p, pk), pk); ll mi = ksm(p, G(n, p) - G(m, p) - G(n - m, p), pk); return fz * fm1 % pk * fm2 % pk * mi % pk; } ll exLucas(ll n, ll m, ll p) { ll tmpp = p, tot = 0, di; for (ll i = 2; i * i <= tmpp; i++) { if (tmpp % i == 0) { tot++; di = 1; while (tmpp % i == 0) { di *= i; tmpp /= i; } a[tot] = di; b[tot] = work_C(n, m, i, di); } } if (tmpp > 1) a[++tot] = tmpp, b[tot] = work_C(n, m, tmpp, tmpp); ll ans = 0;//CRT for (int i = 1; i <= tot; i++) { ll M = p / a[i], V = inv(M, a[i]); ans = (ans + b[i] * M % p * V % p) % p; } return ans; } int main() { scanf("%lld %lld %lld", &n, &m, &p); printf("%lld", exLucas(n, m, p)); return 0; }

__EOF__

本文作者あおいSakura
本文链接https://www.cnblogs.com/Sakura-TJH/p/luogu_P4720.html
关于博主:评论和私信会在第一时间回复。或者直接私信我。
版权声明:本博客所有文章除特别声明外,均采用 BY-NC-SA 许可协议。转载请注明出处!
声援博主:如果您觉得文章对您有帮助,可以点击文章右下角推荐一下。您的鼓励是博主的最大动力!
posted @   あおいSakura  阅读(39)  评论(0编辑  收藏  举报
编辑推荐:
· go语言实现终端里的倒计时
· 如何编写易于单元测试的代码
· 10年+ .NET Coder 心语,封装的思维:从隐藏、稳定开始理解其本质意义
· .NET Core 中如何实现缓存的预热?
· 从 HTTP 原因短语缺失研究 HTTP/2 和 HTTP/3 的设计差异
阅读排行:
· 分享一个免费、快速、无限量使用的满血 DeepSeek R1 模型,支持深度思考和联网搜索!
· 基于 Docker 搭建 FRP 内网穿透开源项目(很简单哒)
· ollama系列01:轻松3步本地部署deepseek,普通电脑可用
· 25岁的心里话
· 按钮权限的设计及实现
点击右上角即可分享
微信分享提示