[BZOJ2242][SDOI2011]计算器

一道非常经典的同余方程入门题

题意

1. 计算pow(a, b)
2. 线性同余方程
3. 高次同余方程

分析

pow(a, b)#

我们可以使用快速幂完成这项操作

Copy
qword qpow(qword a, qword b, qword p) {
	a %= p; qword res = 1 % p; for (; b; b >>= 1, a = a * a % p) if (b & 1) res = res * a % p; return res;
}

线性同余方程#

求解线性同余方程我们一般使用拓展欧几里得算法。
这里稍微证明一下这种算法的正确性并描述一下这种算法

Bezout(贝祖)定理：对于任意整数 $a, b$，存在一对整数 $x, y$，满足 $ax + by = gcd(a, b)$
证明：
当 $b = 0$时， 显然有一对整数 $x = 1, y = 0$ 使得 $a * 1 + 0 * 0 = gcd(a, 0)$
若 $b > 0$，则 $gcd(a, b) = gcd(b, a \bmod b)$。假设存在一对整数 $x, y$ 满足 $b * x + (a \bmod b) * y = gcd(b, a \bmod b)$.
因为 $$bx + (a \bmod b)y = bx + (a - b * [a / b])y = ay - b(x - [a / b]y)$$
所以令 $\mathop{{x}'} = y, \mathop{{y}'} = x - [a / b]y$ 就得到了 $a\mathop{{x}'} + b\mathop{{y}'} = gcd(a, b)$
应用数学归纳法，可知定理成立。

Bezout定理是按照欧几里得算法的流程进行证明的，所以这种能同时计算 $x, y$ 的算法叫做扩展欧几里得算法
例如在本题中，由于 $y, z$ 是给定的参数，我们不妨将其设为 $a, b$
原式变成了 $ax = b + py$ 我们令 $p = -p$ 就有 $ax + py = b$。（这里的 $y$ 是同余方程中设出的另一个解）
易知，当 $gcd(a, p) | b$ 时有解。我们最终只需要解出 $ax + py = gcd(a, p)$ 的值，并且扩大 $\frac{b}{gcd(a,p)}$即可。

Copy
qword exgcd(qword a, qword b, qword &x, qword &y) {
	if (b == 0) { x = 1, y = 0; return a; }
	qword d = exgcd(b, a % b, x, y);
	qword z = x; x = y; y = z - y * (a / b);
	return d;
}

高次同余方程#

求解高次同余方程我们一般使用Baby Step Gaint Step算法
求解形式：$a^x = b \pmod p$, 要求 $a, p$互质
算法复杂度：$O(\sqrt{p})
因为 $a, p$ 互质，所以可以在模 $p$ 意义下执行关于 $a$ 的乘、除运算。
设 $x = i * t - j$ 其中 $t = [\sqrt{p}], 0 \le j \le t - 1$，则方程变为 $a ^{i * t - j} = b \pmod p$,即 $(a^t)^i = b * a ^ j \pmod p$
对于所有的 $j \in [0, t - 1]$ ，把 $b * a^j \pmod p$ 插入一个hash表中
枚举 $i$ 的所有可能取值，计算出 $(a^t)^i$ 在hash表中是否存在对应的 $j$，更新答案即可。

Copy
qword baby_step_gaint_step(qword a, qword b, qword p) {
	mp.clear();
	b %= p;
	qword t = (qword)sqrt(p) + 1;
	for (int j = 0; j < t; ++ j) {
		qword val = (qword)b * qpow(a, j, p) % p;
		mp[val] = j;
	}
	a = qpow(a, t, p);
	if (a == 0) return b == 0 ? 1 : -1;
	for (int i = 0; i <= t; ++ i) {
		qword val = qpow(a, i, p); 
		qword j = mp.find(val) == mp.end() ? -1 : mp[val];
		if (j >= 0 && i * t - j >= 0) return i * t - j;
	}
	return -1;
}

完整代码

Copy
#include <bits/stdc++.h>
using namespace std;

typedef long long qword;

qword qpow(qword a, qword b, qword p) {
	a %= p; qword res = 1 % p; for (; b; b >>= 1, a = a * a % p) if (b & 1) res = res * a % p; return res;
}

qword exgcd(qword a, qword b, qword &x, qword &y) {
	if (b == 0) { x = 1, y = 0; return a; }
	qword d = exgcd(b, a % b, x, y);
	qword z = x; x = y; y = z - y * (a / b);
	return d;
}

map<qword, qword> mp;

qword baby_step_gaint_step(qword a, qword b, qword p) {
	mp.clear();
	b %= p;
	qword t = (qword)sqrt(p) + 1;
	for (int j = 0; j < t; ++ j) {
		qword val = (qword)b * qpow(a, j, p) % p;
		mp[val] = j;
	}
	a = qpow(a, t, p);
	if (a == 0) return b == 0 ? 1 : -1;
	for (int i = 0; i <= t; ++ i) {
		qword val = qpow(a, i, p); 
		qword j = mp.find(val) == mp.end() ? -1 : mp[val];
		if (j >= 0 && i * t - j >= 0) return i * t - j;
	}
	return -1;
}

qword t, op, x, y;

int main() {
	cin >> t >> op;
	for (int i = 1, a, b, p; i <= t; ++ i) {
		cin >> a >> b >> p;
		if (op == 1) cout << qpow(a, b, p) << endl;
		else if (op == 2) {
			qword d = exgcd(a, p, x, y);
			if (b % d) cout << "Orz, I cannot find x!" << endl;
			else {
				x *= b / d;
				cout << (x % p + p) % p << endl;
			}
		} else if (op == 3){
			x = baby_step_gaint_step(a, b, p);
			if (x == -1) cout << "Orz, I cannot find x!" << endl;
			else cout << (x % p + p) % p << endl;
		}
	}
	return 0;
}