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

扩展欧几里得算法可以算出方程 \(ax + by = \gcd(a, b)\)（\(a, b \in \Z\)）的一组整数解。本文假设读者已经知道扩展欧几里得算法的递归实现，下面我们来推导它的非递归实现。

首先回顾欧几里得算法，常见的实现是递归的

// precondition: a，b 非负且不全是零。
int gcd(int a, int b) {
  return b == 0 ? a : gcd(b, a % b);
}

我们考虑上述 gcd 函数的执行过程，设初次调用 gcd 的两个参数是 \((x_1, x_2)\)。若 \(x_2 \ne 0\) 则第二次调用 gcd 的两个参数是 \((x_2, x_1 \bmod x_2)\)，令 \(x_3 = x_1 \bmod x_2\)。这样，在 gcd 执行的过程中就得到序列 \(x_1, x_2, \dots, x_{n-1} = \gcd(x_1, x_2), x_n = 0\)。第一次调用 gcd 的参数是 \((x_1, x_2)\)，第二次调用 gcd 的参数是 \((x_2, x_3)\)，第 \(n - 1\) 次调用 gcd 的参数是 \((x_{n-1}, x_n)\)。序列 \(x\) 满足 \(x_{i} = x_{i - 2} \bmod x_{i - 1}\)（\(i \ge 3\)）。

我们采用递推法。设 \(x_i = p_ix_1 + q_i x_2\)。我们有 \(x_{i} = x_{i - 2} \bmod x_{i - 1}\) 即 \(x_{i} = x_{i-2} - \floor{x_{i-2}/x_{i-1}} x_{i-1}\)，从而有 \[\begin{equation} (p_i, q_i) = (p_{i-2}, q_{i-2}) - \floor{x_{i-2}/x_{i-1}}(p_{i-1}, q_{i-1})。 \end{equation}\] 观察上式，不难看出 \(p, q\) 这两个序列是独立的，我们只考虑序列 \(p\)。我们的目标是求出 \(p_{n - 1}\)。取 \(p_1 = 1, p_2 = 0\)。序列 \(p\) 满足递推 \[\begin{equation} p_i = p_{i - 2} - \floor{x_{i-2}/x_{i-1}} p_{i-1}。 \end{equation}\]

std::pair<int, int> ext_gcd(int x1, int x2) {
  int p1 = 1, p2 = 0;
  while (x2 != 0) {
    int t = x1 / x2;
    p1 -= t * p2;
    x1 -= t * x2;
    swap(p1, p2);
    swap(x1, x2);
  }
  return {p1, x1};
}