【算法笔记】扩展 Euclid 算法
- 本文总计约 8000 字，阅读大约需要 30 分钟。
- 警告！警告！警告！本文有大量的 \(\LaTeX\) 公式渲染，可能会导致加载异常缓慢！
前言
扩展 Euclid 算法是数论中最经典，也是最简单（？）的算法了，但是它的应用又十分的广泛。所以介绍它是十分有必要的。
数论的知识真的是非常奇怪呢（bushi），因为它们感性地理解起来非常轻松，但是要严谨证明就很难 QAQ，笔者在文中加入了大量的数学证明，尽管可能有读者不爱看因为其实我也不爱看，但是基于数学严谨的思想，我还是把它们加了进来，这可能会导致浏览器渲染起来蛮麻烦的，不过还是请读者多多包涵了。
题目引入
给定三个整数 \(a,b,c\)，求关于 \(x,y\) 的二元一次不定方程 \(ax+by=c\) 的一组整数解，注意该方程可能没有整数解。
样例输入 \(1\)：\(a=3,b=2,c=7\)；
样例输出 \(1\)：\(x=1,y=2\)。注意可能有其它的满足要求的解，如 \(x=3,y=-1\)。
样例输入 \(2\)：\(a=2,b=-4,c=1\)；
样例输出 \(2\)：无解，因为 \(2x-4y\) 永远是偶数，但 \(c=1\)，是奇数。
最大公约数及 Euclid 算法
最大公约数的概念及性质
尽管最大公约数在我们小学的时候就已经接触了，但我们还是需要重新介绍一下它：
如果对于两个整数 \(a,b\)，若 \(m\) 是最大的既可以整除 \(a\)（记作 \(m\mid a\)），又可以整除 \(b\) 的正整数，即 \(m\) 是 \(a,b\) 最大的公约数，那么我们称 \(m\) 为 \(a,b\) 的最大公约数（好像很废话的概念 QwQ），记作 \(m=\gcd(a,b)\)。
特别地，若两个整数 \(a,b\) 满足 \(\gcd(a,b)=1\)，则我们称 \(a\) 与 \(b\) 互质，记作 \(a\perp b\)。
显然 \(\gcd\) 函数有如下的性质，此处不加以证明：
- 交换律：\(\gcd(a,b)=\gcd(b,a)\)，特别地，\(\gcd(a,0)=\gcd(0,a)=a\)；
- 可以提出常数：\(\gcd(ka,kb)=k\gcd(a,b)\ (k\in\mathbb{Z}\))；
- 对加减法及取模自由：\(\gcd(a,b)=\gcd(a,b\pm a)=\gcd(a,b\pm ka)=\gcd(a,b \bmod a)\ (k\in\mathbb{Z})\)。
Euclid 算法的内容
现在我们已经知道了最大公约数的基本概念了，那么我们怎么求任意两个给定的整数 \(m,n(m\ge n)\) 的最大公约数呢？
答曰：我会 \(\mathcal{O}(n)\) 暴力枚举！
然而，假如我们令 \(n\le 10^9\) 甚至是 \(n\le 10^{18}\)，那么 \(\mathcal{O}(n)\) 的枚举是一定会超时的 QwQ，我们就需要一种更快的算法来求出上述问题中的 \(\gcd(m,n)\)。于是，Euclid 算法应运而生。
Euclid 算法，又称辗转相除法，顾名思义，是由古希腊数学家 \(\mathcal{Euclid}\) 发明的算法，其步骤如下：
- 对于两个整数 \(m,n\)，如果有 \(n\ne0\)，那么我们求出 \(r=m\bmod n\) 的值；
- 接下来，使 \(m=n,n=r\)；
- 重复上述两个步骤，直到 \(n=0\) 时，\(m\) 的值就是我们所求的最大公约数。
例如，我们要求 \(\gcd(1081,897)\)，步骤如下：
- \(n=897\ne0\)，故我们令 \(m=897\)，\(n=1081\bmod897=184\)；原问题即转化为求 \(\gcd(897,184)\)；
- \(n=184\ne0\)，故我们令 \(m=184\)，\(n=897\bmod184=161\)；原问题即转化为求 \(\gcd(184,161)\)；
- \(n=161\ne0\)，故我们令 \(m=161\)，\(n=184\bmod161=23\)；原问题即转化为求 \(\gcd(161,23)\)；
- \(n=23\ne0\)，故我们令 \(m=23\)，\(n=161\bmod23=0\)；此时 \(n=0\)，所以 \(\gcd(1081,897)=m=23\)。
像这样，不断地用一个数去除以另一个数，辗转不停，故得名为辗转相除法。
代码实现
Euclid 算法的代码非常好写，这里给出迭代求最大公约数的代码：
#include <cstdio>
#include <algorithm> //有用于求 gcd 的标准函数 __gcd(int, int)
using namespace std;
inline int Euclid(int x, int y) { //辗转相除法求 gcd(x, y)，inline 防止无效信息入栈
int r;
while(y) {
r = x % y;
x = y;
y = r;
}
return x;
}
long long m, n;
int main(void) {
scanf("%lld%lld", &m, &n);
printf("%lld", Euclid(m, n));
return 0;
}
//by CaO
当然，C++ 中 <algorithm> 库里有 __gcd(int, int) 函数，功能与上述函数相同，但我不知道 NOIP 中让不让用 QwQ。
时间复杂度分析
（前方数学证明警告！）
有一种非常经典的证明其时间复杂度的方法，如下：
我们构造数列 \(\{U\}\)，其递推式为：
定义 Fibonacci 数列 \(\{F\}\) 递推式为：
显然数列 \(\{U\}\) 中存在一项 \(U_x=0\)。
故有 \(F_0=U_x=0,F_1\ge U_x\)。
又因为 \(U_i \bmod U_{i+1}=U_{i+2}\)，所以 \(U_{i}\ge U_{i+1}+U_{i+2}\)。
利用数学归纳法可得：\(U_i\ge F_{x-i+1}\)，故 \(n=U_1\ge F_{x}\)。
注意到 \(F_x>\dfrac{(\phi+1)^x}{\sqrt{5}}\)，其中 \(\phi=\dfrac{\sqrt{5}-1}{2}\)，所以 \(x<\log_{(\phi+1)}(\sqrt{5}n)\)。
故 Euclid 算法迭代的次数不超过 \(\log_{(\phi+1)}(\sqrt{5}n)\) 次，时间复杂度为 \(\mathcal{O}(\log n)\)。\(\Box\)
这就是说，Euclid 算法相比较于 \(\mathcal{O}(n)\) 的暴力枚举，有了很大程度的优化。
不定方程概念的引入
不定方程
一个含有多个未知数的整系数方程称为不定方程，又称 Diophantine 方程。即，形如：\(\sum_{i=1}^{n}a_ix_i^{p_i}=b\) 的方程即为不定方程，其中 \(a_i,p_i,b\) 均为整数（\(i=1,2,\cdots,n\)）。如果这个方程只有两个未知数，且未知数的次数均为 \(1\)，那么我们称其为二元一次不定方程。
注意，本文在讨论不定方程时，默认指的是二元一次不定方程，即形如 \(ax+by=c\) 的方程。
不定方程有解问题 — — Bézout 定理
正如题目引入中所说的，二元一次不定方程是不一定有解的，例如 \(2x-4y=1\) 就没有整数解。所以，我们应该如何判断，一个不定方程是否有解呢？
针对这个问题，早在 \(17\) 世纪，法国数学家 \(\mathcal{B\acute{e}zout}\) 就给出了结论：
关于 \(x,y\) 的不定方程 \(ax+by=c\) 有整数解，当且仅当 \(\gcd(a,b)\mid c\)。或者用另一种表述：\(ax+by=\gcd(a,b)\) 一定有整数解。这就是著名的 Bézout 定理，或者国内又译为裴蜀定理或贝祖定理。证明如下：（数学证明警告！）
先证充分性，即若 \(\gcd(a,b) \mid c\)，则 \(ax+by=c\) 有整数解。
设 \(S\) 是最小的能用 \(ax+by\) 表示的正整数，则有
\(a\bmod S=a-S\left \lfloor \dfrac{a}{S} \right \rfloor=a-(ax+by)\left \lfloor \dfrac{a}{S} \right \rfloor=a \left(1-\left \lfloor \dfrac{a}{S} \right \rfloor \right)-b\left(\left \lfloor \dfrac{a}{S} \right \rfloor\right)\le S\)。
又因为 \(S\) 是最小的能用 \(ax+by\) 表示的正整数，所以 \(a\bmod S=0\)，即 \(S\mid a\)；
同理 \(S\mid b\)。即 \(S\) 是 \(a,b\) 的公约数，由最大公约数定义得 \(S\le \gcd(a,b)\)。
又因为 \(\gcd(a,b)\mid (ax+by)\)，所以 \(\gcd(a,b)\mid S\)，故 \(\gcd(a,b) \le S\)。
综上 \(S=\gcd(a,b)\)，即 \(ax+by=\gcd(a,b)\) 有整数解。
且对于任意 \(c=kS\ (k\in\mathbb{Z})\)，都存在 \(x'=kx,y'=ky\)，使得 \(ax'+by'=c\)，即 \(ax+by=c\) 存在正整数解。
再证必要性，即若 \(ax+by=c\) 有整数解，则 \(\gcd(a,b) \mid c\)。
假设 \(x_0,y_0\) 是 \(ax+by=\gcd(a,b)\) 的整数解，则任意的 \(ax+by=c\) 的整数解 \(x,y\)，都有 \(x=\dfrac{c}{\gcd(a,b)}x_0,y=\dfrac{c}{\gcd(a,b)}y_0\)。
若 \(\gcd(a,b)\nmid c\)，则 \(x,y\notin\mathbb{Z}\)，与 \(x,y\) 为整数矛盾。
综上所述，不定方程 \(ax+by=c\) 有整数解当且仅当 \(\gcd(a,b)\mid c\)。\(\Box\)
Bézout 定理的推广
当然，即使 Bézout 定理是关于二元一次不定方程的，但我们依旧可以将其推广到 \(n\) 元的不定方程：
关于 \(x_1,x_2,\cdots,x_n\) 的 \(n\) 元一次不定方程 \(\sum_{i=1}^{n}a_ix_i=b\) 有整数解，当且仅当 \(\gcd_{i=1}^{n}(a_i)\mid b\)。证明留给读者作为练习。
不定方程的解
扩展 Euclid 算法的引入
对于不定方程 \(ax+by=c\)，我们可以通过 Euclid 算法求出 \(\gcd(a,b)\) 的值，然后用 Bézout 定理判断其是否存在整数解。
但是，我们只是知道了一个不定方程的可解性，又如何求它的具体解呢？
其实，这个问题通过对 Euclid 算法进行一些小小的拓展、变形，就可以求出来。也就是接下来的扩展 Euclid 算法。
扩展 Euclid 算法的内容
我们现在考虑最简单的一种情况：\(ax+by=\gcd(a,b)\)，应该如何求其整数解呢？
首先，小学的数学知识就告诉我们，商乘除数加余数等于被除数，这一点是永远成立的，所以我们用 \(b\) 除以 \(a\)，有 \(b=a\left \lfloor \dfrac{b}{a} \right \rfloor + (b\bmod a)\)。代入原方程，有 \(ax+\left(a\left \lfloor \dfrac{b}{a} \right \rfloor + (b\bmod a)\right)y=\gcd(a,b)\)。
重新合并，得到了 \(a\left(\left \lfloor \dfrac{b}{a} \right \rfloor y+x\right)+(b\bmod a)y=\gcd(a,b)\)。
右边由 \(\gcd\) 的性质有 \(\gcd(a,b)=\gcd(a,a\bmod b)\)，我们令左边 \(x'=y,y'=\left \lfloor \dfrac{b}{a} \right \rfloor y+x\)，则有 \((a\bmod b)x'+ay'=\gcd(a,a\bmod b)\)，此时，我们就得到了一个新的，系数范围更小的不定方程。
如果我们求出了方程中的 \(x',y'\)，就可以代入回去，得到 \(y=x',x=y'-\left \lfloor \dfrac{b}{a} \right \rfloor x'\)。整个算法就可以递归地进行了。
递归的出口为 \(b=0\)，此时我们有 \(x=1,y=0\)。
像这样，我们利用了系数与系数之间 “辗转相除” 的方式，求出了 \(ax+by=\gcd(a,b)\) 的一组解，同时，我们也可以求出 \(\gcd(a,b)\) 的值，即递归出口处 \(a\) 的值。
这种算法是通过 Euclid 算法变形而来，所以我们称其为扩展 Euclid 算法，简称扩欧。
拓展情形
通过扩欧算法，我们可以求出 \(ax+by=\gcd(a,b)\) 的一组解。但是，根据 Bézout 定理，对于任意的 \(\gcd(a,b)\mid c\)，
代码
#include <cstdio>
using namespace std;
int a, b, c, k, x, y;
int exgcd(int a, int b, int& x, int& y) { //扩欧求不定方程
if(b == 0) {
x = 1;
y = 0;
return a; //递归出口
}
int ans = (b, a % b, y, x);
y -= (a / b * x); //状态转移
return ans;
}
int main(void) {
scanf("%d%d%d", &a, &b, &c);
if(c % exgcd(a, b, x, y) != 0) { //根据裴蜀定理，此时方程无整数解
printf("No Solution.\n");
}
else {
k = c / exgcd(a, b, x, y);
printf("x=%d, y=%d\n", k * x, k * y);
}
return 0;
}
//by CaO
时间复杂度分析
扩欧的算法与 Euclid 算法步骤基本相同，时间复杂度依旧为 \(\mathcal{O}(\log n)\)。
扩展 Euclid 算法的应用
扩欧可以用于同余方程的求解：形如 \(x\equiv a\pmod{m}\) 的方程称为同余方程。我们通过余数的性质，可以将其改写为不定方程的形式：\(x=km+a\)。移项，得 \(x-km=a\)，其中 \(x,k\) 是代求的未知数，我们可以通过扩欧求出来 \(x\)，也就是可以求出同余方程的解了。这种情况下，我们就可以将若干个同余方程连接成同余方程组，并求出满足同余方程组的解— —这就是著名的“中国剩余定理”（CRT）。
同时，我们可以依靠乘法逆元，在取模意义下完成除法，所谓乘法逆元，指的就是使得给定的 \(a,m\) 满足 \(ax\equiv1\pmod{m}\) 的 \(x\) 的值。我们将其可以变形为 \(ax+km=1\) 的形式，于是我们就可以通过扩欧求出 \(x\)，也就是求出 \(a\) 模 \(m\) 意义下的乘法逆元。
例题
本题目列表会持续更新。
浙公网安备 33010602011771号