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

[Kirby Urner]

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