Extended Euclidean Algorithm

Euclidean Algorithm

The original Euclidean Algorithm computes $\gcd(a,b)$ and looks like this:

Copy
ll euclid(ll a, ll b) {
	while (b > 0) {
		ll k = a / b;
		a -= k * b;
		swap(a, b);
	}
	return a;
}

Copy
static long euclid(long a, long b) {
	while (b > 0) {
		long k = a / b;
		a -= k * b;

		long temp = b;
		b = a;
		a = temp;
	}
	return a;
}

Copy
def euclid(a: int, b: int) -> int:
	assert (
		int(a) > 0 and int(b) > 0
	), "Arguments must be positive, non-zero numeric values."

	while b > 0:
		k = a // b
		# subtract multiples of one equation from the other.
		a -= b * k
		a, b = b, a

	return a

Extended Euclidean Algorithm

The extended Euclidean algorithm computes integers $x$ and $y$ such that

We can slightly modify the version of the Euclidean algorithm given above to return more information!

Copy
array<ll, 3> extend_euclid(ll a, ll b) {
	// we know that (1 * a) + (0 * b) = a and (0 * a) + (1 * b) = b
	array<ll, 3> x = {1, 0, a};
	array<ll, 3> y = {0, 1, b};

	// run extended Euclidean algo
	while (y[2] > 0) {
		// keep subtracting multiple of one equation from the other
		ll k = x[2] / y[2];
		for (int i = 0; i < 3; i++) { x[i] -= k * y[i]; }
		swap(x, y);
	}
	return x;  // x[0] * a + x[1] * b = x[2], x[2] = gcd(a, b)
}

Copy
static long[] extendEuclid(long a, long b) {
	// we know that 1 * a + 0 * b = a and 0 * a + 1 * b = b
	long[] x = {1, 0, a};
	long[] y = {0, 1, b};

	// run extended Euclidean algo
	while (y[2] > 0) {
		// keep subtracting multiple of one equation from the other
		long k = x[2] / y[2];
		for (int i = 0; i < 3; i++) { x[i] -= k * y[i]; }

Copy
def extend_euclid(a: int, b: int) -> list[int]:
	assert (
		int(a) > 0 and int(b) > 0
	), "Arguments must be positive, non-zero numeric values."

	# we know that 1 * a + 0 * b = a and 0 * a + 1 * b = b.
	x_arr = [1, 0, int(a)]
	y_arr = [0, 1, int(b)]
	q = -1

Recursive Version

Copy
ll euclid(ll a, ll b) {
	if (b == 0) { return a; }
	return euclid(b, a % b);
}

Copy
static long euclid(long a, long b) {
	if (b == 0) { return a; }
	return euclid(b, a % b);
}

Copy
def euclid(a: int, b: int) -> int:
	"""Recursive Euclidean GCD."""
	return a if b == 0 else euclid(b, a % b)

becomes

Copy
pl extend_euclid(ll a, ll b) {  // returns {x,y}
	if (b == 0) { return {1, 0}; }
	pl p = extend_euclid(b, a % b);
	return {p.s, p.f - a / b * p.s};
}

Copy
static long[] extendEuclid(long a, long b) {  // returns {x,y}
	if (b == 0) { return new long[] {1, 0}; }
	long[] p = extendEuclid(b, a % b);
	return new long[] {p[1], p[0] - a / b * p[1]};
}

Copy
def extend_euclid(a: int, b: int) -> list[int]:
	if not b:
		return [1, 0]
	p = extend_euclid(b, a % b)
	return [p[1], p[0] - (a // b) * p[1]]

The pair will equal to the first two returned elements of the array in the iterative version. Looking at this version, we can prove by induction that when $a$ and $b$ are distinct positive integers, the returned pair $(x,y)$ will satisfy $|x|\le \frac{b}{2\gcd(a,b)}$ and $|y|\le \frac{a}{2\gcd(a,b)}$. Furthermore, there can only exist one pair that satisfies these conditions!

Note that this works when $a,b$ are quite large (say, $\approx 2^{60}$) and we won't wind up with overflow issues.

Application - Modular Inverse

It seems that when multiplication / division is involved in this problem, $n^2 < \texttt{LLONG\_MAX}$.

Copy
ll inv_general(ll a, ll b) {
	array<ll, 3> x = extend_euclid(a, b);
	assert(x[2] == 1);  // gcd must be 1
	return x[0] + (x[0] < 0) * b;
}

Copy
static long invGeneral(long a, long b) {
	long[] x = extendEuclid(a, b);
	assert (x[2] == 1);  // gcd must be 1
	return x[0] + (x[0] < 0 ? 1 : 0) * b;
}

Copy
def inv_general(a: int, b: int) -> int:
	"""Returns the modular inverse of two positive, non-zero integer values."""
	arr = extend_euclid(a, b)
	assert arr[2] == 1, "GCD must be 1."
	return arr[0] + (arr[0] < 0) * b

Application - Chinese Remainder Theorem

For Two Congruences

Explanation

The Chinese Remainder Theorem - also referred to as CRT - yields a unique solution to a system of simultaneous modular congruences with pairwise coprime moduli. More specifically, CRT determines a number $x$ that when divided by the given divisors leaves the given remainders. Mathematically, it solves the following system of modular congruences:

The above system has exactly one solution modulo $M$.

We'll focus on the particular case for $k=2$; for two coprime moduli:

• Let $n_1$ and $n_2$be the modular inverse of $m_1$ modulo $m_2$ and the modular inverse of $m_2$

modulo $m_1$, respecively, thus the properties hold: $n_1 \cdot m_1 \equiv 1 \pmod{m_2}$ and $n_2 \cdot m_2 \equiv 1 \pmod{m_1}$. It's worth mentioning that $n_1$ and $n_2$ always exists because $\gcd(m_1,m_2)=1$. We can define a solution $x$ as:

It satisfies both equations:

• Modulo $m_1$ we have $x \equiv a_1 \cdot m_2 \cdot n_2 \pmod{m_1} \iff x \equiv a_1 \pmod{m_1}$ since $m_2 \cdot n_2 \equiv 1 \pmod{m_1}$
• Modulo $m_2$ we have $x \equiv a_2 \cdot m_1 \cdot n_1 \pmod{m_2} \iff x \equiv a_2 \pmod{m_2}$ since $m_1 \cdot n_1 \equiv 1 \pmod{m_2}$

Now that we have a solution $x \pmod{M}$ it only remains to show that it is indeed the only valid solution modulo $m_1 \cdot m_2$. Assuming that there is a different valid solution $y$ we have: $m_1 | y-x$, $m_2 | y-x$, $m_3 | y-x$, $\ldots$, $m_k | y-x$, then it follows that $m_1 \cdot m_2 \cdots m_k$ divides $y-x$, since $m_1, m_2, \ldots, m_k$ are relatively prime. This means that: $y \equiv x \pmod{m_1 \cdot m_2 \cdots m_k}$.

Implementation

Time Complexity: $\mathcal{O}(T \cdot \log{MN})$

Copy
#include <algorithm>
#include <array>
#include <iostream>
#include <numeric>
#include <vector>

using namespace std;

typedef long long ll;

For Several Congruences

Using the same notations as in $(1)$, consider $n_i = \frac{M}{m_i}$ - the product of all moduli excluding $m_i$ - and $n_i' \equiv n_i^{-1} \pmod{m_i}$. By similar arguments as before, the unique solution mod $M$ is: