经典网文分享: 初等数论Part 2：中国剩余定理

1.这是一个初等数论的入门级别文章

2.适合中学数学水平的读者

3.建议在阅读此篇之前阅读Part1

4.主要内容：中国剩余定理（CRT），贝祖定理，扩展欧几里得算法，逆元

《九章算术》中曾经提到过一个经典的问题

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

翻译一下就是：已知一个正整数模3余2，模5余3，模7余2，求这个数是几？

写成数学语言就是求解同余方程组

$x\equiv 2\pmod{3}\\ x\equiv 3\pmod{5}\\ x\equiv 2\pmod{7}$

在Part1中我们曾经提到过一下，解这种同余方程组的时候我们往往使用中国剩余定理

已知正整数 $n_1, \dots, n_k>1$ 两两互质，定义 $N=\prod_{i=1}^k n_i$ 和 $N_i=\frac{N}{n_i}$
又已知整数 $a_1, \dots, a_k$ ，那么同余方程组
$x\equiv a_1 \pmod{n_1}\\ x\equiv a_2 \pmod{n_2}\\ \vdots\\ x\equiv a_k\pmod{n_k}$
在 $0\le x < N$ 范围内有且仅有一个解
因为 $\gcd(n_i, N_i)=1\quad \forall i\in\{1, 2, \dots, k\}$ ，由贝祖定理我们能够得到对于每一个 $i\in\{1, \dots, k\}$ 存在整数 $M_i, m_i$ 使得
$M_iN_i+m_in_i=1\\$
通过扩展欧几里得算法求出 $M_1, \dots, M_k$ ，这个同余方程组的一个特解就可以表示为
$x_0=\sum_{i=1}^k a_i M_iN_i \\$
通解可以表示为
$x=Nq+x_0\\$
其中 $q\in \mathbb{Z}$ 是任意整数

这段的信息量非常大，我们下面逐步来理解。

首先补充一点必要的背景知识

1.贝祖定理
对于整数 $a, b$ 存在整数 $x, y$ 使得 $\gcd(a, b)=ax+by$
整数 $a, b$ 互质的充分必要条件是存在整数 $x, y$ 使得 $ax+by=1$

已知非零整数 $a, b$ ，记集合 $S=\{ax+by\mid x, y\in\mathbb{Z}\ \ \mathrm{and}\ \ ax+by>0\}$

记整数 $d=ax_0+by_0$ 为 $S$ 的最小元素

我们只需要证明 $d=\gcd(a, b)$ 就证明了贝祖定理，

设 $a=dq+r$ 其中 $q, r$ 是 $a$ 模 $d$ 的商和余数

$r=a-dq=a-(ax_0+by_0)q=a(1-x_0q)+b(y_0 q)\\$

所以 $r\in S\cup\{0\}$ ，我们又知道 $0\le r<d$ 而 $d$ 又是 $S$ 的最小元素，所以 $r=0$

这意味着 $d\mid a$ ，同理可证 $d\mid b$ ，又因为 $\gcd(a, b)\mid d$ ，所以 $d=\gcd(a, b)$
另外，贝祖等式的系数并不唯一，有无穷组系数 $(x, y)$ 都能够满足 $\gcd(a, b)=ax+by$ ，事实上，如果 $(x, y)$ 是一组系数，那么所有系数可以表示为
$\left(x+k\cdot\frac{b}{\gcd(a, b)} ,\ y-k\cdot\frac{a}{\gcd(a, b)}\right)\\$
其中 $k\in\mathbb{Z}$ 是任意整数

2.欧几里得算法（辗转相除法和更相止损术）

在小学可能大家有学过Division-based Euclidean algorithm（辗转相除法）找最大公因数

比如说如果我们想要找18和14的最大公因数

$\begin{align}\gcd(18, 14)&=\gcd(14, 18\bmod 14)=\gcd(14, 4)\\&=\gcd(4, 14\bmod4)=\gcd(4, 2)\\&=\gcd(2, 4\bmod 2)=\gcd(2, 0)=2\end{align}$

def gcd(a, b):     if b == 0:         return a     else:         return gcd(b, a % b)

原理是 $\gcd(a, b)=\gcd(b, a\bmod b)\quad \forall a, b\in \mathbb{Z}$

这个只需要注意到 $a\bmod b = a - b \left\lfloor \frac{a}{b}\right\rfloor$

类似的，我们有Subtraction-based Euclidean algorithm（更相减损术）

$\begin{align}\gcd(18, 14)&=\gcd(14, 18-14)=\gcd(14, 4)\\&=\gcd(4, 14-4)=\gcd(4, 10)\\&=\gcd(4, 10-4)=\gcd(4, 6)\\&=\gcd(4, 6-4)=\gcd(4, 2)\\&=\gcd(2, 4-2)=\gcd(2, 2)=2\end{align}$

def gcd(a, b):     #处理gcd(a, 0)=gcd(0, a)=a的情况     if a * b == 0: return a + b          if a == b:         return a     elif a > b:         return gcd(b, a - b)     else:         return gcd(a, b - a)

原理是 $\gcd(a, b)=\gcd(b, a-b)\quad \forall a, b \in \mathbb{Z}$

辗转相除法，更相减损术和位运算结合还可以优化最大公因数的算法，这里不深入探讨。

4.扩展欧几里得算法

中国剩余定理中有一个重要步骤就是用扩展欧几里得算法求解贝祖等式的系数。

def exgcd(a, b, x0 = 1, x1 = 0, y0 = 0, y1 = 1):     if b == 0:         return a, x0, y0     else:         q, r = divmod(a, b)         return exgcd(b, r, x1, x0 - q * x1, y1, y0 - q * y1)

