欧几里得算法与扩展欧几里得算法

求最大公约数，一般采用gcd算法。http://zh.wikipedia.org/wiki/%E6%AC%A7%E5%87%A0%E9%87%8C%E5%BE%97%E7%AE%97%E6%B3%95

gcd算法简单高效，是对数级别的算法。

下面给出它的递归形式和迭代形式。

def gcd(a, b):
    if b == 0:
        return a
    return gcd(b, a%b)

def gcd1(a, b):
    while b != 0:
        r = a % b
        a = b
        b = r
    return a

print gcd(132, 324)
print gcd1(132, 324)

这次主要是说gcd算法的一个扩展，egcd算法。http://zh.wikipedia.org/wiki/%E6%89%A9%E5%B1%95%E6%AC%A7%E5%87%A0%E9%87%8C%E5%BE%97%E7%AE%97%E6%B3%95

基础：给出任意a, b，必有a*x + b*y = gcd(a, b)。

因为gcd(a, b) = gcd(b, a mod b)，所以一个简单实现是利用gcd算法算出gcd(a, b)再倒回去算 x 和 y 。

下面给出的递归形式，原理如下：

gcd(a, b) = a*x_i + b*y_i
gcd(b, a % b) = b*x_i+1 + (a - [a/b]*b)*y_i+1    （q = [a/b]）
gcd(a, b) = gcd(b, a % b) => a*x_i + b*y_i = a*y_i+1 + b*(x_i+1 - [a/b]*y_i+1)
x_i = y_i+1
y_i = x_i+1 - [a/b]*y_i+1

def extendedGCD1(a, b):
    # a*xi + b*yi = ri
    if b == 0:
        return (1, 0, a)
    (x, y, r) = extendedGCD1(b, a%b)
    """
    gcd(a, b) = a*xi + b*yi
    gcd(b, a %  b) = b*xi+1 + (a - [a/b]*b)*yi+1
    gcd(a, b) = gcd(b, a %  b)   =>   a*xi + b*yi = a*yi+1 + b*(xi+1 - [a/b]*yi+1)
    xi = yi+1
    yi = xi+1 - [a/b]*yi+1
    """
    tmp = x
    x = y
    y = tmp - (a/b) * y
    return (x, y, r)

迭代形式基于下面一个推导：

r₁ = a mod b

r₂ = b mod r₁

r₃ = r₁ mod r₂

...

r_i = r_i-2 mod r_i-1   即（1） ri = r_i-2 - q*r_i-1   （q = [r_i-2 / r_i-1]）

（令r_-1 = a, r₀ = b）

因为 a*x + b*y = gcd(a, b)即有，a*x_i + b*y_i = r_i   （中间的每一步）

代入（1）式有，a*x_i + b*y_i = (a*x_i-2 + b*y_i-2 ) - q * (a*x_i-1 + b*y_i-1 ) = a * (x_i-2 - q*x_i-1) + b * (y_i-2 - q*y_i-1)

则有，x_i  = x_i-2 - q*x_i-1，y_i = y_i-2 - q*y_i-1

def extendedGCD(a, b):
    #a*xi + b*yi = ri
    if b == 0:
        return (1, 0, a)
    #a*x1 + b*y1 = a
    x1 = 1
    y1 = 0
    #a*x2 + b*y2 = b
    x2 = 0
    y2 = 1
    while b != 0:
        q = a / b
        #ri = r(i-2) % r(i-1)
        r = a % b
        a = b
        b = r
        #xi = x(i-2) - q*x(i-1)
        x = x1 - q*x2
        x1 = x2
        x2 = x
        #yi = y(i-2) - q*y(i-1)
        y = y1 - q*y2
        y1 = y2
        y2 = y
    return(x1, y1, a)