[算法][数论] 扩展欧几里得算法
算法 数论 ExGCD
Lastmod: 2020-11-01 周日 19:33:23

扩展欧几里得算法

有时我们需要求解类似于 $ax+by=m$ 的不定方程(丢番图方程),扩展欧几里得算法就是求解这样方程的利器。

除此以外我们可以通过简单地变形处理模逆元和线性同余方程。

裴蜀定理/贝祖定理(Bezout’s Lemma)

$\forall a,b \in \mathbb{Z}$ 且 $a,b$ 不同时为零,方程

$$ax + by = m$$

有整数解,当且仅当

$$(a,b) \mid m$$

证明:

当 $a,b$ 其一为 $0$ 时,假设 $a=0$ 则 $(a,b)=b$ 显然只有 $b \mid m$ 时 $0x + by=m$ 有解,且为无穷多解($x$ 为任意值)。

当 $a,b$ 全不为 $0$ 时,设

$$A = \{xa+yb; (x,y)\in \mathbb{Z}^2\}$$

我们知道 $A \cap \mathbb{N}^+$ 是非空的(其中至少有 $|a|$,$|b|$),由于 $\mathbb{N}^+$ 是良序的,因此我们设 $A \cap \mathbb{N}^+$ 的最小元为

$$d_0 = x_0 a + y_0 b$$

考虑 $A$ 中任意元素 $p=x_1 a + y_1 b$ 对 $d_0$ 的带余除法,则

$$p = qd_0 +r$$

其中 $q>0,0 \le r < d_0$。

此时

$$\begin{align} r =& p - qd_0 \\ =& x_1 a + y_1 b - q(x_0 a + y_0 b) \\ =& (x_1-x_0q )a + (y_1 - y_0q)b \quad \in A \end{align}$$

由 $d_0$ 是 $A$ 中最小正元,而 $0 \le r < d_0$ 且 $r \in A$ 则 $r=0$

也就是说

$$d_0 \mid p$$

由 $p$ 是 $A$ 中任意元素,则 $d_0 \mid a$,$d_0 \mid b$,从而得出

$$d_0 \mid (a,b)$$

又有 $(a,b)\mid a,(a,b) \mid b$,则

$$(a,b) \mid x_0a + y_0b = d_0$$

从而得出

$$d_0 = (a,b)$$

必要性:

当 $ax+by=m$ 有解,则由 $(a,b) \mid a,(a,b) \mid b$ 得出 $(a,b) \mid m=xa+yb$

充分性:

当 $(a,b) \mid m$,此时我们知道 $m=m_0 d_0$,

由 $d_0 = x_0a+y_0b$

则我们至少有解 $x=x_0m_0, y=y_0 m_0$。

如果还想进一步证明解是无数多的,可以看通解形式部分。

通解形式

考虑构造解。 $x$ 增加 $\Delta x$,

$$\begin{align} a(x_0m_0+\Delta x) + b(y_0m_0+ \Delta y) &= m \\ a \cdot \Delta x + b \cdot \Delta y &= 0 \\ \Delta y &= -\frac{a}{b} \Delta x \end{align}$$

由于最终 $x,y \in \mathbb{Z}$ ,而我们已知 $x_0m_0,y_0m_0 \in \mathbb{Z}$,也就是说我们需要 $\Delta x, \Delta y \in \mathbb{Z}$。

对于 $\frac{a}{b}$ 我们可以考虑把分子分母的最大公因数 $(a,b) = d_0$ 除掉得到 $\frac{a'}{b'}$。

此时要满足

$$\Delta x, \Delta y = -\frac{a'}{b'} \Delta x \in \mathbb{Z}$$

由 $(a',b') = 1$ 则

$$\Delta x = k b' = k \cdot \frac{b}{d_0}$$

于是我们得到其通解

$$x=x_0m_0+k \cdot \frac{b}{(a,b)}, y=y_0m_0-k \cdot \frac{a}{(a,b)} \quad (k \in Z)$$

都是方程的解,且彼此不同($a,b$ 均不为 $0$ 则显然 $x$ 彼此不同),即解有无穷组。

构造一组特解

上述过程中我们只是证明了解的存在性和解之间的关系,但其实有一个重点问题没有解决,我们如何找到一组 $x_0,y_0$,满足不定方程 $ax+by=(a,b)$。

扩展欧几里得算法算法首次接触难以理解的一个原因是因为这是一个构造算法,因此思维上存在一定的跳跃性。它巧妙地利用了欧几里得算法的过程推导出了一条特解。

首先我们设 $g=(a,b)$

(1)考虑方程 $g \cdot x+0 \cdot y=g$ 我们可以不假思索的给出一组特解

$$x_0=1,y_0=0$$

而对应的系数 $a_0,b_0$ 刚好是我们在辗转相除求 $(a,b)$ 过程中的最后一步的系数。其实 $y_0$ 可以选择任意整数,但是选择 $0$ 我们可以将特解限制在较小的范围,见下。

(2)我们发现对于辗转相除的过程,如果 $a'=b,b'=a%b$ 我们发现对于这两组系数,我们有相同的最大公约数 $g$,则

$$ax+by=g=a'x'+b'y'$$

我们递归知道 $x',y'$ 的值之后,可以反过来构造 $x,y$ 则

