有物不知其数,三三数之剩二,五五数之剩三,七七数之剩二。问物几何?
该问题出自《孙子算经》,具体问题的解答口诀由明朝数学家程大位在《算法统宗》中给出:
三人同行七十希,五树梅花廿一支,七子团圆正半月,除百零五便得知。
\(2 \times 70 + 3 \times 21 + 2 \times 15 = 233 = 2 \times 105 + 23\),故答案为 \(23\)。
定义
中国剩余定理 (Chinese Remainder Theorem, CRT) 可求解如下形式的一元线性同余方程组(其中 \(n_1, n_2, \cdots, n_k\) 两两互质):
\[\begin{cases}
x \bmod p_1 = a_1\\
x \bmod p_2 = a_2\\
\cdots\\
x \bmod p_n = a_n\\
\end{cases}
\]
过程
对于方程
\[\begin{cases}
x \bmod p_1 = a_1\\
x \bmod p_2 = a_2\\
\end{cases}
\]
我们可以得知 \(x = k_1 p_1 + a_1 = k_2 p_2 + a_2\),进而推出 \(k_1 p_1 - k_2 p_2 = a_2 - a_1\)
我们观察这个式子,是不是跟 \(xa + yg = g\) 很像?
我们可以用扩展欧几里得算法来做。
扩展欧几里得算法代码
int exgcd(int a, int b, int &x, int &y) {
if (!b) {
x = 1, y = 0;
return a;
}
int xp, yp;
int g = exgcd(b, a % b, xp, yp);
x = yp, y = xp - yp * (a / b);
return g;
}
中国剩余定理代码
void solve(int p1, int a1, int p2, int a2, int &p, int &a) {
// x % p1 = a1
// x % p2 = a2
// x % p = a
// x = k1 * p1 + a1 = k2 * p2 + a2
// k1 * p1 - k2 * p2 = a2 - a1
int g, k1, k2;
g = exgcd(p1, p2, k1, k2);
// k1 * p1 + k2 * p2 = g
k2 = -k2;
// k1 * p1 - k2 * p2 = g
if ((a2 - a1) % g != 0) {
p = -1, a = -1;
}
else {
int k = (a2 - a1) / g;
k1 = k1 * k, k2 = k2 * k;
// k1 * p1 - k2 * p2 = a2 - a1
int x = k1 * p1 + a1;
p = p1 / g * p2;
a = (x % p + p) % p;
}
}
通过观察,你会发现,这种解法不需要 \(p_1 \perp p_2\),因此扩展中国剩余定理也是这个解法,这是一种通解。
题目
【模板】扩展中国剩余定理(EXCRT)
卡 __int128。
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N = 1e5 + 5;
int n;
ll a[N], b[N];
ll qtimes(ll x, ll y, ll p) {
if (y == 0) return 0;
ll z = qtimes(x, y / 2, p);
z = (z + z) % p;
if (y & 1) {
z = (1ll * z + x) % p;
}
return z;
}
ll exgcd(ll a, ll b, ll &x, ll &y) {
if (b == 0) {
x = 1, y = 0;
return a;
}
ll xp, yp;
ll g = exgcd(b, a % b, xp, yp);
x = yp, y = xp - yp * (a / b);
return g;
}
void CRT(ll p1, ll a1, ll p2, ll a2, ll &p, ll &a) {
ll g, k1, k2;
g = exgcd(p1, p2, k1, k2);
k2 = -k2;
if ((a2 - a1) % g != 0) {
p = -1, a = -1;
}
else {
ll k = (a2 - a1) / g;
__int128 K1 = k1 * k, K2 = k2 * k;
__int128 x = K1 * p1 + a1;
p = p1 / g * p2;
a = (x % p + p) % p;
}
}
int main() {
ios::sync_with_stdio(false);
cin.tie(0), cout.tie(0);
cin >> n;
for (int i = 1; i <= n; ++ i) {
cin >> a[i] >> b[i];
}
ll mod = a[1], rest = b[1];
for (int i = 2; i <= n; ++ i) {
CRT(mod, rest, a[i], b[i], mod, rest);
}
cout << rest % mod << '\n';
return 0;
}
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N = 1e5 + 5;
int n;
ll a[N], b[N];
ll qtimes(ll x, ll y, ll p) {
if (y == 0) return 0;
ll z = qtimes(x, y / 2, p);
z = (z + z) % p;
if (y & 1) {
z = (1ll * z + x) % p;
}
return z;
}
ll exgcd(ll a, ll b, ll &x, ll &y) {
if (b == 0) {
x = 1, y = 0;
return a;
}
ll xp, yp;
ll g = exgcd(b, a % b, xp, yp);
x = yp, y = xp - yp * (a / b);
return g;
}
void CRT(ll p1, ll a1, ll p2, ll a2, ll &p, ll &a) {
ll g, k1, k2;
g = exgcd(p1, p2, k1, k2);
k2 = -k2;
if ((a2 - a1) % g != 0) {
p = -1, a = -1;
}
else {
ll k = (a2 - a1) / g;
__int128 K1 = k1 * k, K2 = k2 * k;
__int128 x = K1 * p1 + a1;
p = p1 / g * p2;
a = (x % p + p) % p;
}
}
int main() {
ios::sync_with_stdio(false);
cin.tie(0), cout.tie(0);
cin >> n;
for (int i = 1; i <= n; ++ i) {
cin >> a[i] >> b[i];
}
ll mod = a[1], rest = b[1];
for (int i = 2; i <= n; ++ i) {
CRT(mod, rest, a[i], b[i], mod, rest);
}
cout << rest % mod << '\n';
return 0;
}