数论补完计划 Part4 中国剩余定理

中国剩余定理#

《孙子算经》有云：

今有物不知其数，三三数之剩二，五五数之剩三，七七数之剩二，问物几何？

大意为：有个物品，满足如下线性同余方程组：

刚入坑编程的估计在“循环结构”那一章里就写过求解这个方程的枚举代码了。这种方式不好的一点就是码量随方程数的增多而膨胀（~~但是你可以写一个生成计算代码的程序，然后运行编译后的程序~~）。我们急需一种通解，在码量不增的情况下实现对这个方程组的求解。

我们考虑拆开这个方程组，分别令满足：

两个方程组的写法其实是等价的，只是把同余转化为了带余除法罢了。我们惊奇地发现，如果再加上几个约束条件，原方程组的解就可以经过表示出来：

此时，因为，其他同理。

而我们的中国剩余定理断言道，对于模数均互质的元线性同余方程组：

它的解，其中。因此只需要计算同余号右侧的部分即可，由于方程组模数互质并不等于每个方程的模数是质数，因此需要用到扩展欧几里得算法来求解逆元。

1
int n;
2
ll m[N], a[N];
3
ll M = 1;
4

5
ll exgcd(ll a, ll b, ll &x, ll &y) {
6
    if (!b) {
7
        x = 1, y = 0;
8
        return a;
9
    }
10
    ll d = exgcd(b, a % b, y, x);
11
    y -= a / b * x;
12
    return d;
13
}
14

15
ll inv(ll x, ll p) {
16
    ll ans, tmp;
17
    exgcd(x, p, ans, tmp);
18
    return (ans % p + p) % p;
19
}
20

21
ll CRT() {
22
    ll res = 0;
23
    for (int i = 1; i <= n; i++) M *= m[i];
24
    for (int i = 1; i <= n; i++) {
25
        ll Mi = M / m[i];
26
        res += a[i] * Mi * inv(Mi, m[i]);
27
    }
28
    return res % M;
29
}

扩展中国剩余定理#

经历过数论大风大浪的应该都知道，马上要开始推导模数不全互质的情况了。

还是那个方程组：

此时不保证两两互质，求出方程的最小正整数解。此时因为失去了互质，构造通解的方式已经不能用了。因此我们转而寻找另外的方法来简化这个方程组。

对大量事实的观测表明，两个线性同余方程在合并后遵循如下性质：

模数变为原两个方程模数的最小公倍数
合并得到的方程与原方程同型
合并得到的方程不一定有解

也即：

事实上，假设我们有方程和。那么转化成带余方程之后就有：

观察到这是二元一次不定方程的形式，根据贝祖定理，有解的充要条件是，因此无解情况就可以筛去了。

约去最大公倍数得：。用扩展欧几里得解出方程的特解，那么可以用如下关系得到：

于是。用去代结果是一样的。

那么对于方程

仅当时有整数解。那么方程的通解可以表示成的形式，这样就能解释合并后方程的模数为什么为最小公倍数了。

我们假设，并假设存在关系：

做差可得：

得到关系。注意到两个数的最小公倍数一定满足。但是对于非负整数来说，是一定成立的，且在时取得等号，那么当且仅当在时才存在最小公倍数的整除关系。因此导出式是一定正确的。

总结一下，对于方程组：

有解的充要条件是。若已知有解，假设二元一次不定方程的特解是，并记，那么原方程等价于：

1
#include <bits/stdc++.h>
2

3
#define N 100010
4
using namespace std;
5

6
typedef long long ll;
7

8
struct Equation {
9
    ll a, m;
10
};
11

12
ll a[N], b[N];
13
stack <Equation> stk;
14

15
ll exgcd(ll a, ll b, ll &x, ll &y) {
16
    // 扩展欧几里得过程，求解特解 t1 和 t2
17
    if (!b) {
18
        x = 1, y = 0;
19
        return a;
20
    }
21
    ll d = exgcd(b, a % b, y, x);
22
    y -= a / b * x;
23
    return d;
24
}
25

26
ll mul(ll a, ll b, ll p) {
27
    // 著名的龟速乘，就是把快速幂中的乘号改成加号，用来在不爆范围的情况下计算乘积取模的结果
28
    ll res = 0;
29
    bool flag = (b < 0) ^ (a < 0); // 处理负数情况，同号相乘为正，异号为负
30
    a = abs(a) % p;
31
    b = abs(b) % p;
32
    // 注意到：a * (-b) = -(a * b)，数字按绝对值相乘，负号最后处理
33
    while (b) {
34
        if (b & 1) res = (res + a) % p;
35
        a = (a + a) % p;
36
        b >>= 1;
37
    }
38
    return res * (flag ? -1 : 1) % p; // 处理负号
39
}
40

41
ll exCRT() {
42
    while (stk.size() > 1) {
43
        // 不断合并方程直至剩下最后一个
44
        Equation cur = stk.top();
45
        stk.pop();
46
        Equation nxt = stk.top();
47
        stk.pop();
48

49
        ll t1, t2;
50
        ll D = gcd(cur.m, nxt.m);
51
        ll S = lcm(cur.m, nxt.m);
52
        if ((nxt.a - cur.a) % D != 0) return -1; // 无解判断
53
        exgcd(cur.m, nxt.m, t1, t2); // 求解特解
54
        ll X = ((mul(mul(cur.m, t1, S), (nxt.a - cur.a) / D, S)) + cur.a) % S;
55
        stk.push((Equation) {X, S});
56
    }
57
    Equation res = stk.top();
58
    return (res.a % res.m + res.m) % res.m; // 转为正整数
59
}
60

61
int main() {
62
    ios::sync_with_stdio(false);
63
    cin.tie(nullptr);
64
    cout.tie(nullptr);
65

66
    int n;
67
    cin >> n;
68
    for (int i = 1; i <= n; i++) cin >> a[i] >> b[i];
69
    for (int i = 1; i <= n; i++) stk.push((Equation) {b[i], a[i]});
70
    cout << exCRT() << endl;
71

72
    return 0;
73
}

音乐

JustPureH₂O

穷方圆平直之情，尽规矩准绳之用

音乐

数论补完计划 Part4 中国剩余定理

中国剩余定理#

扩展中国剩余定理#

支持与分享

评论区

音乐

目录

音乐

JustPureH2O

穷方圆平直之情，尽规矩准绳之用

音乐

数论补完计划 Part4 中国剩余定理

中国剩余定理#

扩展中国剩余定理#

支持与分享

评论区

音乐

目录

JustPureH₂O