# [Edu-sig] More useless Python (CS/number theory)

Kirby Urner pdx4d@teleport.com
Tue, 28 Aug 2001 09:57:54 -0700

```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
(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

```