Greatest Common Divisor (GCD)

GCD of two integers is the largest number which evenly divides both the numbers. For example: GCD(24, 16) is 8 and GCD(-4, -8) is 4.

A simple algorithm to find GCD can be:

int gcd(int a, int b) {
  for (int i = min(a, b); i >= 1; i--) {
    if (a % i == 0 && b % i == 0) return i;
  }
  return 1;
}

However, there exist a significantly faster algorithm to find the same called Euclidean algorithm.

Euclidean Algorithm#

The algorithm states:

1. If B = 0, then GCD(A, B) = A.
2. Otherwise GCD(A, B) = GCD(B, A % B).

This algorithm reduces the problem to a simpler one in each iteration. Overall the time complexity is log(max(A, B)). The code for the same is:

int gcd(int a, int b) {
  if (b == 0) return a;
  return gcd(b, a % b);
}

Proof#

1. If B = 0, then GCD(A, B) = A.

Since 0 is divisible by all integers (including A). The GCD(A, B) will simply be A.

2. GCD(A, B) = GCD(B, A % B)

To prove this we will first prove the relation GCD(A, B) = GCD(B, A - B).

Let C = A - B.

We know A = x * GCD(A, B) and B = y * GCD(A, B) => C = A - B = (x - y) * GCD(A, B)

GCD(A, B) is a common divisor of B and C (which is smaller or equal to GCD(B, C)).

Similarly, B = x * GCD(B, C) and C = y * GCD(B, C) => A = B + C = (x + y) * GCD(B, C)

GCD(B, C) is a common divisor of A and B (which is smaller or equal to GCD(A, B)).

Since GCD(A, B) <= GCD(B, C) and GCD(B, C) <= GCD(A, B). Hence, GCD(A, B) = GCD(B, C).

Since the order of operation does not matter, GCD(A, B) = GCD(A - B, B).

We can repeatedly apply GCD(A, B) = GCD(A - B, B) to get GCD(A - Q * B, B) where A = Q * B + R

Hence, GCD(A, B) = GCD(A % B, B) = GCD(A, A % B)

Extended Euclidean Algorithm#

In addition to the greatest common divisor of integers a and b, this algorithm also finds x and y such that a * x + b * y = GCD(a, b) (x and y is useful for finding the inverse of a integer).

Base Case: If b = 0 then GCD(a, b) = a. Hence x and y would be (1, 0) respectively.

If we are able to find (x, y) for GCD(a, b) from (x', y') for GCD(b, a % b) then we can easily be able to propagate those values up the recursion from the base case.

We have b * x' + (a - (a / b) * b) * y' = g and a * x + b * y = g => y' * a + b * (x' - y' * (a / b)) = a * x + b * y

Comparing both sides we get, x = y' and y = x' - y' * (a / b)

Code#

int extended_gcd(int a, int b, int& x, int& y) {
  if (b == 0) {
    x = 1, y = 0;
    return a;
  }
  int x1, y1;
  int res = extended_gcd(b, a % b, x1, x2);
  x = x1, y = x1 - y1 * (a / b);
  return res;
}

Modular Multiplicative Inverse#

Modular multiplicative inverse of an integer a wrt to modulus m is an integer x such that (a * x) % m = 1. In other words remainder of (a * x) when divided by m is 1.

This implies that GCD(a * x, m) = 1.

Also, an inverse only exists when a and m are coprime i.e. GCD(a, m) = 1.

Code:

int mod_inv(int a, int m) {
  int x, y;
  extended_gcd(a, m, x, y);
  return (x % m + m) % m;
}