題目傳送門：洛谷 P5345。
很榮幸為 X Round 1 貢獻了自己的一題。
題意簡述：
給定 \(n\) 組 \(k_i,g_i,r_i\)（\(0\le k_i,r_i<g_i\le 10^7\)）。
求關於 \(x\) 的方程組 \(\left\{\begin{matrix}k_1^x\equiv r_1\pmod{g_1}\\k_2^x\equiv r_2\pmod{g_2}\\\vdots\\k_n^x\equiv r_n\pmod{g_n}\end{matrix}\right.\)（定義 \(0^0=1\)）在 \([0,10^9]\) 內的最小整數解，或判斷在這個范圍內無解。
題解：
分為兩個部分解決。
第一部分：分別求出每個方程的解。
首先觀察 \(k^j\bmod g\) 的規律，以 \(g=495616\)，\(k=124\) 為例：
定義 \(f(i)=ki\bmod g\)，則把每個在 \([0,g)\) 之間的整數看作一個節點后，形成了一個基環內向森林。
從 \(1\bmod g\) 開始走唯一的出邊，必然形成一個 \(\rho\) 形結構，讓我們模擬這個過程：
可以看到，出現了長度為 \(5\) 的循環。
-
如果 \(r=245760\)，則 \(x\equiv 7\pmod{5}\)，且 \(x\ge 7\)。
-
如果 \(r=12544\)，則 \(x=4\)。
-
如果 \(r=2\)，則無解。
以上就是每個方程的 \(3\) 種解，無限解，唯一解或無解。
那么，問題就是，如何區分這 \(3\) 種解，以及如何求出解。
首先，需要將 \(\rho\) 形的“尾巴”和循環分開考慮，不難證明“尾巴”的長度不超過 \(\log_2g\)。
而且可以發現，若在“尾巴”上找到了解，那么就只有一組解。所以需要判斷一個值是否在尾巴上。
其實這是很簡單的，如果 \(\gcd(x,m)\neq\gcd(f(x),m)\)，那么 \(x\) 就在尾巴上，反之就在循環內。
據此，我們區分了在“尾巴”上的解，也就是只有一組解的情況，接下來考慮不在“尾巴”上的情況。
因為解不在“尾巴”上，則如果解不在循環內就是無解了，所以首先判斷解是否在循環內。
假設第一個在循環上的數為 \(c\)，尾巴長度為 \(o\)，則有 \(c\equiv k^o\pmod{g}\)。
令 \(d=\gcd(c,g)\)，則原方程 \(k^x\equiv r\pmod{g}\) 可化為 \(c\cdot k^{x-o}\equiv r\pmod{g}\) 且 \(x\ge o\)，這是因為我們假定解一定在循環內。
可以進一步化為 \(k^{x-o}\equiv\frac{r}{d}\left(\frac{c}{d}\right)^{-1}\pmod{\frac{g}{d}}\)，若 \(r\) 不是 \(d\) 的倍數則無解。
這是因為循環內的值均是 \(d\) 的倍數，可以讓方程同除 \(d\)，又因為 \(\frac{c}{d}\perp\frac{g}{d}\)，可以取逆元。
令 \(a=k,b=\frac{r}{d}\left(\frac{c}{d}\right)^{-1},m=\frac{g}{d}\)，則有 \(a\perp m\)。
使用 BSGS 求出 \(a^y\equiv b\pmod{m}\) 的最小自然數解 \(y\) 后，有原方程的最小自然數解為 \(y+o\)。
再使用 BSGS 求出 \(a^z\equiv 1\pmod{m}\) 的最小正整數解 \(z\)，即 \(a\) 對 \(m\) 的階。
則原方程的解為 \(x\equiv y+o\pmod{z}\) 且 \(x\ge y+o\)。
至此，第一部分解決。
第二部分：合並每個方程的解。
首先，若前面的方程出現了無解的情況，則方程組也無解。
若前面的方程出現了只有唯一解的情況，只需要檢查此唯一解是否滿足所有方程即可。
接下來討論以上每個方程的解均形如 \(x\equiv r\pmod{q}\) 且 \(x\ge r\) 的形式。
分開考慮前半部分和后半部分，對於前半部分，顯然是 ExCRT 的形式。
對於后半部分，可以化為 \(x\ge\max r_i\)，令 \(\max r_i=x_0\)，放到最后考慮。
考慮使用 ExCRT 解決前半部分，令前 \(i\) 個方程組的解為 \(x\equiv P_i\pmod{Q_i}\)，有 \(P_0=0,Q_0=1\)。
但是這里出現一個問題，那就是 \(P_i,Q_i\) 可能很大，但是這里其實沒必要使用高精度，要怎么處理這種情況呢？
因為以上方程的解 \(x\equiv r\pmod{q}\) 中，\(q\) 均小於 \(10^7\)，所以不需要擔心這部分。
考慮這樣一種情況：\(Q_{i-1}>10^9\)，而再合並進 \(x\equiv r_i\pmod{q_i}\) 可能會讓 \(Q_i\) 超出 long long 能夠表示的范圍。
這時其實不需要再合並了，只需要判斷 \(P_{i-1}\) 是否滿足第 \(i\) 個方程即可。因為若不滿足，合並后 \(x\) 必然會至少增加一倍的 \(Q_{i-1}\)，這就超出了 \(10^9\) 的范圍，從而可以直接輸出無解。
最后，當方程成功合並完之后，\(Q_n\) 可能不是真實值，但是當 \(Q_n\le 10^9\) 時必然是真實值。
此時可以再考慮 \(x_0\)，真實的解應該為 \(x=\left\lceil\frac{x_0-P_n}{Q_n}\right\rceil\cdot Q_n+P_n\)。
若 \(Q_n> 10^9\)，則必然要滿足 \(P_n\ge x_0\)，否則無解。
據此寫出代碼，復雜度 \(\mathcal{O}\left(\sum_{i=1}^{n}\left(\sqrt{g_i}+\log g_i\right)\right)\)：
#include <cstdio>
#include <algorithm>
#include <cstring>
#include <cmath>
#define Fail return puts("Impossible"), 0
#define mp std::make_pair
typedef long long LL;
typedef std::pair<int, int> pii;
const int MG = 10000005, MS = 3175;
const int U = 1e9;
int gcd(int a, int b) { return b ? gcd(b, a % b) : a; }
template<typename T>
T exgcd(T a, T b, T &x, T &y) {
if (!b) return x = 1, y = 0, a;
int d = exgcd(b, a % b, y, x);
return y -= a / b * x, d;
}
int buk[MG], stk[MS];
inline int BSGS(int a, int b, int m) {
int S = sqrt(m - 1) + 1;
int A = 1, f = -1, t = 0;
for (int i = 0; i < S; ++i) {
buk[stk[++t] = (LL)b * A % m] = i;
A = (LL)A * a % m;
}
int C = 1;
for (int i = 1; !~f && i <= S; ++i)
if (~buk[C = (LL)C * A % m])
f = i * S - buk[C];
while (t) buk[stk[t--]] = -1;
return f;
}
inline pii ExBSGS(int a, int b, int m) {
int o = 0, A = 1 % m, d = 1, nd, x, y;
while (1) {
if (d == (nd = gcd((LL)A * a % m, m))) break;
if (A == b) return mp(o, -1);
++o, A = (LL)A * a % m, d = nd;
}
if (b % d) return mp(-1, -1);
m /= d, b /= d, A /= d;
exgcd(A, m, x, y);
b = (LL)b * (x + m) % m;
x = BSGS(a, b, m);
if (!~x) return mp(-1, -1);
y = BSGS(a, 1 % m, m);
return mp(x % y + o, y);
}
inline bool Combine(LL &a1, LL &m1, LL a2, LL m2) {
LL k1, k2, g = exgcd(m1, m2, k1, k2);
if ((a2 - a1) % g) return 0;
a1 += (k1 * ((a2 - a1) / g) % m2 + m2) * m1 % (m1 / g * m2);
return a1 %= m1 *= m2 / g, 1;
}
const int MN = 1005;
int N;
int k[MN], r[MN], g[MN];
int x[MN], q[MN], X, MaxX;
int main() {
memset(buk, -1, sizeof buk), X = -1;
scanf("%d", &N);
for (int i = 1; i <= N; ++i) {
scanf("%d%d%d", &k[i], &g[i], &r[i]);
pii ans = ExBSGS(k[i] % g[i], r[i] % g[i], g[i]);
x[i] = ans.first, q[i] = ans.second;
if (!~x[i]) Fail;
if (!~q[i]) X = x[i];
MaxX = std::max(MaxX, x[i]);
}
if (~X) {
for (int i = 1; i <= N; ++i) {
if (~q[i] && (X < x[i] || (X - x[i]) % q[i])) Fail;
if (!~q[i] && X != x[i]) Fail;
}
return printf("%d\n", X), 0;
}
LL P = 0, Q = 1;
for (int i = 1; i <= N; ++i) {
if (Q > U && (P - x[i]) % q[i]) Fail;
if (Q <= U && !Combine(P, Q, x[i] % q[i], q[i])) Fail;
}
if (P < MaxX) P += ((MaxX - P - 1) / Q + 1) * Q;
if (P > U) Fail;
printf("%lld\n", P);
return 0;
}