
Here's a version of the Extended Euclidean Algorithm, with related inverse and Chinese Remainder Theorem functions. Very old hat stuff in CS, but that's why important. Note I use the new //, which is available in 2.2a2 w/o import from __future__. The EEA returns a 3-tuple, the 0th element being the gcd(a,b), the next two being s,t such that gcd(a,b) = s*a + t*b, i.e. two integers which bring the initial a,b within the gcd of one another. Example:
eea(17,25) # gcd(a,b) = 1 [1, 3, -2] 3*17 - 2*25 # s = 3, t = -2 1
eea(29834,8282) # gcd(a,b) = 2 [2, -88, 317] -88*29834 + 317*8282 # s = -88, t = 317 2
With EEA, you can also find the inverse of x mod n, provided gcd(x,n)=1. The inverse is that number which, multiplying x, gives a remainder of 1 when divided by n. Example: find x such that (x * 7) mod 4 = 1. Computing: >>> inverse(7,4) 3 i.e. (3*7) mod 4 = 21 mod 4 = integer remainder of 21/4 = 1. The Chinese Remainder Theorem states that if I give you a smattering of divisors, all coprime to each other, and tell you what the remainders for each of them is, then you can tell me a number which meets my specifications. Example: Your divisors are 7,11 and 15. The respective remainders are 2, 3 and 0. What's a number that works? Computing: >>> crt([7,11,15],[2,3,0]) 135 Check: >>> 135 % 7 2 >>> 135 % 11 3 >>> 135 % 15 0 Works. The also EEA comes up in CS in connection with RSA. You need it to find secret pair d, matched with e, such that (d*e) mod (p-1)(q-1) = 1, where p,q are the two humongous primes used to generate p*q = n, the public key. There's a more generalized form of the CRT which allows the divisors to not necessarily be coprime to start with (but with another stipulation about the remainders in that case). See Knuth, Art of Computer Programming, Vol 3, pg. 292. Tweaking the code below to accommodate this case might be an exercise. HINT: you'd probably want an lcm function. Kirby ===================== Here's the code: from operator import mod def eea(a,b): """Extended Euclidean Algorithm for GCD""" v1 = [a,1,0] v2 = [b,0,1] while v2[0]<>0: p = v1[0]//v2[0] # floor division v2, v1 = map(sub,v1,[p*vi for vi in v2]), v2 return v1 def inverse(m,k): """ Return b such that b*m mod k = 1, or 0 if no solution """ v = eea(m,k) return (v[0]==1)*(v[1] % k) def crt(ms,as): """ Chinese Remainder Theorem: ms = list of pairwise relatively prime integers as = remainders when x is divided by ms (ai is 'each in as', mi 'each in ms') The solution for x modulo M (M = product of ms) will be: x = a1*M1*y1 + a2*M2*y2 + ... + ar*Mr*yr (mod M), where Mi = M/mi and yi = (Mi)^-1 (mod mi) for 1 <= i <= r. """ M = reduce(mul,ms) # multiply ms together Ms = [M/mi for mi in ms] # list of all M/mi ys = [inverse(Mi, mi) for Mi,mi in zip(Ms,ms)] # uses inverse,eea return reduce(add,[ai*Mi*yi for ai,Mi,yi in zip(as,Ms,ys)]) % M