谈谈扩展欧几里得算法

❖扩展欧几里得算法

欧几里得算法是数论中用于求得两个整数 $a$ 和 $b$ 的最大公约数 $c=(a,b)$ 的算法。我之前的文章中已经介绍了这个算法的正确性来源于 $(a,b)=(a-k b), k\in \Z$ 。而“扩展欧几里得算法”是一个可以得到

中的系数 $s$ 和 $t$ 的算法。

我第一次看到的这个算法实现，来自于一位精通算法竞赛的同学，他给出的实现类似于

1int exgcd(int a,int b,int &x,int &y)
2{
3    if(b==0)
4    {
5        x = 1; y = 0; return a;
6    }
7    int r = exgcd(b, a % b, x, y);
8    int t = y;
9    y = x - (a / b) * y;
10    x = t;
11    return r;
12}

这个实现把参数x和y当作传递结果的寄存器，有比较浓厚的命令式风格。

用比较易于理解的函数式风格实现如下：

1let rec ext_gcd a b = 
2  if b = 0 then (1, 0, a)
3  else 
4    let quo = a / b
5    let r = a % b
6    let (s, t, gcd) = ext_gcd b r
7    (t, s - quo * t, gcd)

实际上，上面的两种实现都是基于一个手工的代换过程。以下述欧几里得算法过程为例：

将每个除法式都变成无除法的加减法式：

考虑倒数第二个式子

把这个式子里的 $\textcolor{red}{6}$ 用 $20 - 2 \times 7 = \textcolor{red}{6}$ 代换，我们得到：

更一般地说，两个式子为

那么，把他们进行合并代换，就会变为：

显然地，可以把这个过程抽象为一个对三元组的运算，这里用$*$表示，即：

这时再看那个函数式程序，我们会清晰很多：

ext_gcd b r求出的是$s b + tr = gcd$的三元组$(s, t, gcd)$，而此三元组与本次除法定义的$1 \times a + (-quo) \times b = r$之三元组$(1, -quo, r)$正好可以进行$*$运算，而

这正是代码中的返回值(t, s - quo * t, gcd).

我们可以给出直接使用这个运算的程序：

1let F x1 x2 = 
2  let s1, t1, c1 = x1
3  let s2, t2, c2 = x2 
4  (s1 * t2, s2 + t1 * t2, c2)
5
6let rec ext_gcd2 a b =
7  if b = 0 then (1, 0, a)
8  else 
9    let quo = a / b
10    let r = a % b
11    F (1, -quo, r) (ext_gcd2 b r)

❖从递归到循环

❖Wand 的模式

在40多年以前，MITCHELL WAND发表了一篇论文，叫做 Continuation-Based Program Transformation Strategies，他总结了把递归程序变成循环程序的很多通用方法。

当然，我们这里的“循环”和“递归”的真正意义还需要进一步规范，比如说，

1(define (factor n)
2  (let loop ([n n] [cont (lambda (x) x)])
3    (if (= n 0)
4        (cont 1)
5        (loop (- n 1) (lambda (x) (cont (* n x)))))))

这种“循环”不是此处研究的“循环”，我们所说的“循环”必须是“使用固定空间，空间复杂度为$O(1)$”的东西。

通过研究类似于factor这种函数的两种写法：

1(define (factor n)
2  (if (= n 0)
3      1
4      (* n (factor (- n 1)))))

1(define (factor n)
2  (let loop ([n n][acc 1])
3    (if (= n 0)
4        acc
5        (loop (- n 1) (* acc now)))))

WAND总结出了第一个模式：

要求是 $b$ 这个函数是 $A \rightarrow A \rightarrow A$ 的运算，且具有结合性，其中 $l_b$ 是右单位元（意思是说$b(x, l_b)=x$）。

我们是否可以用这个模式对上面的程序进行变换呢？看似是可以的，因为这里的 $p(x)$，$a(x)$，$b$，$F$，$c(x)$，$d(x)$ 都可以明确地给出。

然而，这里有一个问题：上面定义的运算 $*$ 到底是不是结合性的呢？

❖结合性？

一个运算的“结合性”，本质上描述了这个运算在“自组合”的时候能不能改变其运算顺序。考虑一个运算$f$，它具有结合性就是说

然而，我们的运算似乎没有结合性：

这是怎么回事？找个具体例子研究一下：

以前三个为例，如果正着做：

这会得到令人非常吃惊的结果，因为

而我们的组合运算理论上应该得到

问题出在哪里呢？

考虑 $(1, -3, 7) * (1, -2, 2)=(-2, 7, 2)$，等号右侧三元组的意义是：

这个式子不能和

用 $\*$ 进行运算。因为能与 $-2 \times 55 + 7 \times 16 = 2$ 进行 $*$ 运算的式子一定形如

可见，这已经不仅是没有结合性的问题。对于我们定义的 $\*$ 运算来说 $(a\*b)\*c$ 根本没有意义。

这样一来，单纯地用Wang的方法不可能把这个递归变换为循环。

❖如何循环实现？

再次观察上面的式子

我们好像发现了一种新的模式：

如果能够建立上面的关系，那么$s_n$和$t_n$就都是可以计算的：

令人头痛的是开头两个式子的建立。实际上读者如果认真研究过递归版的算法，最后一个式子的建立用了一个有趣的技巧。我们这里也可以运用相同的技巧，但这技巧为何成立就留给读者自行研究吧：

实现如下：

1let rec ext_gcd_loop a b =
2  if b = 0 then (1, 0, a)
3  else 
4    let rec loop a b s_n1 t_n1 s_n2 t_n2 = 
5      let quo = a / b 
6      let r = a % b
7      let s = 1
8      let t = -quo
9      if r = 0 then (s_n1, t_n1, b)
10      else 
11        loop b r (s * s_n2 + t * s_n1) (s * t_n2 + t * t_n1) s_n1 t_n1
12    loop a b 0 1 1 0