Algoritmo de Euclides extendido

Algoritmo de Euclides

Recursos
Fuente	Recurso	Notas
cp-algo	Euclidean

El algoritmo de Euclides original calcula $\gcd(a,b)$ y se ve así:

C++JavaPython


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


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;
}


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

Algoritmo de Euclides extendido

Recursos
Fuente	Recurso	Notas
cp-algo	Extended Euclidean
Wikipedia	Extended Euclidean
cp-algo	Linear Diophantine Equation

El algoritmo de Euclides extendido calcula enteros $x$ e $y$ tales que

ax+by=\gcd(a,b)

¡Podemos modificar ligeramente la versión del algoritmo de Euclides dada arriba para devolver más información!

C++JavaPython


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)
}


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]; }
 
		long[] temp = x;
		x = y;
		y = temp;
	}
	return x;  // x[0] * a + x[1] * b = x[2], x[2] = gcd(a, b)
}


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
 
	while y_arr[2] > 0:  # run extended Euclidean algo.
		q = x_arr[2] // y_arr[2]
		for i in range(3):
			# keep subtracting multiple of one equation from the other.
			x_arr[i] -= y_arr[i] * q
		x_arr, y_arr = y_arr, x_arr
 
	return x_arr  # (x[0] * a) + (x[1] * b) = x[2], x[2] = gcd(a, b)

Versión recursiva

C++JavaPython


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


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


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

se convierte en

C++JavaPython


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};
}


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]};
}


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]]

El par será igual a los primeros dos elementos devueltos del arreglo en la versión iterativa. Mirando esta versión, podemos demostrar por inducción que cuando $a$ y $b$ son enteros positivos distintos, el par devuelto $(x,y)$ cumplirá $|x|\le \frac{b}{2\gcd(a,b)}$ y $|y|\le \frac{a}{2\gcd(a,b)}$ . ¡Además, solo puede existir un par que cumpla estas condiciones!

Observar que esto funciona cuando $a,b$ son bastante grandes (digamos, $\approx 2^{60}$ ) y no tendremos problemas de desbordamiento.

Aplicación - Inverso modular

Recursos
Fuente	Recurso	Notas
cp-algo	Modular Inverse

Hecho	Fuente	Nombre	Dificultad	Tags	Solución
	Kattis	Modular Arithmetic	N/A		en el módulo

Parece que cuando hay multiplicación / división involucrada en este problema, $n^2 < \texttt{LLONG\\_MAX}$ .

C++JavaPython


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;
}


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;
}


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

Aplicación - Teorema Chino del Resto

Recursos
Fuente	Recurso	Notas
Stanford	The Chinese Remainder Theorem	Explicación en profundidad
cp-algo	Chinese Remainder Theorem
CMU	Chinese Remainder Theorem

Hecho	Fuente	Nombre	Dificultad	Tags	Solución
	Kattis	Chinese Remainder	N/A		en el módulo

Para dos congruencias

Explicación

El Teorema Chino del Resto — también llamado CRT — produce una solución única a un sistema de congruencias modulares simultáneas con módulos coprimos dos a dos. Más específicamente, el CRT determina un número $x$ que al ser dividido por los divisores dados deja los restos dados. Matemáticamente, resuelve el siguiente sistema de congruencias modulares:

x \equiv a_1 \pmod{m_1} \\ x \equiv a_2 \pmod{m_2} \\ x \equiv a_3 \pmod{m_3} \\ \vdots\\ x \equiv a_k \pmod{m_k} \\ \tag{1}

M=m_1\cdot m_2 \cdot m_3 \cdots m_k

El sistema de arriba tiene exactamente una solución módulo $M$ .

Nos enfocaremos en el caso particular $k=2$ ; para dos módulos coprimos:

x \equiv a_1 \pmod{m_1} \\ x \equiv a_2 \pmod{m_2} \\ \tag{2}

M=m_1 \cdot m_2

Sean $n_1$ y $n_2$ el inverso modular de $m_1$ módulo $m_2$ y el inverso modular de $m_2$ módulo $m_1$ , respectivamente, de modo que se cumplen las propiedades: $n_1 \cdot m_1 \equiv 1 \pmod{m_2}$ y $n_2 \cdot m_2 \equiv 1 \pmod{m_1}$ . Vale mencionar que $n_1$ y $n_2$ siempre existen porque $\gcd(m_1,m_2)=1$ . Podemos definir una solución $x$ como:

x = a_1 \cdot m_2 \cdot n_2 + a_2 \cdot m_1 \cdot n_1 \pmod{m_1\cdot m_2}

Satisface ambas ecuaciones:

Módulo $m_1$ tenemos $x \equiv a_1 \cdot m_2 \cdot n_2 \pmod{m_1} \iff x \equiv a_1 \pmod{m_1}$ ya que $m_2 \cdot n_2 \equiv 1 \pmod{m_1}$
Módulo $m_2$ tenemos $x \equiv a_2 \cdot m_1 \cdot n_1 \pmod{m_2} \iff x \equiv a_2 \pmod{m_2}$ ya que $m_1 \cdot n_1 \equiv 1 \pmod{m_2}$

Ahora que tenemos una solución $x \pmod{M}$ solo queda mostrar que es en efecto la única solución válida módulo $m_1 \cdot m_2$ . Asumiendo que hay una solución válida distinta $y$ tenemos: $m_1 | y-x$ , $m_2 | y-x$ , $m_3 | y-x$ , $\ldots$ , $m_k | y-x$ , entonces se sigue que $m_1 \cdot m_2 \cdots m_k$ divide a $y-x$ , ya que $m_1, m_2, \ldots, m_k$ son relativamente primos. Esto significa que: $y \equiv x \pmod{m_1 \cdot m_2 \cdots m_k}$ .

Implementación

Complejidad temporal: $\mathcal{O}(T \cdot \log{MN})$

C++


#include <algorithm>
#include <array>
#include <iostream>
#include <numeric>
#include <vector>
 
using namespace std;
 
typedef long long ll;
 
// Extended Euclidean Algorithm
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};
 
	while (y[2] > 0) {
		// keep subtracting the 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)
}
 
// The Modular Inverse
ll inv(ll x, ll y) {
	ll m = extend_euclid(x, y)[0] % y;
	return (m < 0 ? m + y : m);
}
 
// General code for multiple congruences in a system
ll crt(const vector<ll> &remainders, const vector<ll> &moduli) {
	ll MOD = accumulate(moduli.begin(), moduli.end(), 1LL, multiplies<ll>());
	ll x = 0;
	for (int i = 0; i < (int)moduli.size(); i++) {
		ll a = remainders[i] * inv(MOD / moduli[i], moduli[i]) % moduli[i];
		x = (x + a * (MOD / moduli[i])) % MOD;
	}
	return x;
}
 
int main() {
	int tests;
	cin >> tests;
	for (int t = 0; t < tests; t++) {
		vector<ll> r(2), m(2);
		cin >> r[0] >> m[0] >> r[1] >> m[1];
		cout << crt(r, m) << ' ' << m[0] * m[1] << '\n';
	}
}

Para varias congruencias

Usando las mismas notaciones que en $(1)$ , consideremos $n_i = \frac{M}{m_i}$ — el producto de todos los módulos excluyendo $m_i$ — y $n_i' \equiv n_i^{-1} \pmod{m_i}$ . Por argumentos similares a los de antes, la solución única módulo $M$ es:

x = \sum_{i=1}^{n}{a_i \cdot n_i \cdot n_i'} \pmod{M}