CRT and exCRT

📅 2026/7/9 14:46:58 👁️ 阅读次数 📝 编程学习
CRT and exCRT

\(\text{luogu-1495}\)

给定以下方程组:

\[\left\{\begin{matrix} x \equiv b_1 \pmod{a_1} \\ x \equiv b_2 \pmod{a_2} \\ \vdots \\ x \equiv b_n \pmod{a_n} \end{matrix}\right. \]

其中 \(a_i\) 两两互质,求出满足该方程组最小的非负解 \(x\)

\(1 \le n \le 10\)\(0 \le b_i < a_i \le 10^5\)\(1 \le \prod a_i \le 10^{18}\)


CRT 的模板题。

有一个非常重要的性质是 \(a_i\) 两两互质,我们可以考虑构造出一个解。

\(M = \prod a_i\)\(m_i = \frac{M}{a_i}\),那么此时 \(\gcd(m_i, a_i) = 1\)

也就是说 \(m_i^{-1}\) 存在,那么我们可以构造出一个 \(x\)

\[x = \sum_{i=1}^n b_i \cdot m_i \cdot m_i^{-1} \mod M \]

对于 \(i \ne j\) 时,满足 \(a_i \mid m_j\),所以对于所有的 \(i\) 都满足:

\[x \equiv b_i \cdot m_i \cdot m_i^{-1} \pmod{a_i} \]

由于 \(m_i \cdot m_i^{-1} \equiv 1 \pmod{a_i}\),故满足:

\[x \equiv b_i \pmod{a_i} \]

因此,\(x\) 满足所有的同余方程,可以证明这个解在模 \(M\) 下是唯一的。


一些实现细节:

  • 运算过程中可能会超过 \(10^{18}\),用 __int128 或者用龟速乘。
  • 由于 \(a_i\) 不一定为质数,要用 exgcd 计算 \(m_i\) 的逆元。
#include<bits/stdc++.h>
#define db(x) cout << "#" << (#x) << ": " << (x) << "\n"
#define rep(i, a, b) for(int i = (a); i <= (b); ++ i)
#define per(i, a, b) for(int i = (a); i >= (b); -- i)
#define ll long long
#define pii pair<ll, ll>
#define fi first
#define se second
#define pb push_back
using namespace std;
#define N 15
#define i128 __int128ll n, a[N], b[N];void exgcd(ll a, ll b, ll &x, ll &y) {if(!b) { x = 1, y = 0; return; }exgcd(b, a % b, y, x), y -= a / b * x;
}ll inv(ll a, ll p) {ll x, y; exgcd(a, p, x, y);return (x + p) % p;
}void solution() {cin >> n; ll M = 1, res = 0;rep(i, 1, n) cin >> a[i] >> b[i], M *= a[i];rep(i, 1, n) {i128 m = M / a[i], ivm = inv(m, a[i]);res = (res + b[i] * m % M * ivm % M) % M;}cout << res << "\n";
}int main() {ios::sync_with_stdio(0), cin.tie(0), cout.tie(0);ll tas = 1; // cin >> tas;while(tas --) solution();return 0;
}

\(\text{luogu-4777}\)

给定 \(n\) 组非负整数 \(a_i, b_i\) ,求解关于 \(x\) 的方程组的最小非负整数解。

