数论再谈
最大公约数 和 最小公倍数
记 \(a\) 和 \(b\) 最大公约数为 \(\gcd(a, b)\), 记 \(a\) 和 \(b\) 最小公倍数数为 \(lcm(a, b)\)。
设 \(a > b\)。
求最大公约数
\(gcd(a, b) = gcd(b, a \pmod b)\)。
证明
如果 \(a\) 是 \(b\) 的倍数， 那么 \(\gcd(a,b) = b\)，如果 \(b = 0\)， 那么 \(\gcd(a, b) = a\)。
下面讨论 \(a\) 不是 \(b\) 的倍数的情况。
设 \(a = bk + r\)，那么 \(a \pmod b = r\)，就是要证明 \(\gcd(a, b) = \gcd(b, r)\)。
记 \(\gcd(a, b) = x\)， 于是就有 \(a | x\)， \(b | x\)。
把上面的式子变一下，就是 \(r = a - bk\)，同时除以一个 \(x\)，于是 \(\frac{r}{x} = \frac{a}{x} - \frac{bk}{x}\)。
\(\because a | x\), \(b | x\)
\(\therefore a | x, bk | x\)
\(\therefore \frac{a}{x},\frac{bk}{x} \in \mathbb{Z}\)
\(\therefore \frac{a}{x} - \frac{bk}{x} \in \mathbb{Z}\)
及 \(\frac{r}{x} \in \mathbb{Z}\)
\(\therefore r | x\)。
那么 \(a, b\) 的约数就是 \(b, r\) 的约数，就可以得到 \(\gcd(a, b) = \gcd(b, r) = \gcd(b, a \pmod b)\)。
得证。
求最小公倍数
\(lcm(a, b) = a \div gcd(a, b) \times b\)。
证明
由唯一分解可以得到
\(a = p_1^{a_1} \times p_2^{a_2} \times \dots \times p_n^{a_n}\)
\(b = p_1^{b_1} \times p_2^{b_2} \times \dots \times p_n^{b_n}\)。
那么最大公约数就是
\(\gcd(a, b) = p_1^{\min(a_1, b_1)} \times p_2^{\min(a_2, b_2)} \times \dots \times p_n^{\min(a_n, b_n)}\)
最小公倍数就是
\(lcm(a, b) = p_1^{\max(a_1, b_1)} \times p_2^{\max(a_2, b_2)} \times \dots \times p_n^{\max(a_n, b_n)}\)
显然 \(\gcd(a, b) \times lcm(a, b) = a \times b\)。
因此就有 \(lcm(a, b) = a \times b \div \gcd(a, b)\)。
然后还可以整一下，设 \(x = \gcd(a, b)\)
于是 \(a = x \times a'，b = x \times b'\)，于是就有 \(gcd(a', b') = 1\)
那么 \(lcm(a, b) = a' \times x \times b' \times x \div x = a' \times x \times b'\)，可以发现这显然是最优的。
扩展欧几里得算法
这个算法可以求解形如这样的二元一次方程 \(ax + by = gcd(a, b)\)， 其中 \(a，b\in \mathbb{Z}\)
求特解
不难发现有无数组解，我们这样设：
\(ax_1 + by_1 = gcd(a, b) \\ bx_2 + (a\pmod b)y_2 = gcd(a, b \pmod a) \)
由前面的欧几里得算法可以得到 \(gcd(a, b) = gcd(a, b \pmod a)\)， 于是就有 \(ax_1 + by_1 = bx_2 + (a\mod b)y_2\)。
然后又有 \(a\pmod b = a - (\lfloor \frac{a}{b} \rfloor \times b)\)。
于是就有 \(ax_1 + by_1 = bx_2 + (a - (\lfloor \frac{a}{b} \rfloor \times b))y_2\)
于是 \(ax_1 + by_1 = ay_2 + b(x_2 - (\lfloor \frac{a}{b} \rfloor) \times y_2)\)
因为 \(a = a, b = b\)， 所以 \(x_1 = y_2, y_1 = x_2 - (\lfloor \frac{a}{b} \rfloor) \times y_2\)。
可以发现上面设的分别是欧几里得算法前后两种状态，所以在求欧几里得的时候就可以顺便把这组解求了。
求最小解
先给出一个结论:
\(x = x_1 + \frac{b}{\gcd(a, b)}\)
\(y = y_1 - \frac{a}{\gcd(a, b)}\)
可以显然的发现带入式子依然成立，因为加了一个 \(\frac{ab}{\gcd(a, b)}\) 又减了一个 \(\frac{ab}{\gcd(a, b)}\)。
于是就可以得到 \(x\) 的值是有周期性的，然后求最小值就是 \(x1 \mod \frac{b}{\gcd(a, b)}\)。
当 c 是任意的数时
因为扩展欧几里得算法求的是 \(ax + by = \gcd(a, b)\)，现在我们要求 \(ax + by = c\) 的不定方程如何求呢？
显然我们让整个式子乘一个 \(c\) 再除一个 \(\gcd(a,b)\) 就可以了，如果 要求解为整数，那么就要 \(c | \gcd(a, b)\) 来保证整个变换中没有小数的出现。
我们求了一组 \(ax + by = \gcd(a, b)\) 的特解后，左右同时子乘一个 \(c\) 再除一个 \(\gcd(a,b)\) 式子依然成立，于是就有
\(a \times \frac{c}{\gcd(a, b)} \times x + b \times \frac{c}{\gcd(a, b)} \times y = c\)
于是解就是 \(\frac{c}{\gcd(a, b)} \times x\)。
求解模方程
模方程也叫同余方程，是形如这样的方程 \(ax \equiv b\pmod c\)。
求解同余方程
如何求解？我们把这个式子变换一下得到：
\(ax - b = cy\)
可以显然推到，然后移相得到 \(ax - cy = b\)，就变成了一个不定方程，所以我们就可以用扩欧解决。
这里 \(c\) 可以取一个负号就是形如这样的式子了 \(ax + cy = b\)。
逆元
随便说说，如果存在一个数 \(x\) 满足 \(ax\equiv 1 \pmod P\)，我们就称 \(x\) 使 \(a\) 关于 \(P\) 的乘法逆元，记作 \(a^{-1}\)，其实就是为了让模运算中存在除法而创造的一个东西。
然后求法这里提供 2 种
exgcd 求逆元
显然就是求 \(ax + Py = 1\)，满足有解的必要条件是 \(\gcd{a, P} = 1\)，我们就有了一种可以判断是否有解的方法。
直接求就完了
int inv(int a, int b) {
int x, y;
exgcd(a, b, x, y);
return (x + b) % b;
}
费马小定理求逆元
费马小定理是如果 \(P\) 是质数就有
于是就有
所以直接用快速幂求逆元就完了。
int inv(int a, int b) {
return power(a, b - 2);
}
Problems
「NOIP2012」同余方程
题意很清楚了，我们按照上面的东西化简一下就可以得到 \(ax + by = 1\)，
然后题目说一定有解，所以直接上 exgcd 模板就可以了。
i64 a, b;
std::cin >> a >> b;
i64 x, y;
exgcd(a, b, x, y);
std::cout << (x + b) % b << std::endl;
「POJ1061」青蛙的约会
题意很清楚了，求这个东西的最小解 \(x + nk \equiv y + mk \pmod L\)
依旧是推式子，可以得到 \((n - m) k + ly = y - x\)，然后大力带模板就完了.
// k(m - n) - lz = -x + y;
// k(n - m) + lz = x - y;
i64 G = gcd(n - m, l);
if (std::abs(x - y) % G) {
std::cout << "Impossible\n";
return 0;
}
i64 mod = std::abs(l / G);
i64 x1, yy;
exgcd(n - m, l, x1, yy);
x1 = (1ll * x1 * (x - y) / G % mod + mod) % mod;
// x1 = x1 * (x - y) / G;
std::cout << x1 << std::endl;
「POJ2115」C Looooops
求这样的模方程 \(a + cx\equiv b(\pmod2^k)\)
于是就有 \(cx + 2^ky = b - a\)。
带 exgcd 进去就完了
i64 a, b, c, k;
while (std::cin >> a >> b >> c >> k) {
if (!a && !b && !c && !k) return 0;
i64 A = c, B = (1ll << k), C = b - a;
if (A == 0) {
if (std::abs(C) % B) {
std::cout << "FOREVER\n";
continue;
} else {
std::cout << "0\n";
continue;
}
}
i64 x1, yy;
exgcd(A, B, x1, yy);
if (C % gcd(A, B)) {
std::cout << "FOREVER" << "\n";
continue;
}
i64 mod = std::abs(B / gcd(A, B));
i64 ans = (x1 * C / gcd(A, B) % mod + mod) % mod;
std::cout << ans << '\n';
}
「NOI2002」荒岛野人
找到最小的 \(m\) 使这个式子对于任意的 \(i\)，\(j\) 成立
于是就有
然后暴力求 \(m\)，用 exgcd 判断一下就可以了。
其中最小解要满足在这些人寿命之外才会相遇。
int mx = -1;
std::vector<int> c(n), p(n), l(n);
for (int i = 0; i < n; i++) {
std::cin >> c[i] >> p[i] >> l[i];
mx = std::max(mx, c[i]);
}
auto check = [&](int m) {
for (int i = 0; i < n; i++) {
for (int j = i + 1; j < n; j++) {
int A = p[i] - p[j], B = m, C = c[j] - c[i];
int D = gcd(A, B);
if (C % D) {
continue;
}
int x, y;
exgcd(A, B, x, y);
int md = std::abs(B / D);
x = (x * C / D % md + md) % md;
if (x <= l[i] && x <= l[j]) {
return false;
}
}
}
return true;
};
for (; ; mx++) {
if (check(mx)) {
std::cout << mx << std::endl;
return 0;
}
}
BSGS 和 exBSGS
BSGS
BSGS 算法可以在 \(O(\sqrt{p})\) 的时间复杂度内求得满足 \(a^x \equiv b\pmod P\) 的 \(x\) 的值。其中满足 \(a\perp P\)，\(\perp\) 表示互质，你知道就当我没说。
我们设 \(x = i\times \lceil \sqrt{P} \rceil - j\)，于是就可以得到 \(a^{i\times \lceil \sqrt{P} \rceil} \equiv b\times a^j \pmod P\)。
我们暴力求 \(b\times a^j\) 的所有情况，用 map 存，然后再枚举左边那个式子就可以了。
i64 m = std::ceil(std::sqrt(P));
a %= P, b %= P;
std::map<i64, i64> mp, vis;
for (int i = 0; i <= m; i++) {
i64 w = b * power(a, i) % P;
mp[w] = i;
vis[w] = 1;
}
for (int i = 0; i <= m; i++) {
i64 w = power(a, i * m) % P;
if (vis[w]) {
i64 x = mp[w];
if ((i * m - x) % P >= 0 && x >= 0) {
std::cout << (i * m - x) % P << std::endl;
return 0;
}
}
}
std::cout << "no solution" << std::endl;
为什么要求 \(a\perp P\)
如果没有这个条件，那我们就无法左右两边同时乘一个 \(a^j\)。
[SDOI2013] 随机数生成器
做这个题顺便就把等比数列学了。。。
等比数列就是满足这个条件的数列， \(a_i = a_1\times p^{i - 1}\)。
然后有一个求和公式啊，这里我们简单的推导一下：
同时乘上一个 \(p\)，就有：
上下两个式子相减得到：
额，这就出来了：
然后再反观这个题。
随便写前面几项，找找有没有规律：
显然有规律，直接表示出来就是 \(x_i=x_1\times a^{i-1}+b\times \sum_{j=0}^{i-2}a^j\)，于是就可以用我们的等比数列了。
显然 \(x_i=x_1\times a^{i-1}+b\times \frac{1-a^{i-1}}{1-a}\)，我们不是要求 \(t\) 吗，带进去就完了。
\(t\equiv x_1\times a^{i-1}+b\times \frac{1-a^{i-1}}{1-a}\pmod{p}\)
下面就是愉快的推式子时间：
用 BSGS 求就完了，然后答案再加一个 1。
exBSGS
如果没有互质的条件，我们如何求？
我们可以尝试创造这个条件，然后再用 BSGS 求解。
可以这样，设 \(d=\gcd(a,P)\)，现在我们变换一下原式子：
\(\frac{a}{d}\times a^{x-1} \equiv \frac{b}{d}\pmod{\frac{p}{d}}\)
现在令 \(a'=a^{x-1}，b'=\frac{b}{d},p'=\frac{p}{d}\)。就这样一直算，直到 \(d=1\)，就是 \(a, p\) 互质了，然后再套上 BSGS 就完了，注意这里要记得求 \(\frac{a}{d}\) 的逆元哦。
最后化简成这个形式
\(\frac{a^{cnt}}{d} \times a^{x-cnt} \equiv b' \pmod{p'}\)
只要满足了BSGS的条件套上就能用，再搞一下就是 \(a^{x-cnt} \equiv b' \times \frac{a^{cnt}}{d} \pmod{p'}\)
如果 \(\frac{a^{cnt}}{d} = b'\),那么显然最小的答案就是现在的 \(cnt\)。
using i64 = long long;
void exgcd(i64 a, i64 b, i64 &x, i64 &y) {
if (!b) {
x = 1, y = 0;
return ;
}
exgcd(b, a % b, x, y);
i64 t = x;
x = y;
y = t - a / b * y;
}
i64 inv(i64 a, i64 P) {
i64 x, y;
exgcd(a, P, x, y);
return (x % P + P) % P;
}
i64 gcd(i64 a, i64 b) {
return b ? gcd(b, a % b) : a;
}
i64 P, a, b;
i64 power(i64 a, i64 b) {
i64 ans = 1;
while (b) {
if (b & 1) ans = ans * a % P;
a = a * a % P;
b >>= 1;
}
return ans;
}
i64 BSGS(i64 a, i64 b, i64 P) {
i64 m = std::ceil(std::sqrt(P));
a %= P, b %= P;
std::unordered_map<i64, i64> mp;
for (int i = 0; i <= m; i++) {
i64 w = b * power(a, i) % P;
mp[w] = i;
}
for (int i = 0; i <= m; i++) {
i64 w = power(a, i * m) % P;
if (mp.find(w) != mp.end()) {
i64 x = mp[w];
if (i * m - x >= 0) {
return (i * m - x + P) % P;
}
}
}
return -1ll;
}
i64 exBSGS(int a, int b, int P) {
if (b == 1 || P == 1) {
return 0;
}
int d, t = 1, cnt = 0;
while ((d = gcd(a, P)) ^ 1) {
if (b % d) {
return -1;
}
P /= d, b /= d; cnt++;
t = (1ll * t * (a / d)) % P;
if (t == b) {
return cnt;
}
}
i64 ans = BSGS(a, b * inv(t, P) % P, P);
if (ans == -1) return -1ll;
return ans + cnt;
}
中国剩余定理和扩展中国剩余定理
中国剩余定理
求满足条件的最小的 \(x\)。其中 \(b\) 两两互质。
先给出一个结论，令 \(B = \prod^n_{i} b_i\)，答案就是 \(\sum^{n}_i a_i \times \frac{B}{b_i} \times inv(\frac{B}{b_i})\)。
\(inv\) 表示逆元。
证明
老祖宗的智慧。
显然 \(\forall i, j\in{[1,n]},i\neq j\) 都满足 \(\gcd(b_i, b_j)=1\)于是有 \(\gcd(b_i, \frac{B}{b_i})=1\)
于是就存在 \(\frac{B}{b_i}\) 关于 \(b_i\) 的乘法逆元， 记作 \(t_i\)，记 \(\frac{B}{b_i}\) 为 \(B_i\)。
于是就可以从 \(a_i \equiv a_i \pmod {b_i}\) 推到 \(a_i \times B_i \times t_i \equiv a_i \pmod {b_i}\)。
因为 \(m\) 相互互质，于是 \(\forall i, j\in{[1,n]},i\neq j\) ，有 \(a_i\times B_i\times t_i \equiv 0 \pmod{m_j}\)，
于是对于 \(x=\sum^{n}_i a_i \times \frac{B}{b_i} \times inv(\frac{B}{b_i})\) 满足
\(\forall i \in [1, n]，x = a_i \times t_i \times B_i + \sum_{j\neq i}{a_j\times t_j \times B_j}=a^i+\sum_{j\neq i} 0 \equiv a_i \pmod {m_i}\) 。得证。
扩展中国剩余定理
和中国剩余定理没什么鸟关系，，，
没想到了扩展中国剩余定理本质是一个暴力啊。
直接讲了，就是两个两个的合并模方程就完了， 推导一下
显然通过我前面写的某些东西可以得到这样的方程组：
然后写成一个等式就是 \(a_1 + b_1 \times k_1 = a_2 + b_2 \times k_2\) 再表示一下就是
\(b_1\times k_1+(-b_2)\times k_2=a_2-a_1\)。
设 \(d = \gcd(b_1, b_2)\)，上面的式子同时除以一个 \(d\)
\(\frac{b_1}{d}\times k_1 - \frac{b_2}{d}\times k_2= \frac{a_2-a_1}{d}\)
于是又可以化成一个模方程
\(\frac{b_1}{d}\times k_1 \equiv \frac{a_2-a_1}{d} \pmod{\frac{b_2}{d}}\)
然后两边同时乘以个 \(\frac{b_1}{d}\) 的在模 \(\frac{b_2}{d}\) 意义下的逆元，就是
\(k_1 \equiv inv(\frac{b_1}{d}) \times \frac{a_2-a_1}{d} \pmod{\frac{b_2}{d}}\),
接着化成等于的形式就是：（不知道为什么我划化来化去的
\(k_1 = inv(\frac{b_1}{d})\times \frac{a_2-a_1}{d} + \frac{b_2}{d} \times y\)，然后我们再把 \(k_1\) 带进原式中。
\(x = a_1 + b_1 \times (inv(\frac{b_1}{d})\times \frac{a_2-a_1}{d} + \frac{b_2}{d} \times y)\)
\(x = a_1 + b_1\times inv(\frac{b_1}{d})\times \frac{a_2-a_1}{d} + \frac{b_1\times b_2}{d}\times y\)，最后的最后就是
\(x \equiv a_1 + b_1\times inv(\frac{b_1}{d})\times \frac{a_2-a_1}{d} \pmod{\frac{b_1\times b_2}{d}}\)
可以参照之前的一些东西 \(\frac{b_1\times b_2}{\gcd(b_1, b_2)} = lcm(b_1, b_2)\)。
\(x \equiv a_1 + b_1\times inv(\frac{b_1}{d})\times \frac{a_2-a_1}{d} \pmod{lcm(b_1, b_2)}\)
就合并成功了。 用 exgcd 求解就好了。
int n;
i64 a[100010], b[100010];
i64 exCRT() {
i64 a1 = a[0], b1 = b[0];
i64 x, y;
for (int i = 1; i < n; i++) {
__int128 a2 = a[i], b2 = b[i];
__int128 d = gcd(b1, b2), g = a2 - a1;
exgcd(b1, b2, x, y);
if (g % d) {
return -1ll;
}
__int128 md = b2 / d;
x = (__int128)(x * g / d % md + md) % md;
a1 = __int128(x * b1 + a1);
b1 = lcm(b1, b2);
a1 %= b1;
}
return a1 % b1;
}
