Modular Multiplicative Inverse

Copy from this link http://www.doc.ic.ac.uk/~mrh/330tutor/ch03.html

http://rosettacode.org/wiki/Modular_inverse

https://comeoncodeon.wordpress.com/2011/10/09/modular-multiplicative-inverse/

This can be done by running an extended version of Euclidean algorithm which tracks x when computing the gcd. In the extended Euclidean algorithm, we first initialise x1 =0 and x2 =1, then in the following steps, compute xi = xi-2xi-1qi-2 where qi-2 is the quotient computed in step i-2. When the remainder becomes 0, continue the calculation of x for one more round. The final x is the inverse. Here is an example that shows how to find the inverse of 15 when the modulus is 26:

  • step 1: 26÷15, quotient q1= 1, remainder = 11, x1 = 0
  • step 2: 15÷11, quotient q2 = 1, remainder = 4, x2 = 1
  • step 3: 11÷4, quotient q3 = 2, remainder = 3, x3 = x1-x2q1 = 0- 1*1 = -1
  • step 4: 4÷3, quotient q4 = 1, remainder = 1, x4 = x2-x3q2 = 1- (-1)*1 = 2
  • step 5: 3÷1, quotient q5 = 3, remainder = 0, x5 = x3-x4q3 = -1- 2*2 = -5
  • step 6: x6 = x4-x5q4 = 2- (-5)*1 = 7

To verify, 157 = 105 = 426+1, so 15*7 ≡ 1 mod 26, which means 7 is the multiplicative inverse of 15 under modulo 26.

#include <iostream>
 using namespace std;
 
int mul_inv(int a, int b)
{
    int b0 = b, t, q;
    int x0 = 0, x1 = 1;
    if (b == 1) return 1;
    while (a > 1) {
        q = a / b;
        t = b, b = a % b, a = t;
        t = x0, x0 = x1 - q * x0, x1 = t;
    }
    if (x1 < 0) x1 += b0;
    return x1;
}
 
int main(void) {
    cout<<mul_inv(42, 2017)<<endl;
    return 0;
}


1. Brute Force
We can calculate the inverse using a brute force approach where we multiply a with all possible values x and find a x such that ax \equiv 1 \pmod{m}. Here’s a sample C++ code:

1
2
3
4
5
6
int modInverse(int a, int m) {
    a %= m;
    for(int x = 1; x < m; x++) {
        if((a*x) % m == 1) return x;
    }
}

The time complexity of the above codes is O(m).

2. Using Extended Euclidean Algorithm
We have to find a number x such that a·x = 1 (mod m). This can be written as well as a·x = 1 + m·y, which rearranges 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.

In number theory, Bézout’s identity for two integers a, b is an expression ax + by = d, where x and y are integers (called Bézout coefficients for (a,b)), such that d is a common divisor of a and b. If d is the greatest common divisor of a and b then Bézout’s identity ax + by = gcd(a,b) can be solved using Extended Euclidean Algorithm.

The Extended Euclidean Algorithm is an extension to the Euclidean algorithm. Besides finding the greatest common divisor of integers a and b, as the Euclidean algorithm does, it also finds integers x and y (one of which is typically negative) that satisfy Bézout’s identity
ax + by = gcd(a,b). The Extended Euclidean Algorithm is particularly useful when a and b are coprime, since x is the multiplicative inverse of a modulo b, and y is the multiplicative inverse of b modulo a.

We will look at two ways to find the result of Extended Euclidean Algorithm.

Iterative Method
This method computes expressions of the form ri = axi + byi for the remainder in each step i of the Euclidean algorithm. Each successive number ri can be written as the remainder of the division of the previous two such numbers, which remainder can be expressed using the whole quotient qi of that division as follows:
r_i = r_{i-2} - q_i r_{i-1}.\,
By substitution, this gives:
r_i = (ax_{i-2} + by_{i-2}) - q_i (ax_{i-1} + by_{i-1}),\,which can be written
r_i = a(x_{i-2} - q_i x_{i-1}) + b(y_{i-2} - q_i y_{i-1}).\,
The first two values are the initial arguments to the algorithm:
r_1 = a = a\times1 + b\times0\,
r_2 = b = a\times0 + b\times1.\,
So the coefficients start out as x1 = 1, y1 = 0, x2 = 0, and y2 = 1, and the others are given by
x_i = x_{i-2} - q_i x_{i-1},\,
y_i = y_{i-2} - q_i y_{i-1}.\,
The expression for the last non-zero remainder gives the desired results since this method computes every remainder in terms of a and b, as desired.

So the algorithm looks like,

  1. Apply Euclidean algorithm, and let qn(n starts from 1) be a finite list of quotients in the division.
  2. Initialize x0, x1 as 1, 0, and y0, y1 as 0,1 respectively.
    1. Then for each i so long as qi is defined,
    2. Compute xi+1 = xi-1qixi
    3. Compute yi+1 = yi-1qiyi
    4. Repeat the above after incrementing i by 1.
  3. The answers are the second-to-last of xn and yn.
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
/* This function return the gcd of a and b followed by
    the pair x and y of equation ax + by = gcd(a,b)*/
pair<int, pair<int, int> > extendedEuclid(inta, intb) {
    intx = 1, y = 0;
    intxLast = 0, yLast = 1;
    intq, r, m, n;
    while(a != 0) {
        q = b / a;
        r = b % a;
        m = xLast - q * x;
        n = yLast - q * y;
        xLast = x, yLast = y;
        x = m, y = n;
        b = a, a = r;
    }
    returnmake_pair(b, make_pair(xLast, yLast));
}
intmodInverse(inta, intm) {
    return(extendedEuclid(a,m).second.first + m) % m;
}
Advertisements