欧几里得算法及其扩展

欧几里得算法及其扩展

前言

整除：对于整数 \(a(a\ne 0)\) 和 \(b\)，如果 \(\exists q\in Z\)，使得 \(b=a\times q\)，则称 \(a\) 能整除 \(b\)，记作 \(a\mid b\)。否则，称 \(a\) 不能整除 \(b\)，记作 \(a\nmid b\)。

最大公约数

公约数：如果一个整数是集合中任意一个元素的约数，则这个整数是该集合中的最大公约数。

最大公约数：集合的公约数中最大的一个，通常用 \(\gcd\) 表示。

最大公约数经典的算法有：欧几里得算法（辗转相除法）、更相减损术。

欧几里得算法

“欧几里得算法”又名“辗转相除法”

过程

已知两个正整数 \(a\) 和 \(b\)，其中 \(a>b\)

若 \(b\mid a\)，则 \(b\) 就是二者的最大公约数

若 \(b\nmid a\)，\(\gcd(a,b)=\gcd(b,a \bmod b)\)

“辗转相除法”与“更相减损术”都基于事实：

当 \(a \le kb\) 时，\(\gcd(a,b) = \gcd(a-kb,b)=\gcd(b,a-kb)\)，其中 \(k\) 为正整数。
当 \(k\) 取最大值 \(\lfloor \frac{a}{b} \rfloor\)，又可以写成 \(\gcd(a,b) = \gcd(a \bmod b,b) = \gcd(b, a \bmod b)\)

证明：

设 \(a\) 与 \(b\) 的最大公约数是 \(g\)，则 \(\exists x,y\in Z\) 使得 \(a=xg,b=yg\)

由于 \(g\) 是最大公约数，所以 \(x\) 与 \(y\) 中不存在大于 \(1\) 的公约数，即 \(x\) 与 \(y\) 互质

现在要证： \(g\) 也是 \(a-kb\) 与 \(b\) 的最大公约数
\(a-kb=xg-kyg=(x-ky)g\)
即证 \(x-ky\) 与 \(y\) 互质

假设 \(x-ky\) 与 \(y\) 不互质，它们的最大公约数为 \(p\)，则 \(x-ky=up, y=u'p\)

\(x=up+yk=up+u'pk=(u+u'k)p\)

此时 \(x\) 与 \(y\) 不互质，与已证条件不符，因此 \(x-ky\) 与 \(y\) 互质

证毕

实现

/**
 * 通过递归实现的欧几里得算法求解最大公约数
 * @param a 一个正整数
 * @param b 一个正整数
 * @return a与b的最大公约数
 */
int gcd(int a, int b) {
    int rest = a % b;
    if (rest == 0) {
        return b;
    }
    return gcd(b, rest);
}

可以写成迭代的形式

/**
 * 通过迭代实现的欧几里得算法求解最大公约数
 *
 * @param a 一个正整数
 * @param b 一个正整数
 * @return a与b的最大公约数
 */
int gcd(int a, int b) {
    if (a < b) {
        return gcd(b, a);
    }
    while (true) {
        int rest = a % b;
        if (rest == 0) {
            return b;
        }
        a = b;
        b = rest;
    }
}

最小公倍数

公倍数：如果一个整数能被集合中任意一个数整除，则这个整数是该集合的公倍数。

最小公倍数：集合的公倍数中最小的一个，通常用 \(lcm\) 表示。

对于正整数 \(a\) 和 \(b\)，由算术基本定理得：\(a=p_1^{x_1}p_2^{x_2}\cdots p_n^{x_n}\)，\(b=p_1^{y_1}p_2^{y_2}\cdots p_n^{y_n}\)

则由最大公约数与最小公倍数的定义得：

\[\begin{aligned} \gcd(a,b) &= p_{1}^{min(x_1,y_1)} p_{2}^{min(x_2,y_2)} \cdots p_{n}^{min(x_n,y_n)}\\ lcm(a,b) &= p_{1}^{max(x_1,y_1)} p_{2}^{max(x_2,y_2)} \cdots p_{n}^{max(x_n,y_n)} \end{aligned} \]

由于 \(x_k+y_k=min(x_k,y_k)+max(x_k,y_k)\)

所以 \(\gcd(a,b)\times lcm(a,b)=a\times b\)

所以 \(lcm(a,b)=\dfrac{a\times b}{\gcd(a,b)}\)

扩展欧几里得算法

常用于求 \(ax+by=\gcd(a,b)\) 的一组可行解

过程

\[\begin{aligned} ax+by&=\gcd(a,b)\\ &=\gcd(b,a \bmod b)&(1)\\ &=bx^{\prime}+(a \bmod b)y^{\prime}&(2)\\ &=bx^{\prime}+(a-b\cdot \big\lfloor\dfrac{a}{b}\big\rfloor)y^{\prime}&(3)\\ &=ay^{\prime}+b(x^{\prime}-\big\lfloor\dfrac{a}{b}\big\rfloor\cdot y^{\prime})&(4) \end{aligned} \]

