求解模線性方程組

問題描述

給定了 \(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];
}

本文來自部落格園，作者：蒟蒻orz，轉載請註明原文連結：https://www.cnblogs.com/orzz/p/15583848.html