【刷题日记】求解同余方程（扩展欧几里得算法）

扩展欧几里得算法

引入

欧几里得定理告诉我们：\(\gcd(a, b) = \gcd(b, a \mod b)\)。

扩展欧几里得算法，常用于求\(ax + by = \gcd(a, b)\)的一组可行解。

当\(b = 0\)时，解为\(x = 1, y = 0\)。

否则：

推理

设
\(ax_1 + by_1 = \gcd(a,b)\)

\(bx_2 + (a \mod b)y_2 = \gcd(b, a \mod b)\)

由欧几里得定理可知：\(\gcd(a,b) = \gcd(b, a \mod b)\)

所以 \(ax_1 + by_1 = bx_2 + (a \mod b)y_2\)

又因为 \(a \mod b = a - \left(\left\lfloor \frac{a}{b} \right\rfloor \times b\right)\)

所以 \(ax_1 + by_1 = bx_2 + \left(a - \left\lfloor \frac{a}{b} \right\rfloor \times b\right)y_2\)

\(ax_1 + by_1 = ay_2 + bx_2 - \left\lfloor \frac{a}{b} \right\rfloor \times by_2 = ay_2 + b\left(x_2 - \left\lfloor \frac{a}{b} \right\rfloor y_2\right)\)

因为 \(a = a,b = b\)，

所以 \(x_1 = y_2, y_1 = x_2 - \left\lfloor \frac{a}{b} \right\rfloor y_2\)

将x2， y2不断代入递归求解，直到回到x = 1, y = 0的基础情况。

代码实现

int Exgcd(int a, int b, int &x, int &y) {
  if (!b) {
    x = 1;
    y = 0;
    return a;
  }
  int d = Exgcd(b, a % b, x, y);
  int t = x;
  x = y;
  y = t - (a / b) * y;
  return d;
}

函数返回值为gcd。x， y为全局变量，运行结束后均已完成修改。

原方程的解有无数个，有的解会爆long long，而扩展欧几里得算法给出的解必有\(|x|<=b\), \(|y|<=a\)。

第一题：同余方程

题目链接

洛谷 P1082

题目描述

求关于 $ x$ 的同余方程 $ a x \equiv 1 \pmod {b}$ 的最小正整数解。

输入格式

一行，包含两个整数 \(a,b\)，用一个空格隔开。

输出格式

一个整数 \(x_0\)，即最小正整数解。输入数据保证一定有解。

输入输出样例

输入

3 10

输出

7

数据规模与约定

• 对于 \(40\%\) 的数据，\(2 ≤b≤ 1,000\)；
• 对于 \(60\%\) 的数据，\(2 ≤b≤ 50,000,000\)；
• 对于 \(100\%\) 的数据，\(2 ≤a, b≤ 2,000,000,000\)。

思路

\[ax \equiv 1 (mod b) \]

得

\[ax + by = 1 \]

其中\(y\)为整数，作为辅助变量，我们要求的是上式中的最小正整数\(x\)。

对于\(ax + by = m\)，这个方程有解的必要条件是 \(m \mod \gcd(a, b) = 0\).证明如下：

因为\(a\)是\(\gcd(a, b)\)的倍数，\(b\)是\(\gcd(a, b)\)的倍数，\(x\)和\(y\)均为整数，所以\(ax + by\)是\(\gcd(a, b)\)的倍数。

在这道题中，\(m = 1\)，则\(1\)是\(\gcd(a, b)\)的倍数才能有解。题目保证有解，可知：\(\gcd(a, b) = 1\)。

那么我们要求解\(ax + by = 1\)，就是求解：

\(ax + by = \gcd(a, b)\)

直接套用扩展欧几里得算法。

代码

#include <bits/stdc++.h>
using namespace std;
using i64 = long long;

i64 x, y;

void exgcd(i64 a, i64 b) {
	if (b == 0) {
		x = 1;
		y = 0;
		return;
	}
	exgcd(b, a % b);
	i64 t = x;
	x = y;
	y = t - a / b * y;
	return;
}

int main() {
	i64 a, b;
	cin >> a >> b;
	exgcd(a, b);
	x = (x % b + b) % b;	//将x处理成有效的最小正整数
	cout << x;
	return 0;
}

第二题：有理数取余

题目链接

洛谷 P2613

题目描述

给出一个有理数 \(c=\frac{a}{b}\)，求 \(c \bmod 19260817\) 的值。

这个值被定义为 \(bx\equiv a\pmod{19260817}\) 的解。

输入格式

一共两行。

第一行，一个整数 \(a\)。
第二行，一个整数 \(b\)。

输出格式

一个整数，代表求余后的结果。如果无解，输出 Angry!。

输入输出样例

输入

233
666

输出

18595654

说明/提示

对于所有数据，保证 \(0\leq a \leq 10^{10001}\)，\(1 \leq b \leq 10^{10001}\)，且 \(a, b\) 不同时是 \(19260817\) 的倍数。

思路

求\(\frac{a}{b} \mod p\)的值，假设答案为\(x\)。

\[\frac{a}{b} \mod p = x \]

两边同乘以\(b\)得：

\[bx \equiv a (\mod p) \]

不妨先计算

\[bx_1 \equiv 1 (\mod p) \]

那么，\(x = (x_1 * a) \mod p\)。

那么$bx_1 \equiv 1 (\mod p) $怎么计算呢？与第一道题完全一致。

现在思考在什么时候会无解？

当\(b\)是\(p\)的倍数时，\(a\)不为\(p\)的倍数（这是题目中说的，两者不能全为\(p\)的倍数）。这时不可能有解。

小技巧：用快读处理大数，边读边取模。在这道题中 $ bx \equiv a (\mod p) $ 等价于 $ (b \mod p)x \equiv (a \mod p) (\mod p) $

代码

#include <bits/stdc++.h>
using namespace std;
using i64 = long long;

const i64 p = 19260817;
i64 x, y;

i64 read() {
	char c = getchar();
	i64 res = 0;
	while (!isdigit(c) and c != EOF)
		c = getchar();
	while (isdigit(c)) {
		res = (res << 3) + (res << 1) + (c - '0');
		res %= p;
		c = getchar();
	}
	return res;
}

void exgcd(i64 a, i64 b) {
	if (b == 0) {
		x = 1;
		y = 0;
		return;
	}
	exgcd(b, a % b);
	i64 t = x;
	x = y;
	y = t - a / b * y;
	return;
}

int main() {
	i64 a, b;
	a = read();
	b = read();
	if (b == 0) {
		cout << "Angry!";
		return 0;
	}
	exgcd(b, p);
	x = (x % p + p) % p;	//将x1处理成有效的最小正整数
	x = (x * a) % p;	//将x1转化为x
	cout << x;
	return 0;
}