求解模线性方程组

问题描述

给定了 \(n\) 组除数 \(m_i\) 和余数 \(r_i\) ，通过这 \(n\) 组
\((m_i,r_i)\) 求解一个 \(x\) ，使得 \(x \bmod m_i = r_i\)
这就是 模线性方程组 。

数学形式表达就是 ：

\(\begin{cases} x \equiv r_1 \pmod {m_1}\\ x \equiv r_2 \pmod {m_2}\\ \vdots \\ x \equiv r_n \pmod {m_n} \end{cases}\)

求解一个 \(x\) 满足上列所有方程。

算法解决

中国剩余定理

中国剩余定理提供了一个较为通用的解决方法。（似乎是唯一一个以中国来命名的定理）

如果 \(m_1, m_2, \dots , m_n\) 两两互质，则对于任意的整数
\(r_1, r_2 , \dots , r_n\) 方程组 \((S)\)
有解，并且可以用如下方法来构造：

令 \(\displaystyle M = \prod_{i=1} ^ {n} m_i\) ，并设
\(\displaystyle M_i = \frac{M}{m_i}\) ，令 \(t_i\) 为
\(M_i\) 在模 \(m_i\) 意义下的逆元（也就是 \(t_i \times M_i \equiv 1 \pmod {m_i}\) ）。

不难发现这个 \(t_i\)
是一定存在的，因为由于前面要满足两两互质，那么 \(M_i \bot m_i\)
，所以必有逆元。

那么在模 \(M\) 意义下，方程 \((S)\) 有且仅有一个解 ：
\(\displaystyle x = (\sum_{i=1}^{n} r_i t_i M_i) \bmod M\)
。

不难发现这个特解是一定满足每一个方程的，只有一个解的证明可以考虑特解
\(x\) 对于第 \(i\) 个方程，下个继续满足条件的 \(x\) 为 \(x + m_i\) 。那么对于整体来说，下个满足条件的数就是 \(x + \mathrm{lcm} (m_1, m_2, \dots, m_n) = x + M\) 。

严谨数学证明请见 百度百科。

利用扩欧合并方程组

我们考虑合并方程组，比如从 \(n=2\) 开始递推。

\(\begin{cases} x \equiv r_1 \pmod {m_1}\\ x \equiv r_2 \pmod {m_2} \end{cases}\)

也就等价于

\(\begin{cases} x = m_1 \times k_1 + r_1 \\ x = m_2 \times k_2 + r_2 \\ \end{cases}\)

此处 \(k_1, k_2 \in \mathbb{N}\) 。联立后就得到了如下一个方程：

\(\begin{aligned} m_1 \times k_1 + r_1 &= m_2 \times k_2 + r_2 \\ m_1 \times k_1 - m_2 \times k_2 &= r_2 - r_1 \end{aligned}\)

我们令 \(a = m_1, b = m_2, c = r_2 - r_1\) 就变成了 \(ax+by=c\)
的形式，用之前讲过的 扩展欧几里得 ，可以求解。

首先先判断有无解。假设存在解，我们先解出一组 \((k_1,k_2)\)
，然后带入解出 \(x=x_0\) 的一个特解。

我们将这个可以扩展成一个解系：

\(x = x_0 + k \times \mathrm{lcm} (m_1,m_2) , k \in \mathbb{N}\)

由于前面不定方程的结论， \(k_1\) 与其相邻解的间距为
\(\displaystyle \frac{m_2}{\gcd(m_1,m_2)}\) ，又有 \(x=m_1 \times k_1 + r_1\) 。所以 \(x\) 与其相邻解的距离为
\(\displaystyle \frac{m_1 m_2}{\gcd(m_1,m_2)} = \mathrm{lcm}(m_1,m_2)\) 。

所以我们令 \(M=\mathrm{lcm}(m_1, m_2), R = x_0\) 则又有新的模方程
\(x \equiv R \pmod M\) 。

然后我们就将两个方程合并成一个了，只要不断重复这个过程就能做完了。

这个比 中国剩余定理 优秀的地方就在于它这个不需要要求 \(m\)
两两互质，并且可以较容易地判断无解的情况。

代码实现

注意有些细节，比如求两个数的 \(\mathrm{gcd}\) 的时候，一定要先除再乘，防止溢出！！

ll mod[N], rest[N];
ll Solve() {
    For (i, 1, n - 1) {
        ll a = mod[i], b = mod[i + 1], c = rest[i + 1] - rest[i], gcd = __gcd(a, b), k1, k2;
        if (c % gcd) return - 1;
        a /= gcd; b /= gcd; c /= gcd;
        Exgcd(a, b, k1, k2);

        k1 = ((k1 * c) % b + b) % b;
        mod[i + 1] = mod[i] / __gcd(mod[i], mod[i + 1]) * mod[i + 1] ;
        rest[i + 1] = (mod[i] * k1 % mod[i + 1] + rest[i]) % mod[i + 1];
    }
    return rest[n];
}