Modular Multiplicative Inverse

March 20, 2009

Now that we are all proficient with modular exponentiation thanks to my previous post, it is time to tackle a more complicated issue, still in the realm of modular arithmetic.
While addition, subtraction and multiplication are “compatible with the congruence relation” introduced by modular arithmetic, the same doesn’t happen with division. This means that solving a simple equation such as a·x = 1 (mod m) is anything but trivial. In fact, determining x = a-1 (mod m) that fulfills that equation is the whole object of this post…

There is, as always, the brute force approach, which I always like to visit first, so that I can later congratulate myself on how much the efficiency has improved with my deep insight. In this particular case, brute force means trying every number between 1 and m-1, and see if the m-modulo of its product with a is 1. :

def modInv1(a,m) :
    Computes the modular multiplicative inverse of a modulo m,
    using brute force
    a %= m
    for x in range(1,m) :
        if a*x%m == 1 :
        return x
    return None

Do notice that the possibility of no multiplicative inverse existing is contemplated in the code. This will actually be the case if a and m have any common factors, i.e. if they are not coprime. So the previous code could serve as a moron’s approach to determining if two numbers are relatively prime. Not that I recommend it, though, as it should be clear that the code I have just produced will be very inefficient for large values of m.

The Extended Euclidean Algorithm
As stated, we are after a number x that fulfills a·x = 1 (mod m). This can of course be written as well as a·x = 1 + m·y, which rearranges neatly into a·x – m·y = 1. Since x and y need not be positive, we can write it as well in the standard form, a·x + m·y = 1. If a and m are coprime, this is a particular instance of a·x + b·y = gcd(a,b), where gcd(a,b) is the greatest common divisor of a and b. This equation is known as Bézout’s identity, and can be efficiently solved using the extended Euclidean algorithm. To find a solution to our particular equation, this algorithm would proceed as follows:

  1. Rewrite the equation as an·xn + mn·yn = 1, starting with n = 0
  2. Compute the quotient, qn, and the remainder, rn, of an divided by mn, and rewrite the original equation as (mn·qn + rn)·xn + mn·yn =1
  3. Rearrange the terms in the equation to get mn·(qn·xn+yn) + rn·xn = 1
  4. Do the changes of variables
    • xn+1 = qn·xn+yn,
    • yn+1 = xn,
    • an+1 = mn, and
    • mn+1 = rn
  5. Repeat the process, until an = 1, and mn = 0.
  6. Such equations has the trivial solution xn = 1, yn = 0.
  7. Work your way back to x0 and y0, noting that
    • xn-1 = yn,
    • yn-1 = xn – qn-1·yn
  8. x0 is the number we are after.

There is a certain component of faith in believing that, eventually, an = 1 and mn = 0. It is in fact not always the case, as if a and m are not coprime, then the values we will arrive at are mn = 0, but an = gcd(a,m). As a matter of fact, the gcd is normally computed using the Euclidean algorithm, dropping the adjective extended, which basically does the same as above, but without the back-tracking of the xn‘s and the yn‘s. With all of this in mind, we can write two functions, one calling the other, as follows:

