Modular Multiplicative Inverse

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 equatiin 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 a_n·x_n + m_n·y_n = 1, starting with n = 0
2. Compute the quotient, q_n, and the remainder, r_n, of a_n divided by m_n, and rewrite the original equation as (m_n·q_n + r_n)·x_n + m_n·y_n =1
   
3. Rearrange the terms in the equation to get m_n·(q_n·x_n+y_n) + r_n·x_n = 1
   
4. Do the changes of variables
   
   • x_n+1 = q_n·x_n+y_n,
     
   • y_n+1 = x_n,
     
   • a_n+1 = m_n, and
     
   • m_n+1 = r_n
     
5. Repeat the process, until a_n = 1, and m_n = 0.
   
6. Such equations has the trivial solution x_n = 1, y_n = 0.
   
7. Work your way back to x₀ and y₀, noting that
   
   • x_n-1 = y_n,
     
   • y_n-1 = x_n - q_n-1·y_n
     
8. x₀ is the number we are after.
   

There is a certain component of faith in believing that, eventually, a_n = 1 and m_n = 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 m_n = 0, but a_n = 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 x_n's and the y_n'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)²), i.e. time required to compute the modular multiplicative inverse is proportional to log(m)².