Extended Euclid Bug

简介

辗转相除法是经典的求最大公约数的算法，也是已知的最古老的算法，最近学到了扩展的辗转相除法，也就是对于$a$和$b$不仅返回$\gcd(a,b)$,而且返回$x$和$y$使得$xa + yb=\gcd(a,b)$.

在实现的过程中遇到了我觉得还算有意思的一点bug，就记录了下来. 用的Python.

算法

先来看一下算法，我们gcd(a, b)函数需要返回三个值g, x, y，g是a和b的最大公约数，而x和y满足xa+yb=g，依旧使用经典的递归实现。

Bug

先看两个有Bug的代码

import numpy as np
def gcd_pos(a, b):
    if b == 0:
        return a, 1, 0
    g, x, y = gcd(b, a % b)
    return g, y, x + y * (- a // b)

def gcd(a, b):
    g, x, y = gcd_pos(abs(a), abs(b))
    return g, x * np.sign(a), y * np.sign(b)

当时发现这个程序有Bug的时候，我一时没有回过神来，不过这其实是一个Python语法相关的Bug，因为gcd_pos保证了参数是正数，那么(- a // b)就是一个负数对一个正数求整数除法商，而Python中的整数除法是通过去掉小数部分完成的，并不是向下取整完成的。

所以这个看起来和递归公式无比契合的代码其实是错的.它表达的意思是$x + \lceil -a/b\rceil$.

import numpy as np
def gcd_pos(a, b):
    if b == 0:
        return a, 1, 0
    g, x, y = gcd(b, a % b)
    return g, y, x - y * a // b

def gcd(a, b):
    g, x, y = gcd_pos(abs(a), abs(b))
    return g, x * np.sign(a), y * np.sign(b)

这个Bug则更加难于发现，Python的表达式乘除计算是从左到右的，所以y * a // b实际上是(y * a) // b。

那么这个代码表达的是$x-\lfloor (ya)/b\rfloor$.根据数论知识，我明白了这是与递归公式不同的。

正确实现

import numpy as np
def gcd_pos(a, b):
    if b == 0:
        return a, 1, 0
    g, x, y = gcd(b, a % b)
    return g, y, x - y * (a // b)
def gcd(a, b):
    g, x, y = gcd_pos(abs(a), abs(b))
    return g, x * np.sign(a), y * np.sign(b)

总结

代码与数学公式之间的对应关系是很微妙的，以后编程时要注意这一点。