1. 欧几里得定理
2. \(\gcd(a,b)=ax+by\)，\(b\)代入\(a\)，\(a \bmod b\)带入\(b\)
3. 正整数取模定义，\(a \bmod b= a-\big\lfloor\dfrac{a}{b}\big\rfloor\cdot b\)
4. 合并同类项

因此，\(ax+by=ay^{\prime}+b(x^{\prime}-\big\lfloor\dfrac{a}{b}\big\rfloor\cdot y^{\prime})\)

对比系数得：

\[\left\{\begin{array}{l} x=y^{\prime} \\ y=x^{\prime}-\left\lfloor\frac{a}{b}\right\rfloor y^{\prime} \end{array}\right. \]

当 \(b=0\) 时， \(a=\gcd(a,b)\)，由 \(ax+by=\gcd(a,b)\) 得 \(gcd\cdot x+0\cdot y=\gcd\)，则 \(x=1,y\in Z\) 是此时 \((a,b)\)的 一个整数解

作为人类，通常会选择令最后一层的 \(y=0\)

因为 \(y\) 在回代时会增长很快，无论正负

再从下至上递归求出 最开始的 \(ax+by=\gcd(a,b)\) 的一组整数解 \((x_0,y_0)\)

此时可以得到方程的通解为 \(x=x_0+\dfrac{b}{\gcd(a,b)}\cdot t\)，\(y=y_0-\dfrac{a}{\gcd(a,b)}\cdot t\)，其中 \(t\in Z\)

若要求得 \(x\) 的最小正整数解，可以通过 \((x_0\%b+b)\%b\) 求得

值域分析

由于 \(ax+by=\gcd(a,b)\) 的解有无数个，那求出的解如果趋近于无穷呢？

若 \(b\ne 0\)，扩展欧几里得算法求出的可行解必有 \(|x|\le b,|y|\le a\)

详见：值域分析 - OI Wiki

实现

int x = 0, y = 0;
int exgcd(int a, int b) {
    if (b == 0) {
        x = 1;
        y = 0;
        return a;
    }
    int gcd = exgcd(b, a % b), t = x;
    x = y;
    y = t - a / b * y;
    return gcd;
}

迭代方法

欧几里得算法可以通过迭代的方式求得最大公约数，那扩展欧几里得是否也能迭代呢？

答案是肯定的。

若 \(x,y\) 对应当前层，\(x^{\prime},y^{\prime}\) 对应下一层

初始情况，当 \(x=1\)，\(y=0\)，\(x^{\prime}=0\)，\(y^{\prime}=1\)时，则 \(ax+by=a\)，\(ax^{\prime}+by^{\prime}=b\)，方程成立。

由 \(\gcd(a,b)=\gcd(b,a-\left\lfloor\frac{a}{b}\right\rfloor b)\)，存储两层，顺次向下一层递推

\[\begin{aligned} ax+by&=a&(1)\\ a x^{\prime}+b y^{\prime}&=b&(2) \end{aligned} \]

\((1)-\left\lfloor\frac{a}{b}\right\rfloor\times(2)\) 得 \((3)\)：\(a\left(x-\left\lfloor\frac{a}{b}\right\rfloor x^{\prime}\right)+b\left(y-\left\lfloor\frac{a}{b}\right\rfloor y^{\prime}\right) =a-\left\lfloor\frac{a}{b}\right\rfloor b\)

接下来 \((2)\) 变 \((1)\)，\((3)\) 变 \((2)\)，继续递推，直至求出 \(gcd\)，即 \(b=0\)

int x = 0, y = 0;
int exgcd(int a, int b) {
    x = 1;
    y = 0;
    int nx = 0, ny = 1, q, t;
    while (b != 0) {
        q = a / b;
        
        x -= q * nx;
        t = x; x = nx; nx = t;
        
        y -= q * ny;
        t = y; y = ny; ny = t;
        
        a -= q * b;
        t = b; b = a; a = t;
    }
    // 结束后
    // a 是 gcd(a,b)
    // 如果 gcd(a,b)!=1 则 无整数解
    // 否则 x,y 是 ax+by=gcd(a,b) 的一组整数解
    return a;
}

int x, y;
int exgcd(int a, int b) {
    x = 1, y = 0;
    int nx = 0, ny = 1, q;
    while (b != 0) {
        q = a / b;
        std::swap(x -= q * nx, nx);
        std::swap(y -= q * ny, ny);
        std::swap(a -= q * b, b);
    }
    // 结束后
    // a 是 gcd(a,b)
    // 如果 gcd(a,b)!=1 则 无整数解
    // 否则 x,y 是 ax+by=gcd(a,b) 的一组整数解
    return a;
}

二元一次不定方程

关于 \(x,y\) 的形如 \(ax+by=c\) 的方程是二元一次不定方程

裴蜀定理

对于 \(\forall a,b\in Z,\ \exists x,y\in Z\) 使得 \(ax+by=\gcd(a,b)\)

证明（By other）：

首先，显然有 \(\operatorname{gcd}(a, b)|a\) ， \(\operatorname{gcd}(a, b)| b\) 。由整除的性质可知 $\forall x, y \in Z, g c d(a, b) \mid(a x+b y) $。设 \(s\) 为 \(a x+b y\) 的最小正值，故有 \(\operatorname{gcd}(a, b) \mid s\)。
设 \(k=\left\lfloor\frac{a}{s}\right\rfloor, r=a \bmod s\) ，则有 \(r=a-k(a x+b y)=a(1-k x)+b(-k y)\) ，发现 \(r\) 也为 \((a, b)\) 的线性组合，而 \(r\) 又在模 \(s\) 剩余系中，故 \(r \in[0, s-1]\)，又由于 \(s\) 为最小正值，故 \(r=0\)。
所以 \(s \mid a\)，同理 \(s \mid b\)，根据两个数的公约数必为这两个数的最大公约数的约数，有 \(s \mid \operatorname{gcd}(a, b)\)，而上文已证明 \(\operatorname{gcd}(a, b) \mid s\)，故 \(s=\operatorname{gcd}(a, b)\)。
进而得证，故 \(a x+b y=\gcd(a, b)\) 这个方程一定有解

判断方程有整数解的条件

方程 \(ax+by=c\) 有整数解的充要条件是 \(\gcd(a,b)\mid c\)

证明（By other）：

考虑将 \(ax+by=\operatorname{gcd}(a, b)\) 变为 \(ax+by=c\)
为了区分这两个方程，我们将前者改为 \(a^{\prime} x+b^{\prime} y=\operatorname{gcd}\left(a^{\prime}, b^{\prime}\right)\)
令 \(k=\operatorname{gcd}\left(a^{\prime}, b^{\prime}\right)\)，方程两边同除 \(k\)，得: \(\frac{a^{\prime}}{k} x+\frac{b^{\prime}}{k} y=1\)

方程两边同乘c得：\(\frac{a^{\prime} c}{k} x+\frac{b^{\prime} c}{k} y=c\)

只要通过换元，即令 \(a=\frac{a^{\prime} c}{k}, b=\frac{b^{\prime} c}{k}\)，即可将方程转化为 \(a x+b y=c\) 的形式
考虑到 \(\operatorname{gcd}\left(\frac{a^{\prime} c}{k}, \frac{b^{\prime} c}{k}\right)=c=\operatorname{gcd}(a, b)\)，且对于一般的方程: \(a t x+b t y=c^{\prime}\left(c^{\prime}=c t, t \in Z\right)\)
故只有 \(c\) 为 \(\operatorname{gcd}(a, b)\) 倍数时原方程有解，即 \(\operatorname{gcd}(a, b) \mid c\)，从而得证。

求解二元一次不定方程

对于解方程 \(ax+by=c\)，可以先解方程 \(ax+b y=\operatorname{gcd}\left(a, b\right)\)

当我们得到了 \(a x+b y=\operatorname{gcd}\left(a, b\right)\) 的一组解 \(x_0^{\prime},y_0^{\prime}\)，即 \(ax_0^{\prime}+by_0^{\prime}=gcd(a,b)\)

两边同时除以 \(gcd(a,b)\)，再乘以 \(c\)，即 \(a\dfrac{c\times x_0^{\prime}}{gcd(a,b)}+b\dfrac{c\times y_0^{\prime}}{gcd(a,b)}=c\)

则 \(x_0=\dfrac{c\times x_0^{\prime}}{\gcd(a,b)},\ y_0=\dfrac{c\times y_0^{\prime}}{\gcd(a,b)}\) 是 \(ax+by=c\) 的一组特解

此时可以得到方程的通解为 \(x=x_0+\dfrac{b}{\gcd(a,b)}\cdot t\)，\(y=y_0-\dfrac{a}{\gcd(a,b)}\cdot t\)，其中 \(t\in Z\)

若要求得 \(x\) 的最小正整数解，可以通过 \((x_0\%b+b)\%b\) 求得

线性同余方程

关于 \(x\) 的形如 \(ax\equiv c(mod\ b)\) 的方程被称为线性同余方程。

与 二元一次不定方程 等价

对于方程 \(ax\equiv c(mod\ b)\)，等价于 \(ax+by=c\)，其中 \(y=\left\lfloor\frac{c}{b}\right\rfloor-\left\lfloor\frac{a x}{b}\right\rfloor\)

证明：

\[\begin{aligned} ax&\equiv c(mod\ b)\\ ax\ mod\ c&= b\ mod\ c&(1)\\ ax-b\times \left\lfloor\frac{a x}{b}\right\rfloor&=c-b\times\left\lfloor\frac{c}{b}\right\rfloor&(2)\\ ax+(\left\lfloor\frac{c}{b}\right\rfloor-\left\lfloor\frac{a x}{b}\right\rfloor)b&=c&(3)\\ \end{aligned} \]

1. 同余号转等号
2. 取模定义 \(a\ mod\ b=a-b\times \big\lfloor\dfrac{a}{b}\big\rfloor\)
3. 合并同类项

求解线性同余方程

转化为 二元一次不定方程，求解该不定方程，即是求解线性同余方程

参考资料

最大公约数 - OI Wiki

求解一类方程的方法 - Nickel_Angel's nest