【题解】Luogu-P5345 【XR-1】快乐肥宅

P5345 【XR-1】快乐肥宅

Preface

好题！！！

当然也很毒瘤。

突然感到屠龙勇士在这题面前

就是逊内！！！

Description

给定高次同余方程组
{ k 1 x ≡ r 1 ( m o d g 1 ) k 2 x ≡ r 2 ( m o d g 2 ) ⋯ k n x ≡ r n ( m o d g n ) \begin{cases} k_1^x\equiv r_1\pmod {g_1}\\ k_2^x\equiv r_2\pmod {g_2}\\ \cdots\\ k_n^x\equiv r_n\pmod {g_n} \end{cases} ⎩⎪⎪⎪⎨⎪⎪⎪⎧k1x≡r1(modg1)k2x≡r2(modg2)⋯knx≡rn(modgn)
请求出该方程组的最小非负整数解，若无解或 最小解 ≥ 1 0 9 \ge 10^9 ≥109，输出 Impossible。

对于 100 % 100\% 100% 的数据， 1 ≤ n ≤ 1 0 3 1 \le n \le 10^3 1≤n≤103， 1 ≤ k i , r i ≤ g i ≤ 1 0 7 1 \le k_i, r_i \le g_i \le 10^7 1≤ki,ri≤gi≤107。

Solution

前置芝士：

exBSGS
exCRT

一、求出单个同余方程的解

对于 a x m o d p a^x\bmod p axmodp，结果应该分为一段“尾巴”和一段循环部分，例如 12 4 x m o d 495616 124^x\bmod 495616 124xmod495616：

感谢兔队的图

这时对于一个正整数 b b b， a x ≡ b ( m o d p ) a^x\equiv b\pmod p ax≡b(modp) 的解有 3 3 3 种情况：

唯一解：即 x x x 在尾巴上。例如 b = 15376 b=15376 b=15376 时， x = 2 x=2 x=2。
无穷解：即 x x x 在循环上。例如 b = 241664 b=241664 b=241664 时， x ≡ 3 + 5 ( m o d 5 ) x\equiv 3+5\pmod 5 x≡3+5(mod5) 且 x ≥ 3 + 5 x\ge 3+5 x≥3+5，其中 3 3 3 表示 b b b 是循环中的第 3 3 3 个， + 5 +5 +5 表示尾巴的长度，模数为 5 5 5 表示循环的长度，不等式是为了保证 x x x 在循环上。
无解：例如 b = 114514 b=114514 b=114514 时。

在 exBSGS 中，除到 gcd ⁡ ( a , p ) = 1 \gcd(a,p)=1 gcd(a,p)=1 时，先用 BSGS 求出 a x ≡ b ( m o d p ) a^x\equiv b\pmod p ax≡b(modp) 的最小非负整数解，再求出 a a a 模 p p p 的阶（ gcd ⁡ ( a , p ) = 1 \gcd(a,p)=1 gcd(a,p)=1 时阶一定存在），那么阶就是循环长度。

下面的代码返回一个 pair ⁡ \operatorname{pair} pair，第一个数表示最小解（如果无解则为 − 1 -1 −1），第二个数表示阶（如果不存在则为 − 1 -1 −1）。

与 exBSGS 模板一样，需要特判 b = 1 b=1 b=1 或 p = 1 p=1 p=1 的情况。

pair<int, int> exBSGS(int a, int b, int p)
{if (p == 1){return make_pair(0, 1);}a %= p, b %= p;int g = exgcd(a, p), f = 1, k = 0;while (g > 1){if (b % g){if (b == 1){return make_pair(0, -1);}return make_pair(-1, -1);}k++;b /= g, p /= g;f = (ll)f * (a / g) % p;g = exgcd(a, p);if (f == b){if (g > 1){return make_pair(k, -1);}int y = BSGS(a, 1, p);return make_pair(k, y);}}int x = BSGS(a, (ll)b * inv(f, p) % p, p);if (x == -1){return make_pair(-1, -1);}int y = BSGS(a, 1, p);return make_pair(x % y + k, y);
}

在主函数中：

pair<int, int> res = exBSGS(k, r, g);
b[i] = res.first, a[i] = res.second;
if (b[i] == -1) // x 无解就肯定无解
{puts("Impossible");return 0;
}
else if (a[i] == -1) // 阶不存在，用一个 X 记录唯一的解
{X = b[i];continue;
}
mn = max(mn, b[i]); // x ≥ mn

二、合并解

如果出现了唯一解（即上面代码中的 X ≠ − 1 X\ne -1 X=−1）的情况，直接将这个 X X X 带入剩下的方程中检验即可。

否则直接用 exCRT 合并即可，注意要加上 lcm ⁡ ( a 1 , a 2 , … , a n ) \operatorname{lcm}(a_1,a_2,\dots,a_n) lcm(a1,a2,…,an) 的倍数使得答案不小于上面代码中的 m n mn mn。

我们设前 ( i − 1 ) (i-1) (i−1) 个方程组合并后的解为 x ≡ B ( m o d A ) x\equiv B\pmod A x≡B(modA)，第 i i i 个方程为 x ≡ b ( m o d a ) x\equiv b\pmod a x≡b(moda)。

问题来了， 1 0 3 10^3 103 个 1 0 7 10^7 107，合并后 A , B A,B A,B 都有可能爆 long long。

观察题目中的要求：

若无解或 最小解 ≥ 1 0 9 \ge 10^9 ≥109，输出 Impossible。

当 A > 1 0 9 A>10^9 A>109 时，直接检测 B B B 是否满足 B ≡ b ( m o d a ) B\equiv b\pmod a B≡b(moda) 即可，如果不满足那么一定要加上若干个 A A A，那么一定超过了 1 0 9 10^9 109，不满足要求。

全部合并完之后记得也要判断 B B B 是否在 1 0 9 10^9 109 内。

Code

//18 = 9 + 9 = 18.
#include <iostream>
#include <cstdio>
#include <map>
#include <cmath>
#define Debug(x) cout << #x << "=" << x << endl
typedef long long ll;
using namespace std;ll x, y;ll exgcd(ll a, ll b)
{if (!b){x = 1, y = 0;return a;}ll Gcd = exgcd(b, a % b);ll tmp = x;x = y;y = tmp - a / b * y;return Gcd;
}ll A, B;int merge(ll a, ll b)
{ll Gcd = exgcd(A, a), c = B - b;if (c % Gcd){return -1;}c /= Gcd;x = (x * c % a + a) % a;ll Lcm = A / Gcd * a;B = (B - A * x % Lcm + Lcm) % Lcm;A = Lcm;return 1;
}int div_ceil(int a, int b)
{return (a - 1) / b + 1;
}const int MAXN = 1e3 + 5;
const int D = 1e9;int a[MAXN], b[MAXN];int exCRT(int n, int mn)
{A = a[1], B = b[1];for (int i = 2; i <= n; i++){if (A > D && (B - b[i]) % a[i]){return -1;}else if (A <= D && merge(a[i], b[i]) == -1){return -1;}}if (B < mn){B += div_ceil(mn - B, A) * A;}if (B > D){return -1;}return B;
}int inv(int a, int p)
{exgcd(a, p);x = (x % p + p) % p;return x;
}int phi(int p)
{int ans = p;for (int i = 2; i * i <= p; i++){if (p % i == 0){ans = ans / i * (i - 1);while (p % i == 0){p /= i;}}}if (p != 1){ans = ans / p * (p - 1);}return ans;
}int BSGS(int a, int b, int p)
{map<int, int> hash;int t = ceil(sqrt(phi(p))), aj = 1;for (int j = 0; j < t; j++){int val = (ll)b * aj % p;hash[val] = j;aj = (ll)aj * a % p;}a = aj;int ai = 1;for (int i = 1; i <= t; i++){ai = (ll)ai * a % p;if (hash.find(ai) != hash.end()){int j = hash[ai];if (i * t - j >= 0){return i * t - j;}}}return -1;
}pair<int, int> exBSGS(int a, int b, int p)
{if (p == 1){return make_pair(0, 1);}a %= p, b %= p;int g = exgcd(a, p), f = 1, k = 0;while (g > 1){if (b % g){if (b == 1){return make_pair(0, -1);}return make_pair(-1, -1);}k++;b /= g, p /= g;f = (ll)f * (a / g) % p;g = exgcd(a, p);if (f == b){if (g > 1){return make_pair(k, -1);}int y = BSGS(a, 1, p);return make_pair(k, y);}}int x = BSGS(a, (ll)b * inv(f, p) % p, p);if (x == -1){return make_pair(-1, -1);}int y = BSGS(a, 1, p);return make_pair(x % y + k, y);
}int main()
{int n;scanf("%d", &n);int mn = 0, X = -1;for (int i = 1; i <= n; i++){int k, g, r;scanf("%d%d%d", &k, &g, &r);pair<int, int> res = exBSGS(k, r, g);b[i] = res.first, a[i] = res.second;if (b[i] == -1){puts("Impossible");return 0;}else if (a[i] == -1){X = b[i];continue;}mn = max(mn, b[i]);}if (X != -1){for (int i = 1; i <= n; i++){if ((a[i] == -1 && X != b[i]) || (a[i] != -1 && (X < b[i] || (X - b[i]) % a[i]))){puts("Impossible");return 0;}}printf("%d\n", X);return 0;}int ans = exCRT(n, mn);if (ans == -1){puts("Impossible");}else{printf("%d\n", ans);}return 0;
}

References

[1] 小粉兔：洛谷 P5345: 【XR-1】快乐肥宅