[BZOJ2242][SDOI2011]计算器
一道非常经典的同余方程入门题
题意
- 计算pow(a, b)
- 线性同余方程
- 高次同余方程
分析
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 \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)}\)即可。
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;
}
