[BZOJ2242][SDOI2011]计算器

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

题意

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

分析

pow(a, b)

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

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 mod b))。假设存在一对整数 (x, y) 满足 (b * x + (a mod b) * y = gcd(b, a mod b)).
因为 $$bx + (a mod b)y = bx + (a - b * [a / b])y = ay - b(x - [a / b]y)$$
所以令 (mathop{{x}'} = y, mathop{{y}'} = x - [a / b]y) 就得到了 (amathop{{x}'} + bmathop{{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)})即可。

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)，更新答案即可。

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

完整代码

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