$$ \begin{align} ax+by &= a'x'+b'y' \\ ax+by &= bx'+(a\%b) y' \\ ax+by &= bx'+(a-\lfloor \frac{a}{b} \rfloor \cdot b) y' \\ ax+by &= bx'+ay'-\lfloor \frac{a}{b} \rfloor \cdot b y' \\ ax+by &= a \cdot y' + b (x' - \lfloor \frac{a}{b} \rfloor \cdot y') \end{align} $$

通过比较系数,我们可以得到一组 $x,y$

$$\begin{align} x &= y' \\ y &= x' - \lfloor \frac{a}{b} \rfloor \cdot y' \end{align} $$

因此我们只需要递归求解 $(a,b)$ 然后回溯时求解 $x,y$ 就可以得出原不定方程的一组解了。

特解的范围

这里我们可以进一步研究一下特解的范围,因为反推的过程如果构造的特解过大,就可能因为溢出而无法使用。

我们发现当 $a=0$ 或 $b=0$ 时递归只会进行不超过两层,这样最终解的值为 $x=0,y=1(a=0)$ 或 $x=1, y=0(b=0)$。

同样如果输入 $a=b$ 时递归只会进行两层,这样最终得出的值为 $x=0,y=1$。

如果输入的 $a,b$ 不同,则取模操作不会让两个参数在递归的过程中变得相同(总会是一个比另一个更小)这里我们假设输入的 $a>0,b>0,a\neq b$,来看一下对应的 $a,b,x,y$ 的值:

No. a b x y
0 $g$ $0$ $1$ $0$
1 $kg \quad (k > 1)$ $g$ $0$ $1$
2 $jk g+g \quad (k>1,j \ge 1)$ $kg$ $1$ $-j$
3 $i(jkg+g) + kg=(ij+1)kg+ig \quad (k>1,j \ge 1, i \ge 1)$ $jk g+g$ $-j$ $1+ij$

我们发现除了第 $0$ 行,总有 $|x|<b,|y|<a$ 成立。

这里可以尝试证明一下,假设 $a',b',x',y'$ 是上一层递归的参数和结果,则假设 $|x'| < b', |y'| < a'$ 成立。

我们有

$$ \begin{align} a' &= b \\ b' &= a\%b = a - \lfloor \frac{a}{b}\rfloor \cdot b \\ x &= y' \\ y &= x' - \lfloor \frac{a}{b} \rfloor \cdot y' \\ \end{align} $$

对于 $x$,有

$$|x| = |y'| < a'=b$$

对于 $y$, 有

$$|y| = |x' - \lfloor \frac{a}{b} \rfloor \cdot y'|$$

如果 $x'$ 与 $y'$ 同号则

$$ \begin{align} |y| &= |x' - \lfloor \frac{a}{b} \rfloor \cdot y'| \\ &\le \left| |x'| + \left|\lfloor \frac{a}{b} \rfloor \cdot y' \right| \right| \quad (考虑 x',y' 同号异号两种情况)\\ &\le |x'| + \left|\lfloor \frac{a}{b} \rfloor \cdot y' \right| \\ &\le |x'| + \lfloor \frac{a}{b} \rfloor \cdot |y'| \\ &< b' + \lfloor \frac{a}{b} \rfloor \cdot a' \\ &< (a - \lfloor \frac{a}{b}\rfloor \cdot b) + \lfloor \frac{a}{b} \rfloor \cdot (b) \\ &< a \end{align} $$

得证。

至此我们终于可以书写代码了。

代码

void exgcd(LL a,LL b,LL &g,LL &x,LL &y)
{
	if(!b)
    {
        g=a; x=1; y=0;
    }
	else
    {
        exgcd(b,a%b,g,y,x); // 巧妙地利用参数顺序,调整好了对应的 x,y
        y-=a/b*x;
    }
}

扩展欧几里得算法的应用

求解不定方程

显然这个算法可以用来求解形如 $ax+by=m$ 的不定方程嘛。

不过做题的过程一般情况是问解是否存在,寻找 $x,y$ 满足某些条件的特解,或者是问 $x,y$ 满足某些条件时的解的个数。

解是否存在的问题其实根据裴蜀定理和 GCD 就已经足够求解了。对于满足某些条件的特解,相当于我们要先将通解的形式写出来,然后求解通解符合条件的方案。

求解线性同余方程

对于线性同余方程 $ax \equiv b \pmod m$ 我们可以简单的变形为

$$ax+my=b$$

之后进行求解。我们也可以知道这个方程有解的条件是 $b|(a,m)$。以及根据我们推导的结论,求出的 $|x|$ 会在 $[0,m)$。

求解逆元

更进一步的,其实我们令上面同余方程中的 $b = 1$,于是根据模下乘法逆元的定义求出的 $x$ 即为模下逆元(注意可能是负的,需要调整到 $[0,m)$)。

LL inv(LL a,LL m)
{
    LL g,x,y;
    exgcd(a,m,g,x,y);
    if(g!=1) return -1; // no inverse element
    x%=m; if(x<0) x+=m;
    return x;
}
Prev: [算法][数论] 最大公约数(GCD)与欧几里得算法
Next: [模板][数论] 扩展欧几里得算法(ExGCD)