\[\left\{\begin{matrix} x \equiv b_1 \pmod{a_1} \\ x \equiv b_2 \pmod{a_2} \\ \vdots \\ x \equiv b_n \pmod{a_n} \end{matrix}\right. \]

\(1 \le n \le {10}^5\)\(1 \le a_i \le {10}^{12}\)\(0\leq b_i \leq 10^{12}\),保证所有 \(a_i\) 的最小公倍数不超过 \({10}^{18}\)

请注意程序运行过程中进行乘法运算时结果可能有溢出的风险。

数据保证有解。


exCRT 模板题。

这道题中不保证 \(a_i\) 两两互质,我们就无法构造出一个解了。

考虑每次把式子的前两项合并,最终合并成一个式子,即可得到最小非负解。

下面考虑如何将两个式子合并,假定要合并以下式子:

\[\left\{\begin{matrix} x \equiv b_1 \pmod{a_1} \\ x \equiv b_2 \pmod{a_2} \end{matrix}\right. \]

\(x = k_1 a_1 + b_1 = k_2 a_2 + b_2\),移项得:

\[k_1 a_1 - k_2 a_2 = b_2 - b_1 \]

根据裴蜀定理,当 \(\gcd(a_1, a_2) \nmid b_2 - b_1\) 时,无解。

考虑用 exgcd 求出下式的一组特解 \(k_1',k_2'\)

\[k_1 a_1 - k_2 a_2 = \gcd(a_1, a_2) \]

可得 \(k_1\) 的通解为:

\[k_1 = k_1' \frac{b_2 - b_1}{\gcd(a_1, a_2)} + t \frac{a_2}{\gcd(a_1, a_2)}, t \in \mathbb{Z} \]

带入可得:

\[\begin{align} x &= \left(k_1' \frac{b_2 - b_1}{\gcd(a_1, a_2)} + t \frac{a_2}{\gcd(a_1, a_2)}\right) a_1 + b_1 \\ &= k_1' \frac{b_2 - b_1}{\gcd(a_1, a_2)} a_1 + t \frac{a_1 a_2}{\gcd(a_1, a_2)} + b_1 \\ &= k_1' \frac{b_2 - b_1}{\gcd(a_1, a_2)} a_1 + t \cdot \text{lcm}(a_1, a_2) + b_1 \\ &\equiv k_1' \frac{b_2 - b_1}{\gcd(a_1, a_2)} a_1 + b_1 \pmod{\text{lcm}(a_1, a_2)} \end{align} \]

这样就把两个式子合并为一个同余式了,合并 \(n-1\) 次即可得到最小非负解。

#include<bits/stdc++.h>
#define db(x) cout << "#" << (#x) << ": " << (x) << "\n"
#define rep(i, a, b) for(int i = (a); i <= (b); ++ i)
#define per(i, a, b) for(int i = (a); i >= (b); -- i)
#define ll long long
#define pii pair<ll, ll>
#define fi first
#define se second
#define pb push_back
using namespace std;
#define N 100005
#define i128 __int128ll n, a[N], b[N];ll exgcd(ll a, ll b, ll &x, ll &y) {if(!b) { x = 1, y = 0; return a; }ll d = exgcd(b, a % b, y, x); y -= a / b * x; return d;
}ll exCRT() {ll res = b[1], M = a[1];rep(i, 2, n) {ll x, y, d = exgcd(M, a[i], x, y); if((b[i] - res) % d) return -1;ll a1 = M; M = M / d * a[i];ll k1 = (i128)x * ((b[i] - res) / d) % M;res = ((i128)k1 * a1 % M + res) % M;}return (res + M) % M;
}void solution() {cin >> n;rep(i, 1, n) cin >> a[i] >> b[i];cout << exCRT() << "\n";
}int main() {ios::sync_with_stdio(0), cin.tie(0), cout.tie(0);ll tas = 1; // cin >> tas;while(tas --) solution();return 0;
}

\(\text{luogu-2480}\)

求解下式的值:

\[g^{\sum\limits_{k \mid n} \binom{n}{k}} \mod 999911659 \]

\(1 \le n, g \le 10^9\)


由于 \(n\) 太大了,指数部分直接爆了,显然需要对指数取个模。

\(P = 999911659\),由于:

\[g^{P-1} \equiv 1 \equiv g^0 \mod P \]

故:

\[g^x \equiv g^{x \bmod (P - 1)} \mod P \]

因此原式可化为:

\[g^{\sum\limits_{k \mid n} \binom{n}{k} \mod (P - 1)} \mod P \]

但是这个指数还是很难搞,模数太大,也没办法用 Lucas 快速计算 \(\binom{n}{k}\)

既然因为模数太大没办法算,那么就把模数拆成小模数就好了。

因为 \(999911658 = 2 \times 3 \times 4679 \times 35617\),拆成这四个模数即可,枚举因数再用 Lucas 求出对这四个模数取模后的余数,记为 \(a_1,a_2,a_3,a_4\),令 \(x = \sum\limits_{k \mid n} \binom{n}{k}\),则可以得到下式:

\[\left\{\begin{matrix} \begin{array}{l} x \equiv a_1 \pmod{2} \\ x \equiv a_2 \pmod{3} \\ x \equiv a_3 \pmod{4679} \\ x \equiv a_4 \pmod{35617} \end{array} \end{matrix}\right. \]

用 CRT 就可以求出 \(x\) 的最小非负解,也就是 \(\sum\limits_{k \mid n} \binom{n}{k} \mod (P - 1)\) 的值。

然后直接用快速幂计算 \(g^x \bmod P\) 即可,注意需要特判 \(P \mid g\) 的情况,输出 \(0\)

#include<bits/stdc++.h>
#define db(x) cout << "#" << (#x) << ": " << (x) << "\n"
#define rep(i, a, b) for(int i = (a); i <= (b); ++ i)
#define per(i, a, b) for(int i = (a); i >= (b); -- i)
#define ll long long
#define pii pair<ll, ll>
#define fi first
#define se second
#define pb push_back
using namespace std;
#define N 40005
#define P 999911659ll n, g, a[5], b[5], inv[5][N], p[5][N], ivp[5][N];void init() {a[1] = 2, a[2] = 3, a[3] = 4679, a[4] = 35617;rep(i, 1, 4) p[i][0] = p[i][1] = inv[i][1] = 1,ivp[i][0] = ivp[i][1] = 1;rep(i, 2, N - 5) rep(j, 1, 4)p[j][i] = p[j][i - 1] * i % a[j], inv[j][i] = a[j] - a[j] / i * inv[j][a[j] % i] % a[j],ivp[j][i] = ivp[j][i - 1] * inv[j][i] % a[j];
}ll C(ll n, ll m, ll k) { if(n < m) return 0;return p[k][n] * ivp[k][n - m] % a[k] * ivp[k][m] % a[k]; 
}ll lucas(ll n, ll m, ll p) {if(!m) return 1;return C(n % a[p], m % a[p], p) * lucas(n / a[p], m / a[p], p) % a[p];
}ll qpow(ll a, ll b, ll p) {ll res = 1;while(b) {if(b & 1) res = res * a % p;a = a * a % p, b >>= 1;}return res;
}ll CRT() {ll M = P - 1, res = 0;rep(i, 1, 4) {ll m = M / a[i], ivm = qpow(m, a[i] - 2, M);res = (res + b[i] * m % M * ivm % M) % M;}return res;
}void solution() {cin >> n >> g; init();if(!(g % P)) { cout << "0\n"; return; }rep(i, 1, sqrt(n)) if(!(n % i)) rep(j, 1, 4) {b[j] = (b[j] + lucas(n, i, j)) % a[j];if(i * i < n) b[j] = (b[j] + lucas(n, n / i, j)) % a[j];}cout << qpow(g, CRT(), P) << "\n";
}int main() {ios::sync_with_stdio(0), cin.tie(0), cout.tie(0);ll tas = 1; // cin >> tas;while(tas --) solution();return 0;
}

\(\text{luogu-4774}\)

\(n\) 只巨龙,第 \(i\) 只的初始生命值为 \(a_i\),并且拥有一个生命恢复值 \(p_i\)。你有 \(m\) 把剑,每把剑有一个攻击力。现在你要用这 \(m\) 把剑依次杀死这 \(n\) 只巨龙。

在攻击第 \(i\) 只巨龙之前,你需要选出攻击力不大于 \(a_i\) 且攻击力最大的一把剑。如果这样的剑不存在,则选出攻击力最小的一把剑。设这把剑的攻击力为 \(k\),接下来:

  • 你会连续攻击巨龙 \(x\) 次,巨龙的生命值变为 \(a_i-xk\)
  • 巨龙会不断恢复生命值,每次恢复 \(p_i\),直到生命值非负;
  • 此时,若巨龙生命值为 \(0\),巨龙死去;
  • 用来攻击的这把剑会消失。与此同时,如果巨龙死去,你会获得一把新的剑,攻击力可能和这一把不同。

\(x\) 的最小非负解,或者指出 \(x\) 不存在并输出 \(-1\)

\(1 \le n, m \le 10^5\)\(p_i \ge 1\)\(a_i,\text{lcm } p_i \le 10^{12}\)


exCRT 综合练习好题。

可以用 multiset 找到攻击每只巨龙使用的剑的攻击力 \(b_i\)

对于每只巨龙都会攻击 \(x\) 次,则:

\[\left\{\begin{matrix} a_1 - b_1x \equiv 0 \pmod{p_1} \\ a_2 - b_2x \equiv 0 \pmod{p_2} \\ \vdots \\ a_n - b_nx \equiv 0 \pmod{p_n} \end{matrix}\right. \]

移向得:

\[\left\{\begin{matrix} b_1x \equiv a_1 \pmod{p_1} \\ b_2x \equiv a_2 \pmod{p_2} \\ \vdots \\ b_nx \equiv a_n \pmod{p_n} \end{matrix}\right. \]

由于有系数,不能直接用 exCRT,考虑把系数处理掉,有下式:

\[b_i x + p_i y = a_i \]

根据裴蜀定理,当 \(\gcd(b_i, p_i) \nmid a_i\) 时,无解。

用 exgcd 求出下式的一组特解 \(x', y'\)

\[b_i x + p_i y = \gcd(b_i, p_i) \]

可得 \(x\) 的通解为:

\[x = x' \frac{a_i}{\gcd(b_i, p_i)} + t \frac{p_i}{\gcd(b_i, p_i)}, t \in \mathbb{Z} \]

因此有:

\[x \equiv x' \frac{a_i}{\gcd(b_i, p_i)} \pmod{\frac{p_i}{\gcd(b_i, p_i)}} \]

这样就把系数处理掉了,剩下的用 exCRT 做就好了。

注意,exCRT 求出来的最小非负解 \(x\) 不一定可以把所有巨龙的生命值都变为非正,要不断加上 \(\text{lcm} \left(\frac{p_i}{\gcd(b_i, p_i)}\right)\),直到 \(x\) 可以使所有巨龙的生命值变为非正。

#include<bits/stdc++.h>
#define db(x) cout << "#" << (#x) << ": " << (x) << "\n"
#define rep(i, a, b) for(int i = a; i <= b; ++ i)
#define per(i, a, b) for(int i = a; i >= b; -- i)
#define ll long long
#define pii pair<ll, ll>
#define fi first
#define se second
#define pb push_back
using namespace std;
#define N 100005
#define i128 __int128ll n, m, a[N], b[N], p[N], t[N], c[N], d[N], mx;
multiset<ll> s;ll exgcd(ll a, ll b, ll &x, ll &y) {if(!b) { x = 1, y = 0; return a; }ll g = exgcd(b, a % b, y, x);y -= a / b * x; return g;
}ll exCRT() {ll M = d[1], res = c[1];rep(i, 2, n) {ll x, y, g = exgcd(M, d[i], x, y);if((c[i] - res) % g) return -1;ll d1 = M; M = M / g * d[i];ll k1 = (i128)x * ((c[i] - res) / g) % M;res = ((i128)k1 * d1 % M + res) % M;}res = (res + M) % M;if(res < mx) res += (mx - res + M - 1) / M * M;return res;
}void solution() {cin >> n >> m; s.clear(), mx = 0;rep(i, 1, n) cin >> a[i];rep(i, 1, n) cin >> p[i];rep(i, 1, n) cin >> t[i];rep(i, 1, m) { ll x; cin >> x; s.insert(x); }rep(i, 1, n) {auto it = s.upper_bound(a[i]);if(it != s.begin()) it --;b[i] = *it, s.erase(it), s.insert(t[i]);mx = max(mx, (a[i] + b[i] - 1) / b[i]);}rep(i, 1, n) {ll x, y, g = exgcd(b[i], p[i], x, y);if(a[i] % g) { cout << "-1\n"; return; }d[i] = p[i] / g, c[i] = (i128)x * (a[i] / g) % d[i];}cout << exCRT() << "\n";
}int main() {ios::sync_with_stdio(0), cin.tie(0), cout.tie(0);ll tas = 1; cin >> tas;while(tas --) solution();return 0;
}