组合数取模 - Lucas/exLucas - LibreOJ #181. 二项式系数

文章目录

Lucas定理
证明
拓展问题
一些优化（LibreOJ #181. 二项式系数）
- 预处理部分
- 复杂度分析
- 大佬的玄学优化（取模）
- AC代码

证明和代码分开，可以根据自己的需要跳转。

去我的Blog看，可能阅读效果更好

Lucas定理

(nm)modp=(⌊np⌋⌊mp⌋)(n%pm%p)modp\begin{aligned} {n \choose m} \mod p = {\lfloor \frac{n}{p} \rfloor \choose \lfloor \frac{m}{p} \rfloor} {{n \% p} \choose {m \% p}} \mod p \end{aligned} (mn)modp=(⌊pm⌋⌊pn⌋)(m%pn%p)modp
上式当ppp为素数时成立。

此时，若ppp为一个较小的素数，而n,mn,mn,m为一个较大的数（指不能直接通过预处理阶乘及其逆元的大数，比如1e181e181e18），那么我们就可以递归地求解该结果。

inline ll func(ll n, ll m, ll p, ll res = 1LL) { // n<=mif (n < m) return 0;while (m) res = res * C(n % p, m % p) % p, n /= p, m /= p;return res;
}

根据上述代码，我们只要预处理出ppp以内的阶乘及其逆元即可，时间复杂度O(p+log⁡m)O(p + \log m)O(p+logm)。

证明

其实在上面的递归过程已经提示了一件事：ppp进制分解。

考虑n,mn,mn,m在ppp进制的表达式n=∑i=0kakpk,m=∑i=0kbkpkn={\sum_{i=0}^{k}{a_kp^k}},m={\sum_{i=0}^{k}{b_kp^k}}n=∑i=0kakpk,m=∑i=0kbkpk，

那么Lucas定理实际就是：
(nm)≡(a0b0)(a1b1)…(akbk)(modp){n \choose m}\equiv{a_0 \choose b_0}{a_1 \choose b_1}\dots{a_k \choose b_k}\pmod{p}(mn)≡(b0a0)(b1a1)…(bkak)(modp)

现在先来看一个简单的式子：

(a+b)p=∑i=0p(pi)aibp−i≡ap+bp(modp)(a+b)^p={\sum_{i=0}^{p}{p \choose i}{a^ib^{p-i}}}\equiv{a^p+b^p}\pmod{p}(a+b)p=∑i=0p(ip)aibp−i≡ap+bp(modp)
因为其余的部分都有因子ppp

那么我们可以得到：
(1+x)n≡(1+x)∑aipi≡∏(1+xpi)ai(modp){(1+x)^n}\equiv{(1+x)^{\sum{a_ip^i}}}\equiv{\prod(1+x^{p^i})^{a_i}}\pmod{p}(1+x)n≡(1+x)∑aipi≡∏(1+xpi)ai(modp)
现在考虑xmx^mxm的系数，左边是(nm){n \choose m}(mn)，右边是∏(aibi)\prod{a_i \choose b_i}∏(biai)，已经是Lucas的表达式啦。

拓展问题

如果上式中ppp为合数，同时n,mn,mn,m依旧很大，怎么求解呢？

由中国剩余定理可以知道，