算法返回的 $\mathtt{(a, x0, y0)}$ 中 $\mathtt{a}$ 是最大公因数， $\mathtt{(x0, y0)}$ 是贝祖等式的系数

而且 $\mathtt{(x0, y0)}$ 是两组数值最小的系数之一，也就是说他们满足

$|\mathtt{x0}|\le \left |\frac{b}{\gcd(a, b)}\right |\quad\quad |\mathtt{y0}|\le \left | \frac{a}{\gcd(a, b)}\right |\\$

取等号的条件是 $a\mid b$ 或者 $b\mid a$

为了证明扩展欧几里得算法的正确性，可以用以下递推式描述它

$\begin{alignat}{2} &r_{i+1}=r_{i-1}\ - r_i &&\left\lfloor \frac{r_{i-1}}{r_i}\right\rfloor\\ &x_{i+1}=x_{i-1} - x_i&&\left\lfloor \frac{r_{i-1}}{r_i}\right\rfloor\\ &y_{i+1}=y_{i-1} \ - y_i&&\left\lfloor \frac{r_{i-1}}{r_i}\right\rfloor\\ \end{alignat}\\$

初始值 $(r_0, r_1, x_0, x_1, y_0, y_1)=(a, b, 1, 0, 0, 1)$ ，可以证明，当 $(r_0, r_1)=(a, b)$ 时

$a x_i+b y_i=r_i\quad\forall i\in\mathbb{N}\\$

这个通过归纳法即可

$\begin{align}r_{i+1}&=a x_{i-1}+b y_{i-1} -(a x_i+b y_i) \left\lfloor \frac{r_{i-1}}{r_i}\right\rfloor \\&=a\left(x_{i-1} - x_i\left\lfloor \frac{r_{i-1}}{r_i}\right\rfloor\right)+b\left(y_{i-1} \ - y_i\left\lfloor \frac{r_{i-1}}{r_i}\right\rfloor\right)\\&=ax_{i+1}+by_{i+1}\end{align}\\$

容易验证 $r_0=ax_0+by_0$ 所以算法的正确性就得到证明了

我们可以从「插值」的角度来理解中国剩余定理。

首先提出这么一个问题：如何找到一个多项式方程 $f:\mathbb{R}\rightarrow\mathbb{R}$ 使得

$f(a_1)=b_1\\ f(a_2)=b_2\\ \vdots\\ f(a_k)=b_k$

其中实数 $a_1, \dots, a_k$ 和实数 $b_1, \dots, b_k$ 满足 $a_1<a_2<\cdots<a_k$

这里我们就需要拉格朗日多项式的帮助了

（推荐阅读知乎用户 @马同学的关于拉格朗日多项式的回答）

3.经过点 $(a_1, b_1), (a_2, b_2), \dots, (a_k, b_k)$ 的拉格朗日多项式函数为
$\mathcal{L}(x)=\sum_{i=1}^k b_i \delta_i(x)\\$
其中
$\delta_i(x)=\prod_{1 \le j\le k \atop j\ne i}\frac{x-a_j}{a_i-a_j}=\frac{x-a_1}{a_i-a_1}\cdots \frac{x-a_{i-1}}{a_i-a_{i-1}}\cdot \frac{x-a_{i+1}}{a_i-a_{i+1}}\cdots \frac{x-a_k}{a_i-a_k}\\$

我们容易发现 $\delta_{i}(a_j)=\left\{ \begin{matrix}1 \quad \quad i=j \\ 0 \quad\quad i\ne j\end{matrix} \right.$ ，所以代入一下

$\begin{align}\mathcal{L}(a_i)&=b_1\delta_1(a_i)+b_2\delta_2(a_i)+\cdots+b_i\delta_i(a_i)+\cdots+b_{k-1}\delta_{k-1}(a_i)+b_k\delta_k(a_i)\\&=0+0+\cdots+b_i+\cdots+0+0\\&=b_i\end{align}\\$

用拉格朗日多项式做类比，求解同余方程组（序言中的定义）

$x\equiv a_1 \pmod{n_1}\\ x\equiv a_2 \pmod{n_2}\\ \vdots\\ x\equiv a_k\pmod{n_k}$

我们的特解可以写成

$x_0=\sum_{i=1}^k a_i M_iN_i \\$

因为 $M_iN_i=1-m_in_i$ 我们能知道 $M_iN_i\equiv\left\{ \begin{matrix}1 \quad \quad i=j \\ 0 \quad\quad i\ne j\end{matrix} \right. \pmod{n_j}$ ，代入

$\begin{align}x_0&\equiv a_1M_1N_1+a_2M_2N_2+\cdots+a_iM_iN_i+\cdots+a_{k-1}M_{k-1}N_{k-1}+a_kM_kN_k \pmod{n_i}\\&\equiv0+0+\cdots+a_i+\cdots+0+0\\&\equiv a_i\end{align}\\$

回到一开始的例子，求解同余方程组

$x\equiv 2\pmod{3}\\ x\equiv 3\pmod{5}\\ x\equiv 2\pmod{7}$

首先，注意到

$\begin{alignat}{3} &5\times 7\times\color{red}{ (-1)}&&+3\times\color{red}{ 12}&&=1\\ &3\times7\times\color{red}{ 1}&&+5\times\color{red}{(-4)}&&=1\\ &3\times 5\times\color{red}{ 1}&&+7\times\color{red}{(-2)}&&=1 \end{alignat}\\$