非递归的扩展欧几里得算法

$DeclareMathOperator{extgcd}{extgcd}$

作者按：写这篇随笔是为了解释 tourist 的逆元模板

template <typename T>
T inverse(T a, T m) {
    T u = 0, v = 1;
    while (a != 0) {
        T t = m / a;
        m -= t * a; swap(a, m);
        u -= t * v; swap(u, v);
    }
    assert(m == 1);
    return u; // 注意：u 可能为负数
}

扩展欧几里德算法是求不定方程 $ax + by = gcd(a,b)$ 的一组解 $(x, y)$；若 $a$ 和 $b$ 互素，$x$ 即 $a$ 在模 $b$ 下的逆元。注意到方程 $ax + by = gcd(a, b)$ 和方程 $frac{a}{gcd(a, b)} x + frac{b}{gcd(a, b)} y = 1$ 同解。以下只考虑 $a, b$ 互素的情形。常见的求逆元算法实现如下

int extgcd(int a, int b, int& x, int& y) {
  int d = a;
  if (b != 0) {
    d = extgcd(b, a % b, y, x);
    y -= (a / b) * x;
  } else {
    x = 1; y = 0;
  }
  return d;
}

int mod_inverse(int a, int m) {
    int x, int y;
    extgcd(a, m, x, y);
    return (m + x % m) % m;
}

求 $a$ 在模 $m$ 下的逆元时，对于方程 $ax + my = 1$，我们并不需要把 $x, y$ 都求出来，而只需要求出 $x$。相比于 tourist 的实现，上述代码是不够高效的，一则含有冗余计算，二则采用递归实现。下面分析 tourist 的非递归实现之原理。

一般地，可以通过下述过程求出 $gcd(a, m)$。
令 $r_0 = m$，$r_1 = a$。
$r_0 = r_1 q_2 + r_ 2 $
$r_1 = r_2 q_3 + r_3$
$r_2 = r_3 q_4 + r_4$
$vdots$
$color{blue}{r_{i - 2} = r_{i - 1} q_i + r_i}$
$r_{i - 1} = r_{i} q_{i + 1} + r_{i + 1}$
$r_{i} = r_{i + 1} q_{i + 2} + r_{i + 2}$
$vdots$
$r_{n - 2} = r_{n - 1} q_n + r_{n}$
$r_{n} = 0$

我们有

$gcd(r_0, r_1) = gcd(r_1, r_2) = dots = gcd(r_{n - 1}, r_{n}) = gcd(r_{n-1}, 0) = r_{n - 1} = 1$

设 $r_i = s_i a + t_i m$，注意到 $r_i = r_{i - 2} - r_{i - 1} q_i $，故有
egin{aligned}
r_i &= s_{i - 2} a + t_{i - 2} m - (s_{i - 1} a + t_{i - 1} m) q_i \
&= (s_{i - 2} - q_i s_{i - 1}) a + (t_{i - 2} - q_i t_{i - 1}) m
end{aligned}
即有递推
egin{aligned}
s_i &= s_{i - 2} - q_i s_{i - 1}, \
t_i &= t_{i - 2} - q_i t_{i - 1}.
end{aligned}

我们所感兴趣的是等式 $r_{n - 1} = s_{n - 1} a + t_{n - 1} m = 1$，易见 $s_{n-1}$ 即我们欲求的 $a$ 在模 $m$ 下的逆元。注意到 $s_i$ 的递推式中未出现 $t$，因此可以不管 $t$，只考虑 $s$。

下面考虑边界条件。
注意到 $r_0 = s_0 a + t_0 m = m$，$r_1 = s_1 a + t_1 m = a$，迳取 $s_0 = 0, s_1 = 1$。
循环开始之前，变量 u 存放 $s_0$，变量 v 存放 $s_1$。

对于 $i ge 0$，第 $i$ 轮循环：

1. 根据等式 $r_{i} = r_{i + 1} q_{i+2} + r_{i + 2} $，算出 $q_{i+2}$（即变量 t）和 $r_{i+2}$；
2. 根据递推式 $s_{i+2} = s_i - q_{i+2} s_{i + 1} $，算出 $s_{i+2}$。

第 $i$ 轮循环结束之后

• 最新的 $r$，即 $r_{i + 2}$，保存在变量 a 中，$r_{i + 1}$ 存放在 m 中；
• $s_{i + 1}$ 存放在 u 中，$s_{i + 2}$ 存放在 v 中。

通过上述分析可以知道当循环结束时，即 a 的值是 $r_n = 0$ 时，m 的值恰是 $r_{n - 1} = gcd(a, m)$。

题外话

我们看到两种实现所依据的都是欧几里得算法。作为补充，下面来讨论扩展欧几里得算法算出的 $ax + by = gcd(a,b)$ 的解的大小。

用归纳法可以证明，若 $ab e 0$，则 $|x| le b$ 且 $|y|le a$。

在 $ b = 0$ 的前一步，即 $a \% b = 0$ 时有 $x = 0$ 且 $y = 1$，结论显然成立。假设调用 extgcd(b, a % b, y', x') 后有 $|x'| le b$ 且 $|y'| le a \% b$ 。在 extgcd(a, b, x, y) 中 $x = x', y = y' - (a / b) x '$，所以有如下不等式成立
$|x| = |x'| le b$，
$|y| = |y' - (a / b)x'| le |y'| + (a /b )|x'| le a \% b + (a / b) imes b = a$

References

https://brilliant.org/wiki/extended-euclidean-algorithm/#extended-euclidean-algorithm