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

中国剩余定理

《孙子算经》有云：

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

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

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

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

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

此时，因为，其他同理。

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

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

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

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 inv(ll x, ll p) {
    ll ans, tmp;
    exgcd(x, p, ans, tmp);
    return (ans % p + p) % p;
}

ll CRT() {
    ll res = 0;
    for (int i = 1; i <= n; i++) M *= m[i];
    for (int i = 1; i <= n; i++) {
        ll Mi = M / m[i];
        res += a[i] * Mi * inv(Mi, m[i]);
    }
    return res % M;
}

扩展中国剩余定理

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

还是那个方程组：

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

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

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

也即：

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

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

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

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

那么对于方程

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

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

做差可得：

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

总结一下，对于方程组：

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

#include <bits/stdc++.h>

#define N 100010
using namespace std;

typedef long long ll;

struct Equation {
    ll a, m;
};

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

ll exgcd(ll a, ll b, ll &x, ll &y) {
    // 扩展欧几里得过程，求解特解 t1 和 t2
    if (!b) {
        x = 1, y = 0;
        return a;
    }
    ll d = exgcd(b, a % b, y, x);
    y -= a / b * x;
    return d;
}

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

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

        ll t1, t2;
        ll D = gcd(cur.m, nxt.m);
        ll S = lcm(cur.m, nxt.m);
        if ((nxt.a - cur.a) % D != 0) return -1; // 无解判断
        exgcd(cur.m, nxt.m, t1, t2); // 求解特解
        ll X = ((mul(mul(cur.m, t1, S), (nxt.a - cur.a) / D, S)) + cur.a) % S;
        stk.push((Equation) {X, S});
    }
    Equation res = stk.top();
    return (res.a % res.m + res.m) % res.m; // 转为正整数
}

int main() {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);
    cout.tie(nullptr);

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

    return 0;
}