i64 qpow(i64 a, i64 b) {
i64 res = 1;
while (b) {
if (b & 1) {
res = res * a;
}
a = a * a;
b >>= 1;
}
return res;
}
i64 qpow(__int128 a, i64 b, i64 c, __int128 res = 1) {
res %= c;
a %= c;
while (b) {
if (b & 1) {
res = res * a % c;
}
a = a * a % c;
b >>= 1;
}
return res;
}
template<class T> struct Random {
mt19937 mt;
Random() : mt(chrono::steady_clock::now().time_since_epoch().count()) {}
T operator()(T L, T R) {
uniform_int_distribution<int64_t> dist(L, R);
return dist(mt);
}
};
Random<i64> rng;
namespace Miller_Rabin {
bool Miller_Rabin(i64 n, const vector<i64> &as) {
i64 d = n - 1;
while (d % 2 == 0) {
d /= 2;
}
i64 e = 1, rev = n - 1;
for (auto &a : as) {
if (n <= a) {
break;
}
i64 t = d;
__int128 y = qpow(a, t, n);
while (t != n - 1 && y != e && y != rev) {
y = y * y % n;
t *= 2;
}
if (y != rev && t % 2 == 0) {
return false;
}
}
return true;
}
bool is_prime(i64 n) {
if (n % 2 == 0) {
return n == 2;
}
if (n <= 1) {
return false;
}
if (n < (1LL << 30)) {
return Miller_Rabin(n, {2, 7, 61});
}
return Miller_Rabin(n, {2, 325, 9375, 28178, 450775, 9780504, 1795265022});
}
} // Miller_Rabin
using Miller_Rabin::is_prime; // 判断质数
namespace Pollard_rho {
template<class T>
struct Random {
mt19937 mt;
Random() : mt(chrono::steady_clock::now().time_since_epoch().count()) {}
T operator()(T L, T R) {
uniform_int_distribution<int64_t> dist(L, R);
return dist(mt);
}
};
Random<i64> rng;
i64 get(i64 n) {
if (n % 2 == 0) {
return 2;
}
if (is_prime(n)) {
return n;
}
i64 R;
auto f = [&](__int128 x) {
return (x * x + R) % n;
};
while (true) {
i64 x, y, ys, q = 1;
R = rng(2, n - 1), y = rng(2, n - 1);
i64 g = 1;
int m = 128;
for (int r = 1 ; g == 1 ; r *= 2) {
x = y;
for (int i = 0 ; i < r ; i++) {
y = f(y);
}
for (int k = 0 ; g == 1 && k < r ; k += m) {
ys = y;
for (int i = 0 ; i < m && i < r - k ; i++) {
y = f(y);
q = (__int128) q * (((x - y) % n + n) % n) % n;
}
g = __gcd(q, n);
}
}
if (g == n) {
do {
ys = f(ys);
g = __gcd(((x - ys) % n + n) % n, n);
} while (g == 1);
}
if (g != n) {
return g;
}
}
return 0;
}
vector<i64> factorize(i64 n) {
if (n <= 1) {
return {};
}
i64 p = get(n);
if (p == n) {
return {n};
}
auto L = factorize(p);
auto R = factorize(n / p);
copy(R.begin(), R.end(), back_inserter(L));
return L;
}
vector<pair<i64, int>> prime_factor(i64 n) {
auto ps = factorize(n);
sort(ps.begin(), ps.end());
vector<pair<i64, int>> res;
for (auto &e : ps) {
if (!res.empty() && res.back().first == e) {
res.back().second++;
} else {
res.emplace_back(e, 1);
}
}
return res;
}
vector<i64> divisors(i64 n) {
auto ps = prime_factor(n);
int cnt = 1;
for (auto &[p, t] : ps) {
cnt *= t + 1;
}
vector<i64> res(cnt, 1);
cnt = 1;
for (auto &[p, t] : ps) {
i64 pw = 1;
for (int i = 1 ; i <= t ; i++) {
pw *= p;
for (int j = 0 ; j < cnt ; j++) {
res[cnt * i + j] = res[j] * pw;
}
}
cnt *= t + 1;
}
return res;
}
} // Pollard_rho
using Pollard_rho::prime_factor; // 所有质因子 (质因子, 个数)
using Pollard_rho::divisors; // 所有因子
namespace Pohlig_Hellman {
i64 BSGS(i64 A, i64 B, i64 C, i64 P) {
int m = ceil(sqrtl(C));
A %= P;
B %= P;
unordered_map<i64, int> vis;
for (int i = 1 ; i <= m ; i++) {
B = B * A % P;
vis[B] = i;
}
i64 temp = qpow(A, m, P);
B = 1;
for (int i = 1 ; i <= m ; i++) {
B = temp * B % P;
if (vis.count(B)) {
return ((__int128) i * m - vis[B] + P) % P;
}
}
return -1;
}
i64 getK(i64 A, i64 B, i64 prime, i64 cnt, i64 phi, i64 P) {
vector<i64> pi;
i64 temp = 1;
for (int i = 0 ; i <= cnt ; i++) {
pi.push_back(temp);
temp *= prime;
}
i64 k = qpow(A, pi[cnt - 1], P);
i64 inv = 0;
temp = 1;
for (int i = cnt - 1 ; i >= 0 ; i--) {
i64 tp = qpow(A, pi[cnt] - inv, P);
i64 tx = (__int128) temp * BSGS(k, qpow((__int128) B * tp % P, pi[i], P), prime, P);
inv += tx;
temp *= prime;
}
return inv;
}
int getMinRoot(i64 P, i64 phi, const vector<pair<i64, int>>& factors ) {
for (int k = 2 ; ; k++) {
bool flag = true;
for (auto &[x, y] : factors) {
if (qpow(k, phi / x, P) == 1LL) {
flag = false;
break;
}
}
if (flag) return k;
}
}
void Exgcd(i64 a, i64 b, i64 &x, i64 &y) {
if (b == 0) {
x = 1;
y = 0;
return;
}
Exgcd(b, a % b, y, x);
y -= a / b * x;
}
i64 CRT(const vector<i64> &a, const vector<pair<i64, int>> &factors) {
// p_i 两两互质
// ans % p_i = a_i
i64 n = 1, ans = 0;
int len = a.size();
vector<i64> p(factors.size());
for (int i = 0 ; i < len ; i++) {
p[i] = qpow(factors[i].first, factors[i].second);
n *= p[i];
}
for (int i = 0 ; i < len ; i++) {
i64 m = n / p[i];
i64 x, y;
Exgcd(m, p[i], x, y);
ans = (ans + (__int128) a[i] * m % n * x % n + n) % n;
}
return ans;
}
i64 getX(i64 A, i64 B, i64 phi, i64 P, const vector<pair<i64, int>> &res) {
vector<i64> k;
for (auto &[x, y] : res) {
i64 z = qpow(x, y);
i64 tA = qpow(A, phi / z, P);
i64 tB = qpow(B, phi / z, P);
k.push_back(getK(tB, tA, x, y, phi, P));
}
return CRT(k, res);
}
i64 solve(i64 A, i64 B, i64 P) {
if (B == 1LL) {
return 0LL;
}
i64 phi = P - 1;
vector<pair<i64, int>> res = Pollard_rho::prime_factor(phi);
int rt = getMinRoot(P, phi, res);
i64 x = getX(A, rt, phi, P, res), y = getX(B, rt, phi, P, res);
i64 a, b;
if (x == 0LL) {
if (y == 0LL) {
return 1LL;
} else if (y == 1LL) {
return 0LL;
}
return -1LL;
}
i64 d = __gcd(x, phi);
if (y % d != 0) return -1;
x /= d;
phi /= d;
y /= d;
Exgcd(x, phi, a, b);
a = ((__int128) a * y % phi + phi) % phi;
return a;
}
} // Pohlig_Hellman