如果合数ppp的唯一分解为∏i=1spiαi\prod_{i=1}^{s}{p_i}^{\alpha_i}∏i=1spiαi，则有：
x≡(nm)(modp)⇒{x≡(nm)(modp1α1)x≡(nm)(modp2α2)⋮x≡(nm)(modpsαs)x \equiv {n \choose m} \pmod{p} \Rightarrow \begin{cases} x & \equiv & {n \choose m} \pmod{{p_1}^{\alpha_1}} \\ x & \equiv & {n \choose m} \pmod{{p_2}^{\alpha_2}} \\ & \vdots & \\ x & \equiv & {n \choose m} \pmod{{p_s}^{\alpha_s}} \\ \end{cases} x≡(mn)(modp)⇒⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧xxx≡≡⋮≡(mn)(modp1α1)(mn)(modp2α2)(mn)(modpsαs)

所以只要求出ppp的唯一分解，再分别解出各个同余式，通过CRT合并答案就能完成全部解答。

现在的问题就转变为，如何求出(nm)modpiαi{n \choose m} \mod{{p_i}^{\alpha_i}}(mn)modpiαi。

我们知道(nm)=n!m!(n−m)!{n \choose m}=\frac{n!}{m!(n-m)!}(mn)=m!(n−m)!n!，如果m!m!m!与pαp^{\alpha}pα互素，那么只要拓展欧几里得求出阶乘逆元就能直接算出结果。但是，先不论n,mn,mn,m的大小是否支持我们直接求出阶乘。它们是否互素呢？显然是否定的。

不过，如果能分离出不互素的部分和互素的部分，我们可以就求逆元得到答案了。

注意到，pαp^{\alpha}pα只有一种质因子ppp，那么阶乘中与之不互素的部分显然都是与ppp不互素的。这种不互素的数的个数也可以轻松算出，即potp(n!)=∑i=1∞⌊npi⌋{pot_p(n!)}={\sum_{i=1}^{\infty}{\lfloor \frac{n}{p^i} \rfloor}}potp(n!)=∑i=1∞⌊pin⌋。然后在模pαp^{\alpha}pα意义下所有数在Zpα\Z_{p^{\alpha}}Zpα中，（注意pαp^{\alpha}pα比较小）也就是说大于pαp^{\alpha}pα的数通过取模必然在[0,p−1][0,p-1][0,p−1]中，且很容易知道这样的数是循环出现的。

所以现在的思路就是，预处理pαp^{\alpha}pα内所有与ppp互素的数的乘积（令为basebasebase），那么n!n!n!中必然有⌊npα⌋\lfloor \frac{n}{p^{\alpha}} \rfloor⌊pαn⌋个这样的乘积，剩余部分为n%pαn\%p^{\alpha}n%pα内与ppp互素的数的乘积。

fac(n)=base⌊npα⌋fac(n%pα)fac(n)=base^{\lfloor \frac{n}{p^{\alpha}} \rfloor}fac(n\%p^{\alpha})fac(n)=base⌊pαn⌋fac(n%pα)

总结一下就是，现在我们把阶乘在O(pα+log⁡pn)O({p^{\alpha}}+\log_{p}{n})O(pα+logpn)时间上，分成了互素数乘积fac(n)和不互素的p^pot两部分。

所以，上代码了。

inline ll pot(ll n, ll p, ll res = 0) {while (n) res += n / p, n /= p;return res;
}inline ll fac(ll n, ll p, ll mod, ll res = 1LL) {if (n == 0) return 1LL;ll base = 1LL;for (ll i = 2; i < mod; i++) if (i % p) base = mul(base, i, mod);while (n) {ll tmp = qpow(base, n / mod, mod), sz = n % mod;for (ll i = 2; i <= sz; i++) if (i % p) tmp = mul(tmp, i, mod); res = mul(res, tmp, mod), n /= p;} return res;
}

那么，现在我们就能得到：
(nm)≡fac(n)×inv(fac(m))×inv(fac(n−m))×ppotp(n)−potp(m)−potp(n−m)(modpα){n \choose m} \equiv fac(n) \times inv(fac(m)) \times inv(fac(n-m)) \times p^{pot_p(n)-pot_p(m)-pot_p(n-m)} \pmod{p^{\alpha}} (mn)≡fac(n)×inv(fac(m))×inv(fac(n−m))×ppotp(n)−potp(m)−potp(n−m)(modpα)
之后就是简单的CRT啦，这个就不单独上代码了。

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

直接给出全部代码吧：

因为最近沉迷压行，可能看起来有点不舒服，但是我舒服就行但是可以复制到本地去格式化，所以我就不改啦。

// #pragma GCC optimize("Ofast")
// #pragma GCC target("avx,avx2,fma")#include <bits/stdc++.h>
using namespace std;
#define fast ios_base::sync_with_stdio(false), cin.tie(0), cout.tie(0)
#define fin freopen("input.txt", "r", stdin), freopen("output.txt", "w", stdout)
#define PII pair<int, int>
#define PLL pair<ll, ll>
#define pb push_back
#define mp make_pair
#define fi first
#define se second
#define ld nd << 1
#define rd nd << 1 | 1typedef long double lld;
typedef long long ll;
typedef unsigned long long ull;
// typedef __int128_t i128;
typedef vector<ll> vec;// const int mod = 1e9 + 7;
// const int mod = 998244353;
const int maxn = 1e6 + 50;
const int inf = 0x3f3f3f3f;            // 1e9;
const ll INF = 0x3f3f3f3f3f3f3f3f; // 1e18;int T;/* inline ll mul(ll a, ll b, ll mod, ll res = 0) {for ( ; b; b >>= 1, a = (a + a) % mod) if (b & 1) res = (res + a) % mod;return res;
} */inline ll mul(ll a, ll b, ll p) { return (a * b - (ll)((lld)a / p * b) * p + p) % p; } inline ll qpow(ll x, ll m, ll mod, ll res = 1) {for ( ; m; m >>= 1, x = mul(x, x, mod)) if (m & 1) res = mul(res, x, mod);return res;
}inline void exgcd(ll a, ll b, ll &g, ll &x, ll &y) {if (b) exgcd(b, a % b, g, y, x), y -= a / b * x;else x = 1, y = 0, g = a;
}inline ll inv(ll x, ll m) {ll a, b, g; exgcd(x, m, g, a, b);return (a % m + m) % m;
}ll r[maxn], m[maxn], tot = 1;
inline ll CRT() {ll M = 1LL, ans = 0LL;for (int i = 1; i < tot; i++) M *= m[i];for (int i = 1; i < tot; i++) { // combinell Mi = M / m[i], Mi_inv = inv(Mi, m[i]);ans = (ans + mul(mul(r[i], Mi, M), Mi_inv, M)) % M;} return ans;
}inline ll pot(ll n, ll p, ll res = 0) {while (n) res += n / p, n /= p;return res;
}inline ll fac(ll n, ll p, ll mod, ll res = 1LL) {if (n == 0) return 1LL;ll base = 1LL;for (ll i = 2; i < mod; i++) if (i % p) base = mul(base, i, mod);while (n) {ll tmp = qpow(base, n / mod, mod), sz = n % mod;for (ll i = 2; i <= sz; i++) if (i % p) tmp = mul(tmp, i, mod); res = mul(res, tmp, mod), n /= p;} return res;
}inline ll exLucas(ll a, ll b, ll p) {for (ll i = 2; i * i <= p; i ++) {if (p % i == 0) {ll t = 1LL;while (p % i == 0) p /= i, t *= i;ll fa = fac(a, i, t), fb = fac(b, i, t), fc = fac(a - b, i, t);ll pk = qpow(i, pot(a, i) - pot(b, i) - pot(a - b, i), t);m[tot] = t, r[tot ++] = mul(pk, mul(fa, mul(inv(fb, t), inv(fc, t), t), t), t);}} if (p ^ 1LL) {ll fa = fac(a, p, p), fb = fac(b, p, p), fc = fac(a - b, p, p);ll pk = qpow(p, pot(a, p) - pot(b, p) - pot(a - b, p), p);m[tot] = p, r[tot ++] = mul(pk, mul(fa, mul(inv(fb, p), inv(fc, p), p), p), p);}return CRT();
}ll a, b, p, ans;int main() {scanf("%lld%lld%lld", &a, &b, &p);ans = exLucas(a, b, p);printf("%lld", ans);return 0;
}

一些优化（LibreOJ #181. 二项式系数）

上面给的代码只能过一些数据比较水的单组测试样例，因为是板子题所以没有优化。

要过一些其它的题比如LibreOJ #181. 二项式系数，这种数据比较毒瘤的是肯定不行的。

以下优化基于模数确定，多组输入的情况。

预处理部分

在求分离阶乘的时候，我们其实可以通过前缀和预处理出，前k个数中与pαp^{\alpha}pα互素的数的乘积sum[k]sum[k]sum[k]，以及basebasebase的λ\lambdaλ次幂base_sum[λ]base\_sum[\lambda]base_sum[λ]。

由于basebasebase与pαp^{\alpha}pα互素，通过欧拉定理可知baseφ(pα)≡1(modpα)base^{\varphi(p^{\alpha})}\equiv 1\pmod{p^{\alpha}}baseφ(pα)≡1(modpα)。

φ(pk)=pk−1φ(p)=pk−1(p−1)\varphi(p^k)=p^{k-1}\varphi(p)=p^{k-1}(p-1)φ(pk)=pk−1φ(p)=pk−1(p−1)

所以预处理时，sumsumsum预处理到pαp^{\alpha}pα，base_sumbase\_sumbase_sum预处理到pk−1(p−1)p^{k-1}(p-1)pk−1(p−1)即可。

再一个是CRT函数。

设m1,m2,…,mkm_1,m_2,\dots,m_km1,m2,…,mk为kkk个两两互素的正整数，m=∏mim=\prod{m_i}m=∏mi
令m=miMim=m_iM_im=miMi，则同余方程组{x≡b1(modm1)x≡b2(modm2)⋮x≡bk(modmk)\begin{cases}x &\equiv & b_1 \pmod{m_1} \\ x &\equiv & b_2 \pmod{m_2} \\ & \vdots & \\ x &\equiv & b_k \pmod{m_k} \end{cases}⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧xxx≡≡⋮≡b1(modm1)b2(modm2)bk(modmk)
有唯一解x≡∑i=1kbiMiMi′(modm)x\equiv {\sum_{i=1}^{k}{b_iM_iM_i^{'}}}\pmod{m}x≡∑i=1kbiMiMi′(modm)，MiMi′≡1(modmi)M_iM_i^{'}\equiv1\pmod{m_i}MiMi′≡1(modmi)

由于需要多次调用CRT，我们可以直接预处理MiMi′M_iM_i^{'}MiMi′这一部分。

再预处理完这些之后，我们就来改写fac函数，使之同时求出pot。

inline int fac(long long x, long long & pot) {int res = pot = 0;while (x) {res = res * base_sum[x / pk % phi] * sum[x % pk];pot += x / p, x /= p; } return res;
}

然后是exLucas函数

inline int exLucas(long long n, long long m) {long long pa, pb, pc;int fa = fac(n, pa);int fb = fac(m, pb);int fc = fac(n - m, pc);return fa * inv(fb) * inv(fc) * qpow(p, pa - pb - pc);
}

复杂度分析

我们来算算复杂度

唯一分解加上预处理O(log⁡P∑(φ(pk)+pk+log⁡pk))O(\log{P}\sum(\varphi(p^k)+p^k+\log{p^k}))O(logP∑(φ(pk)+pk+logpk))
查询O(Tlog⁡P(log⁡p(n)+log⁡p(m)+log⁡p(n−m)+log⁡(pot)))O(T\log{P}(\log_p{(n)}+\log_p{(m)}+\log_p{(n-m)}+\log(pot)))O(TlogP(logp(n)+logp(m)+logp(n−m)+log(pot)))

总共大概1e81e81e8和一个(可能)不太大的常数，然后本地跑了3.6秒，一直TLE。。。

感觉太玄学了就去看了大佬提交的代码，然后就看到了一个很玄学的优化（取模）。

大佬的玄学优化（取模）

来自id: CCF_NOI的提交

struct FastMod {unsigned long long b, m;void in(unsigned long long x) {b = x, m = (unsigned long long)((__uint128_t(1) << 64) / x);}int operator()(unsigned long long a) {unsigned long long q = (unsigned long long)((__uint128_t(m) * a) >> 64), r = a - q * b;return r >= b ? r - b : r;}
};

个人理解：除法的效率很低，用位运算代替除法提高效率的优化（？）

这个是我自己一直在用的，比较常规。

return (a * b - (ll)((ld)a / p * b) * p + p) % p;

我搜索了一下c++中各种运算的效率

这是在这篇博客种看到的，不保证正误，但是至少知道除法很慢就对了。

移位 > 赋值 > 大小比较 > 加法 > 减法 > 乘法 > 取模 > 除法。

由于模数确定，我们处理一下再用位运算右移代替除法会比正常取模快。

AC代码

LibreOJ 提交链接加了快读，以及压行，所以有点丑，但应该不影响阅读。

#include <cstdio>
#include <algorithm>
#include <vector>
#define gc() (p1 == p2 ? (p2 = buf + fread(p1 = buf, 1, 1 << 20, stdin), p1 == p2 ? EOF : *p1++) : *p1++)
#define read() ({ register long long x = 0, f = 1; register char c = gc(); while(c < '0' || c > '9') { if (c == '-') f = -1; c = gc();} while(c >= '0' && c <= '9') x = x * 10 + (c & 15), c = gc(); f * x; })
#define pc(x) (p - puf == 1 << 20 ? (fwrite(puf, 1, 1 << 20, stdout), p = puf) : 0, *p++ = x)
#define print(x, b) ({ pt(x), pc(b); })
char buf[1 << 20], *p1, *p2, puf[1 << 20], *p = puf;
int pt(int x) { return x <= 9 ? pc(x + '0') : (pt(x / 10), pc(x % 10 + '0')); }int P;struct FastMod {unsigned long long b, m;void in(unsigned long long x) {b = x, m = (unsigned long long)((__uint128_t(1) << 64) / x);}int operator()(unsigned long long a) {unsigned long long q = (unsigned long long)((__uint128_t(m) * a) >> 64), r = a - q * b;return r >= b ? r - b : r;}
} M;struct node {int p, pk, o, phi; FastMod f;std::vector<int> sum, base_sum;void exgcd(int a, int b, int & x, int & y) { b ? (exgcd(b, a % b, y, x), y -= a / b * x) : (x = 1, y = 0); } int inv(int a) { int x, y; exgcd(a, pk, x, y); return f(x + pk); }int qpow(int x, int n, int res = 1) { for ( ; n; n >>= 1, x = f(1ll * x * x)) if (n & 1) res = f(1ll * res * x); return res; }void set(int _p, int _pk) {f.in(_pk), p = _p, pk = _pk, phi = _pk - _pk / _p, o = M(P / _pk * 1ll * inv(P / _pk));sum.resize(pk), base_sum.resize(phi), sum[0] = base_sum[0] = 1;for (register int i = 1; i < pk; i++) sum[i] = i % p ? f(1ll * sum[i-1] * i) : sum[i-1];for (register int i = 1; i < phi; i++) base_sum[i] = f(1ll * base_sum[i-1] * sum[pk-1]);}int fac(long long x, long long & pot, int res = 1) {pot = 0;while (x) res = f(f(base_sum[x / pk % phi] * 1ll * res) * 1ll * sum[x % pk]), pot += x / p, x /= p;return res;}int get(long long a, long long b) {long long pa, pb, pc; int fa = fac(a, pa), fb = fac(b, pb), fc = fac(a - b, pc);return f(f(f(1ll * qpow(p, pa - pb - pc) * fa) * 1ll * inv(fb)) * 1ll * inv(fc));}
} pr[11];long long n, m;int main() {register int T = read(), t = P = read(), tot = 0, ans; M.in(P);for (register int i = 2; i <= t / i; i++) // Factorizationif (t % i == 0) {for (ans = 1; t % i == 0; ans *= i, t /= i) continue;pr[++ tot].set(i, ans);}if (t ^ 1) pr[++ tot].set(t, t);while (T -- ) {n = read(), m = read(), ans = 0;for (register int i = 1; i <= tot; i++)ans = M(ans + M(pr[i].get(n, m) * 1ll * pr[i].o));print(ans, '\n');}fwrite(puf, 1, p - puf, stdout);
}