扩展欧几里得算法（C++）

扩展欧几里得算法

前置知识：裴蜀定理

裴蜀定理

\qquad 若 $a,b$ 是整数,且 $(a,b)=d$，那么对于任意的整数 $x,y,ax+by$ 都一定是 $d$ 的倍数。特别地，一定存在整数 $x,y$，使 $ax+by=d$ 成立，（x,y不唯一）。

\qquad 若 $x_0,y_0$ 满足 $ax+by=d$，则满足 $\{\begin{array}{l}x=x_0-\frac{b}{d}\cdot k\\ y=y_0+\frac{a}{d}\cdot k\end{array},k\in \mathbb{Z}$，的 $x$ 和 $y$ 也均满足上式。且这些是 $x,y$ 所有可能的解。

作用

\qquad 扩展欧几里得算法就是用来求解这个特殊的 $x,y$.

推导证明

$$\begin{gathered} amodb=a-\lfloor \frac{a}{b}\rfloor \cdot b \\ d=(a,b)=(b,amodb)=\left(a-\lfloor \frac{a}{b}\rfloor \cdot b\right)x+by \\ d=ax+b\left(y-\lfloor \frac{a}{b}\rfloor x\right) \end{gathered}$$

\qquad 这里的 $y-\lfloor \frac{a}{b}\rfloor x$ 就是下面代码第9行 y -= a / b * x; 的由来。

代码

#include <iostream>
using namespace std;
int exgcd(int a, int b, int &x, int &y) {		//以传址的方式返回两个参数
	if (!b) {
		x = 1, y = 0;
		return a;
	}
	int d = exgcd(b, a % b, y, x);
	y -= a / b * x;
	return d;
}
int main() {
	int a, b, x, y;
	scanf("%d%d", &a, &b);
	exgcd(a, b, x, y);
	printf("%d %d\n", x, y);
	return 0;
}

应用

1.求逆元

前置知识：逆元

\qquad 若整数 $b$，$m$ 互质，并且对于任意的整数 $a$，如果满足 $b|a$，则存在一个整数 $x$，使得 $\frac{a}{b}\equiv a\times x(modm)$，则称 $x$ 为 $b$ 的模 $m$ 乘法逆元，记为 $b^{-1}(modm)$。

也就是说，求得一个 $x$，使得 $b\cdot x\equiv 1(modm)$ 成立。

思路

\qquad 这个问题可以看作求一元线性同余方程的一种特殊情况。
$$\begin{gathered} \exists x\in \mathbb{Z},使得ax\equiv 1(modm)成立 \\ \Updownarrow \\ \exists x,y\in \mathbb{Z},使得ax=my+1成立 \\ \Updownarrow \\ 令y^{\prime }=-y,\exists x,y^{\prime }\in \mathbb{Z},使得ax+m{y}^{\prime }=1成立 \end{gathered}$$
\qquad 根据裴蜀定理可知，当且仅当 $(a,m)=1$，即 $a$ 与 $m$ 互质时，此线性同余方程有解，即 $a$ 有存在逆元。

\qquad 利用扩展欧几里得算法求得的 $x$ 即为 $a$ 的逆元。

代码

#include <iostream>
using namespace std;
int exgcd(int a, int b, int &x, int &y) {		//以传址的方式返回两个参数
	if (!b) {
		x = 1, y = 0;
		return a;
	}
	int d = exgcd(b, a % b, y, x);
	y -= a / b * x;
	return d;
}
int main() {
	int a, m, x, y;
	scanf("%d%d", &a, &m);
	int d = exgcd(a, m, x, y);
	if (d == 1)
		printf("%d\n", x < 0 ? x + m : x);		//为了防止出现负数，等价于：printf("%d\n", (x + m) % m);
	else
		puts("no exist");
	return 0;
}

2.求解线性同余方程

例题

878. 线性同余方程 - AcWing题库

\qquad 给定 $n$ 组数据 $a_i,b_i,m_i$，对于每组数求出一个 $x$ ，使其满足 $a_i\ast x_i\equiv b_i(mod{m}_i)$，如果无解，则输出 impossible.

数据范围

\qquad $1\le n\le 10^5,1\le a_i,b_i,m_i\le 2\times 10^9$

输入样例

2
2 3 6
4 3 5

输出样例

impossible
2

思路

\qquad 这就是上文求逆元方法的扩展和延申。
$$\begin{gathered} \exists x\in \mathbb{Z},使得ax\equiv b(modm)成立 \\ \Updownarrow \\ \exists x,y\in \mathbb{Z},使得ax=my+b成立 \\ \Updownarrow \\ 令y^{\prime }=-y,\exists x,y^{\prime }\in \mathbb{Z},使得ax+m{y}^{\prime }=b成立 \end{gathered}$$
\qquad 根据裴蜀定理可知，记 $(a,m)=d$ 当且仅当 $d|b$ 时，此线性同余方程有解。

\qquad 要求 $x$ ，则需用扩展欧几里得算法算出满足 $a{x}_0+m{y}_0=d$ 式子的 $x_0$，可得 $x=x_0\cdot \frac{b}{d}$ .

举个例子：

要解 $4x\equiv 3(mod5)$ 这一方程：

1. 先判断 $(4,5)=1$ 是 $3$ 的因数，说明有解；
2. 再找出满足 $4x_0+5y_0=1$ 的一组解 $x_0=-1,y_0=1$.
3. 最后计算 $x$：$-1\times \frac{3}{1}=-3$，因此解为 $x\equiv -3(mod5)$.

代码

#include <iostream>
using namespace std;
typedef long long LL;
int exgcd(int a, int b, int &x, int &y) {
	if (!b) {
		x = 1, y = 0;
		return a;
	}
	int d = exgcd(b, a % b, y, x);
	y -= a / b * x;
	return d;
}
int main() {
	int n;
	scanf("%d", &n);
	while (n--) {
		int a, b, m, x, y;
		scanf("%d%d%d", &a, &b, &m);
		int d = exgcd(a, m, x, y);
		if (b % d)
			puts("impossible");
		else{
            LL res = (LL)x * b / d % m;
            printf("%lld\n", res < 0 ? res + m : res);
		}
	}
	return 0;
}

3.求解线性同余方程组

前置知识：中国剩余定理

\qquad 中国剩余定理又称孙子定理或中国余数定理，是数论中的一个关于一元线性同余方程组的定理，说明了一元线性同余方程组有解的准则以及求解方法。

经典题目表述

有物不知其数，三三数之剩二，五五数之剩三，七七数之剩二。问物几何？

\qquad 大意是：一个整数除以三余二，除以五余三，除以七余二，求这个整数。

\qquad 用现代数学的语言来说明，题目为：给出以下一元线性同余方程组，求满足下列同余式的整数 $x$.
$$\{\begin{array}{l}x\equiv a_1(mod{m}_1)\\ x\equiv a_2(mod{m}_2)\\ ⋮\\ x\equiv a_n(mod{m}_n)\end{array}$$
其中，$m_1,m_2,\dots m_n$ 两两互质。

注意：这里两两互质是指任意两个数均互质，而不是这些数的最大公因数为1.

例如：{6，8，9} 这三个数的最大公因数为1，但是他们并非两两互质。

解题方法
$$\begin{gathered} 设M=\prod _{i=1}^n{m}_i,M_i=\frac{M}{m_i} \\ 用M_i^{-1}表示M_i在模m_i意义下的逆元，则解为： \\ x\equiv \sum _{i=1}^n{a}_i{M}_i{M}_i^{-1}(modM) \end{gathered}$$
\qquad 在此处需要用扩展欧几里得定理计算一个数字的逆元。由此我们就可以知道，这种方法需要 $m_1,m_2,\dots m_n$ 两两互质的原因就在于需要求每个 $M_i$ 在模 $m_i$ 下的逆元（逆元存在的充要条件就是 $M_i$ 和 $m_i$ 互质）。

举个例子：
$$\{\begin{array}{l}x\equiv 2(mod3)\\ x\equiv 3(mod5)\\ x\equiv 2(mod7)\end{array}$$
就拿“物不知数”这个经典题举例：

1. 首先判断 $3,5,7$ 之间确实两两互质，可以使用中国剩余定理；

2. 求出 $M=3\times 5\times 7=105$；

3. 分别计算 $M_i$ 和 $M_i^{-1}$ ：
   $$\begin{gathered} M_1=5\times 7=35,M_1^{-1}=2^{-1}(mod3)=2; \\ {M}_2=3\times 7=21,M_2^{-1}=1^{-1}(mod5)=1; \\ {M}_3=3\times 5=15,M_3^{-1}=1^{-1}(mod7)=1; \end{gathered}$$

4. 求解 $x$.
   $$\begin{gathered} x=a_1{M}_1{M}_1^{-1}+a_2{M}_2{M}_2^{-1}+a_3{M}_3{M}_3^{-1} \\ =2\times 35\times 2+3\times 21\times 1+2\times 15\times 1 \\ =233 \\ 因此，x\equiv 23(mod105) \end{gathered}$$
   最小的正整数解为 $23$ 。

注意：

\qquad 这个定理给出了在 $m_1,m_2,\dots m_n$ 两两互质这种特殊情况下的构造解，而并不是一元线性同余方程组的通解。

\qquad 我们可以验证这个定理的正确性：

\qquad 只需验证
对于 ∀ i ∈ { x ∈ N ∗ ∣ 1 ≤ x ≤ n } , 均满足 ∑ j = 1 n a j M j M j − 1 ≡ a i ( m o d    m i ) 对于\forall i \in \{x\in \mathbb{N}^*|1\le x \le n\},均满足\sum\limits_{j=1}^{n}a_jM_jM_j^{-1} \equiv a_i(mod\;m_i) 对于∀i∈{x∈N∗∣1≤x≤n},均满足j=1∑n​aj​Mj​Mj−1​≡ai​(modmi​)
\qquad 计算 $\sum _{j=1}^n{a}_j{M}_j{M}_j^{-1}$ 的值时，可以看作这个式子中 $n$ 项分别取模再相加。对于其中某一项：

1. 当 $j\ne i$ 时，其 $M_j$ 中包含因数 $m_i$，因此该项对 $m_i$ 取模得 $0$.
2. 当 $j=i$ 时，其 $M_j\cdot M_j^{-1}$ 对 $m_i$ 取模得 $1$，故实际为 $a_i$ 对 $m_i$ 取模。

综上所述，
$$\sum _{j=1}^n{a}_j{M}_j{M}_j^{-1}(mod{m}_i)=a_i+0+0+\cdots +0=a_i$$

代码：

#include <iostream>
using namespace std;
typedef long long LL;
LL exgcd(LL a, LL b, LL &x, LL &y) {		//扩展欧几里得算法
	if (!b) {
		x = 1, y = 0;
		return a;
	}
	LL d = exgcd(b, a % b, y, x);
	y -= a / b * x;
	return d;
}
LL M = 1, m1[30], m2[30];
LL a[30], m[30], res;
int main() {
	int n;
	scanf("%d", &n);
	for (int i = 1; i <= n; i++) {
		scanf("%lld%lld", a + i, m + i);
		M *= m[i];
	}
	for (int i = 1; i <= n; i++) {
		m1[i] = M / m[i];
		LL x, y;
		LL d = exgcd(m1[i], m[i], x, y);
		if (d == 1) {
			m2[i] = x < 0 ? x + m[i] : x;
		} else {				//说明不满足两两互质
			puts("error!");
			return 0;
		}
		res += a[i] * m1[i] * m2[i] % M;
	}
	res %= M;
	printf("%lld", res < 0 ? res + M : res);
	return 0;
}

例题

204. 表达整数的奇怪方式 - AcWing题库（略微改动）

\qquad 给定 $2n$ 个整数 $a_1,a_2,\dots ,a_n$ 和 $m_1,m_2,\dots ,m_n$，求一个最小的非负整数 $x$，满足 $\forall i\in [1,n],x\equiv a_i(mod{m}_i)$.

数据范围

\qquad $1\le a_i\le 2^{31}-1,$

\qquad $0\le a_i<m_i$

\qquad $1\le n\le 25$

输入样例

2
7 8			//a1和m1
9 11		//a2和m2

输出样例

31

思路

\qquad 此题没有指明 $a_1,a_2,\dots a_n$ 两两互质，因此不能盲目使用中国剩余定理。

\qquad 采用另一种通解，步骤如下：
$$\{\begin{array}{l}x\equiv a_1(mod{m}_1)(1)\\ x\equiv a_2(mod{m}_2)(2)\\ ⋮\\ x\equiv a_n(mod{m}_n)\end{array}$$
\qquad 先研究(1)(2)式组成的方程组，可以写作：
$$\begin{gathered} x=k_1{m}_1+a_1(k_1\in \mathbb{Z}) \\ =k_2{m}_2+a_2(k_2\in \mathbb{Z}) \end{gathered}$$
\qquad 变形为：
$$k_1{m}_1-k_2{m}_2=a_2-a_1(3)$$
\qquad 根据裴蜀定理，要使 $k_1,k_2$ 有解，即 $(m_1,m_2)|(a_2-a_1)$.

\qquad 记 $d=(m_1,m_2)$ ，利用扩展欧几里得算法可以得到 $k_1{m}_1-k_2{m}_2=d$ 的一组解，记作 $k_1^{\prime },k_2^{\prime }$。将其扩大至成为(3)式的解：
$$\{\begin{array}{l}{k}_{10}=k_1^{\prime }\cdot \frac{a_2-a_1}{d}\\ k_{20}=k_2^{\prime }\cdot \frac{a_2-a_1}{d}\end{array}$$
\qquad 那么所有的解为 $\{\begin{array}{l}{k}_1=k_{10}+\frac{m_2}{d}\cdot p\\ k_2=k_{20}+\frac{m_1}{d}\cdot p\end{array},p\in \mathbb{Z}$，则
$$\begin{gathered} x=k_1{m}_1+a_1 \\ =\left(k_{10}+\frac{m_2}{d}\cdot p\right)m_1+a_1 \\ =m_1{k}_{10}+a_1+p\cdot \frac{m_1{m}_2}{d} \\ =m_1{k}_{10}+a_1+p\cdot [m_1,m_2] \end{gathered}$$
\qquad 中括号表示两个数的最小公倍数。

\qquad 这样来看，我们把两个式子融合成了一个形式相同的式子：
$$\begin{gathered} 令M=[m_1,m_2],A=m_1{k}_{10}+a_1 \\ 则:x=p\cdot M+A \\ 即:x\equiv A(modM) \end{gathered}$$

代码

#include <iostream>
using namespace std;
typedef long long LL;
LL exgcd(LL a, LL b, LL &x, LL &y) {		//扩展欧几里得算法
	if (!b) {
		x = 1, y = 0;
		return a;
	}
	LL d = exgcd(b, a % b, y, x);
	y -= a / b * x;
	return d;
}
int main() {
	int n;
	scanf("%d", &n);

	bool flag = true;		//是否有解
	LL a1, m1;

	scanf("%lld%lld", &a1, &m1);

	for (int i = 1; i < n; i++) {
		LL a2, m2;
		scanf("%lld%lld", &a2, &m2);

		LL k1, k2;
		LL d = exgcd(m1, m2, k1, k2);
		if ((a2 - a1) % d) {
			flag = false;		//无解
			break;
		}

		k1 *= (a2 - a1) / d;
		LL t = m2 / d;
		k1 = (k1 % t + t) % t;		//更新k1为符合要求的最小正整数

		a1 = m1 * k1 + a1;			//转化成一个新的式子
		m1 = m1 / d * m2;
	}
	if (flag)
		printf("%lld\n", (a1 % m1 + m1) % m1);
	else
		puts("-1");
	return 0;
}