def extEuclideanAlg(a, b) :
    Computes a solution  to a x + b y = gcd(a,b), as well as gcd(a,b)
    if b == 0 :
        return 1,0,a
    else :
        x, y, gcd = extEuclideanAlg(b, a % b)
        return y, x - y * (a // b),gcd

def modInvEuclid(a,m) :
    Computes the modular multiplicative inverse of a modulo m,
    using the extended Euclidean algorithm
    x,y,gcd = extEuclideanAlg(a,m)
    if gcd == 1 :
        return x % m
    else :
        return None

The gains in speed are tremendous with this new approach. Again, write your own testing function as an exercise, but mine spits out the following results…

>>> testModInv(1234,5,10000)
Got 4 using method 1 in 0.0160000324249 sec.
Got 4 using Euclidean algorithm in 0.0309998989105 sec.
>>> testModInv(1234,56789,10000)
Got 31800 using method 1 in 59.4000101089 sec.
Got 31800 using Euclidean algorithm in 0.047000169754 sec.

So while brute force is good enough, or even the best choice, for very small values of m, as soon as it grows, the benefits of the more elaborate algorithm are all too obvious. According to wikipedia, the algorithm we have coded has complexity O(log(m)2), i.e. time required to compute the modular multiplicative inverse is proportional to log(m)2.


Modular Exponentiation

March 20, 2009

I have a post cooking on wheel factorization, yet another sieving algorithm for the computation of all primes below a given number n. It turns out that while implementing it, one comes across a number of not so obvious preliminary questions. It is a rule of this blog not to resort to any non standard library, and I feel it would be breaking the spirit of that rule to post code which hasn’t been thoroughly explained in advance. So in order not to make that soon-to-be post on prime numbers longer than anyone could bear, I find myself in need of covering a couple of preliminaries, and so here comes this thrilling post on modular exponentiation

What we are after is determining ab (mod m), or to put it in plain English, the remainder of a to the power of b, divided by m. Doesn’t seem much, does it? We could just go ahead and write:

def modExp1(a, b, m) :
    Computes a to the power b, modulo m, directly
    return a**b % m

The problems begin if a and b are both big. Python’s integers don’t overflow, and we all love it for it, but the memory juggling it has to do to accommodate a huge number in memory does come at a speed cost. Specially if m is small, it definitely seems a little overkill to keep a zillion digit number in memory, given that the answer will be smaller than m…

We can, and will, use modular arithmetic to our advantage. See, wikipedia says that “the congruence relation (introduced by modular arithmetic) on the integers is compatible with the multiplication operation of the ring of integers,” and if it’s on the internet it has to be true, right? Translated to plain English, the idea is that if we want to take the product of two numbers, a and b, modulo m, we can first multiply them, and then take the m-modulo of the result, or we can take the m-modulo of a and/or b, prior to multiplication, and the m-modulo of the result will still be the right answer. So the least a conscious programmer should do is to rewrite the previous function as:

def modExp2(a, b, m) :
    Computes a to the power b, modulo m, a little less directly
    return (a%m)**b % m

Of course, if b is large enough, you may still find yourself in the realm of zillion digit numbers, which is a place we don’t really want to go. It turns out then that it may be a good option to remember that exponentiation is just repeated multiplication. If we multiply a by itself b times, and take the m-modulo of the result after every multiplication, the larger number we will ever have in memory will be m2. It doesn’t really seem right to put a loop of b iterations into the code, since what we are worried about are large b’s, but lets do it anyway:

def modExp3(a, b, m) :
    Computes a to the power b, modulo m, by iterative multiplication
    ret = 1
    a %= m
    while b :
        ret *= a
        ret %= m
        b -= 1
    return ret

The only reason why this last function deserves a place here, is because it helps introduce the really cool algorithm that this whole post really is about, binary exponentiation. In John D. Cook’s blog there is a great explanation on how it could be implemented iteratively, based on the binary representation of the exponent, but I’ll go for a recursive version, which is essentially the same, but which I find easier for my brain’s wiring…

So what we are about to do is, when asked to compute ab, check if b is even, in which case we will compute it as the square of ab/2. If b is odd, we will compute it as a·ab-1. If b is 0, we will of course return 1. If we also throw in all the modulo m stuff I have been discussing previously, we finally arrive at THE function to do modular exponentiation:

def modExp(a, b, m) :
    Computes a to the power b, modulo m, using binary exponentiation
    a %= m
    ret = None
    if b == 0 :
        ret = 1
    elif b%2 :
        ret = a * modExp(a,b-1,m)
    else :
        ret = modExp(a,b//2,m)
        ret *= ret
    return ret%m

Now that we have all four versions of the same function written, we can do a little testing, I will leave writing the testing code as an exercise for the reader (I’ve been wanting to say this since my first year of university…), but when computing 1234 to the power 5678, modulo 90, a 1000 times with each function, these is the performance I get:

>>> testModExp(1234,5678,90,1000)
Got 46 using method 1 in 3.73399996758 sec.
Got 46 using method 2 in 0.546999931335 sec.
Got 46 using method 3 in 1.25 sec.
Got 46 using THE method in 0.0150001049042 sec.

Results which really don’t require much commenting…

EDIT: There is some very interesting information in the comments: python comes with a built-in pow function, which can take two or three arguments, the third being, precisely, the modulo. It probably implements an algorithm similar to what has been presented, but is more than an order of magnitude faster.