Entry

Euclid's algoritm for greatest common denominator

Jul 5th, 2000 10:01
Nathan Wallace, Hans Nowak, Snippet 208, Tim Peters

"""
Packages: maths
"""
#!/usr/local/bin/python
#!c:\python\python.exe
import sys
import os
import string
import math
def exgcd ( u, v, pr = 0 ) : 
    """An implementation of the Extended Euclidean algorithm for finding 
    greatest common denominator, called `The Pulverizer', or 'kuttaka'.  
    (See Knuth, V2 (3d Ed.), 4.5.2.  Here is Knuth's description of the 
    Pulverizer:
    This extension of Euclid's algorithm can be described conveniently in 
    vector notation:
    Algorithm X (Extended Euclid's algorithm). Given nonnegative integers u
    and v, this algorithm determines a vector (u1, u2, u3) such that uu1 +
    uu2 = u3 = gcd(u, v). The computation makes use of auxiliary vectors
    (v1, v2, v3), (t1, t2, t3); all vectors are manipulated in such a way
    that the relations 
    ut1 + vt2 = t3; uu1 + vu2 = u3; uv1 + vv2 = v3; 
    hold throughout the calculation. 
    X1. [Initialize.]
    Set (u1, u2, u3)  (1, 0, u), (v1, v2, v3)  (0, 1, v). 
    X2. [Is v3 = 0?] If v3 = 0, the algorithm terminates. 
    X3. [Divide, subtract.] Set q  u3/v3, and then set 
    (t1, t2, t3)  (u1, u2, u3) - (v1, v2, v3)q, 
    (u1, u2, u3)  (v1, v2, v3), (v1, v2, v3)  (t1, t2, t3). 
    Return to step X2. 
    For example, let u = 40902, v = 24140. At step X2 we have 
    q    u1    u2    u3      v1   v2     v3
    0     1    0  40902      0    1  24140
    1     0    1  24140      1   -1  16762
    1     1   -1  16762     -1    2   7378
    2    -1    2   7378      3   -5   2006
    3     3   -5   2006    -10   17   1360
    1   -10   17   1360     13  -22    646
    2    13  -22    646    -36   61     68
    9   -36   61     68    337 -571     34
    2   337 -571     34   -710 1203      0
    The solution is therefore
    337 · 40902 -571 · 24140 = 34 = gcd(40902, 24140). 
    If the pr argument is non-zero, then exgcd() will print the intermediate 
    values for each step."""
    u1, u2, u3 = 1, 0, u
    v1, v2, v3 = 0, 1, v
    q = 0
    if pr != 0 :
        print "q %2d   u1 %4d u2 %4d u3 %6d   v1 %4d v2 %4d v3 %6d" % \
	    ( q, u1, u2, u3, v1, v2, v3 )
    while v3 > 0 :
        q = u3 / v3
        t1, t2, t3 = ( u1 - ( v1 * q ) ), ( u2 - ( v2 * q ) ), ( u3 - ( v3 * q ) )
        u1, u2, u3 = v1, v2, v3
        v1, v2, v3 = t1, t2, t3
	if pr != 0 :
            print "q %2d   u1 %4d u2 %4d u3 %6d   v1 %4d v2 %4d v3 %6d" % \
	        ( q, u1, u2, u3, v1, v2, v3 )
    return u3, u1, u2
if __name__ == "__main__" :
    u = 40902
    v = 24140
    if len ( sys.argv ) > 1 :
	if len ( sys.argv ) > 2 :
	    u = string.atoi ( sys.argv[ 1 ] )
	    v = string.atoi ( sys.argv[ 2 ] )
	else :
	    print "Two integer arguments required."
	    sys.exit ( -1 )
    u3, u1, u2 = exgcd ( u, v, 1 )
    print "The answer is %d" % ( u3 )
    print "u1 %d * u %d        u2 %d * v %d" % ( u1, u, u2, v )
    print " = %d                = %d" % ( u1 * u, u2 * v )
    print "         a = %d" % ( ( u1 * u ) + ( u2 * v ) )
    print "u3 %d u1 %d u2 %d" % ( u3, u1, u2 )