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

求最大公约数,一般采用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)
Python

这次主要是说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*xi + b*yi
gcd(b, a % b) = b*xi+1 + (a - [a/b]*b)*yi+1    (q = [a/b])
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

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)
Python

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

r1 = a mod b

r2 = b mod r1

r3 = r1 mod r2

...

ri = ri-2 mod ri-1   即(1) ri = ri-2 - q*ri-1   (q = [ri-2 / ri-1])

(令r-1 = a, r0 = b)

因为 a*x + b*y = gcd(a, b)即有,a*xi + b*yi = ri   (中间的每一步)

代入(1)式有,a*xi + b*yi = (a*xi-2 + b*yi-2 ) - q * (a*xi-1 + b*yi-1 ) = a * (xi-2 - q*xi-1) + b * (yi-2 - q*yi-1)

则有,xi  = xi-2 - q*xi-1,yi = yi-2 - q*yi-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)
Python

 

posted @ 2013-11-04 14:50  7hat  阅读(4590)  评论(0编辑  收藏  举报