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
前置芝士：
exBSGSexCRT
一、求出单个同余方程的解
对于 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】快乐肥宅
扫一扫
316
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
