#line 1 "library/Template/template.hpp" #include using namespace std; #define rep(i, a, b) for (int i = (int)(a); i < (int)(b); i++) #define rrep(i, a, b) for (int i = (int)(b)-1; i >= (int)(a); i--) #define ALL(v) (v).begin(), (v).end() #define UNIQUE(v) sort(ALL(v)), (v).erase(unique(ALL(v)), (v).end()) #define SZ(v) (int)v.size() #define MIN(v) *min_element(ALL(v)) #define MAX(v) *max_element(ALL(v)) #define LB(v, x) int(lower_bound(ALL(v), (x)) - (v).begin()) #define UB(v, x) int(upper_bound(ALL(v), (x)) - (v).begin()) using uint = unsigned int; using ll = long long int; using ull = unsigned long long; using i128 = __int128_t; using u128 = __uint128_t; const int inf = 0x3fffffff; const ll INF = 0x1fffffffffffffff; template inline bool chmax(T &a, T b) { if (a < b) { a = b; return 1; } return 0; } template inline bool chmin(T &a, T b) { if (a > b) { a = b; return 1; } return 0; } template T ceil(T x, U y) { assert(y != 0); if (y < 0) x = -x, y = -y; return (x > 0 ? (x + y - 1) / y : x / y); } template T floor(T x, U y) { assert(y != 0); if (y < 0) x = -x, y = -y; return (x > 0 ? x / y : (x - y + 1) / y); } template int popcnt(T x) { return __builtin_popcountll(x); } template int topbit(T x) { return (x == 0 ? -1 : 63 - __builtin_clzll(x)); } template int lowbit(T x) { return (x == 0 ? -1 : __builtin_ctzll(x)); } template ostream &operator<<(ostream &os, const pair &p) { os << "P(" << p.first << ", " << p.second << ")"; return os; } template ostream &operator<<(ostream &os, const vector &vec) { os << "{"; for (int i = 0; i < vec.size(); i++) { os << vec[i] << (i + 1 == vec.size() ? "" : ", "); } os << "}"; return os; } template ostream &operator<<(ostream &os, const map &map_var) { os << "{"; for (auto itr = map_var.begin(); itr != map_var.end(); itr++) { os << "(" << itr->first << ", " << itr->second << ")"; itr++; if (itr != map_var.end()) os << ", "; itr--; } os << "}"; return os; } template ostream &operator<<(ostream &os, const set &set_var) { os << "{"; for (auto itr = set_var.begin(); itr != set_var.end(); itr++) { os << *itr; ++itr; if (itr != set_var.end()) os << ", "; itr--; } os << "}"; return os; } #ifdef LOCAL #define show(...) _show(0, #__VA_ARGS__, __VA_ARGS__) #else #define show(...) true #endif template void _show(int i, T name) { cerr << '\n'; } template void _show(int i, const T1 &a, const T2 &b, const T3 &...c) { for (; a[i] != ',' && a[i] != '\0'; i++) cerr << a[i]; cerr << ":" << b << " "; _show(i + 1, a, c...); } #line 2 "library/Utility/fastio.hpp" #include namespace fastio { static constexpr uint32_t SZ = 1 << 17; char ibuf[SZ]; char obuf[SZ]; char out[100]; // pointer of ibuf, obuf uint32_t pil = 0, pir = 0, por = 0; struct Pre { char num[10000][4]; constexpr Pre() : num() { for (int i = 0; i < 10000; i++) { int n = i; for (int j = 3; j >= 0; j--) { num[i][j] = n % 10 | '0'; n /= 10; } } } } constexpr pre; inline void load() { memmove(ibuf, ibuf + pil, pir - pil); pir = pir - pil + fread(ibuf + pir - pil, 1, SZ - pir + pil, stdin); pil = 0; if (pir < SZ) ibuf[pir++] = '\n'; } inline void flush() { fwrite(obuf, 1, por, stdout); por = 0; } void rd(char &c) { do { if (pil + 1 > pir) load(); c = ibuf[pil++]; } while (isspace(c)); } void rd(string &x) { x.clear(); char c; do { if (pil + 1 > pir) load(); c = ibuf[pil++]; } while (isspace(c)); do { x += c; if (pil == pir) load(); c = ibuf[pil++]; } while (!isspace(c)); } template void rd_real(T &x) { string s; rd(s); x = stod(s); } template void rd_integer(T &x) { if (pil + 100 > pir) load(); char c; do c = ibuf[pil++]; while (c < '-'); bool minus = 0; if constexpr (is_signed::value || is_same_v) { if (c == '-') { minus = 1, c = ibuf[pil++]; } } x = 0; while ('0' <= c) { x = x * 10 + (c & 15), c = ibuf[pil++]; } if constexpr (is_signed::value || is_same_v) { if (minus) x = -x; } } void rd(int &x) { rd_integer(x); } void rd(ll &x) { rd_integer(x); } void rd(i128 &x) { rd_integer(x); } void rd(uint &x) { rd_integer(x); } void rd(ull &x) { rd_integer(x); } void rd(u128 &x) { rd_integer(x); } void rd(double &x) { rd_real(x); } void rd(long double &x) { rd_real(x); } template void rd(pair &p) { return rd(p.first), rd(p.second); } template void rd_tuple(T &t) { if constexpr (N < std::tuple_size::value) { auto &x = std::get(t); rd(x); rd_tuple(t); } } template void rd(tuple &tpl) { rd_tuple(tpl); } template void rd(array &x) { for (auto &d : x) rd(d); } template void rd(vector &x) { for (auto &d : x) rd(d); } void read() {} template void read(H &h, T &...t) { rd(h), read(t...); } void wt(const char c) { if (por == SZ) flush(); obuf[por++] = c; } void wt(const string s) { for (char c : s) wt(c); } void wt(const char *s) { size_t len = strlen(s); for (size_t i = 0; i < len; i++) wt(s[i]); } template void wt_integer(T x) { if (por > SZ - 100) flush(); if (x < 0) { obuf[por++] = '-', x = -x; } int outi; for (outi = 96; x >= 10000; outi -= 4) { memcpy(out + outi, pre.num[x % 10000], 4); x /= 10000; } if (x >= 1000) { memcpy(obuf + por, pre.num[x], 4); por += 4; } else if (x >= 100) { memcpy(obuf + por, pre.num[x] + 1, 3); por += 3; } else if (x >= 10) { int q = (x * 103) >> 10; obuf[por] = q | '0'; obuf[por + 1] = (x - q * 10) | '0'; por += 2; } else obuf[por++] = x | '0'; memcpy(obuf + por, out + outi + 4, 96 - outi); por += 96 - outi; } template void wt_real(T x) { ostringstream oss; oss << fixed << setprecision(15) << double(x); string s = oss.str(); wt(s); } void wt(int x) { wt_integer(x); } void wt(ll x) { wt_integer(x); } void wt(i128 x) { wt_integer(x); } void wt(uint x) { wt_integer(x); } void wt(ull x) { wt_integer(x); } void wt(u128 x) { wt_integer(x); } void wt(double x) { wt_real(x); } void wt(long double x) { wt_real(x); } template void wt(const pair val) { wt(val.first); wt(' '); wt(val.second); } template void wt_tuple(const T t) { if constexpr (N < std::tuple_size::value) { if constexpr (N > 0) { wt(' '); } const auto x = std::get(t); wt(x); wt_tuple(t); } } template void wt(tuple tpl) { wt_tuple(tpl); } template void wt(const array val) { auto n = val.size(); for (size_t i = 0; i < n; i++) { if (i) wt(' '); wt(val[i]); } } template void wt(const vector val) { auto n = val.size(); for (size_t i = 0; i < n; i++) { if (i) wt(' '); wt(val[i]); } } void print() { wt('\n'); } template void print(Head &&head, Tail &&...tail) { wt(head); if (sizeof...(Tail)) wt(' '); print(forward(tail)...); } void __attribute__((destructor)) _d() { flush(); } } // namespace fastio using fastio::flush; using fastio::print; using fastio::read; inline void first(bool i = true) { print(i ? "first" : "second"); } inline void Alice(bool i = true) { print(i ? "Alice" : "Bob"); } inline void Takahashi(bool i = true) { print(i ? "Takahashi" : "Aoki"); } inline void yes(bool i = true) { print(i ? "yes" : "no"); } inline void Yes(bool i = true) { print(i ? "Yes" : "No"); } inline void No() { print("No"); } inline void YES(bool i = true) { print(i ? "YES" : "NO"); } inline void NO() { print("NO"); } inline void Yay(bool i = true) { print(i ? "Yay!" : ":("); } inline void Possible(bool i = true) { print(i ? "Possible" : "Impossible"); } inline void POSSIBLE(bool i = true) { print(i ? "POSSIBLE" : "IMPOSSIBLE"); } /** * @brief Fast IO */ #line 3 "sol.cpp" #line 2 "library/Math/fastdiv.hpp" struct FastDiv { using u64 = unsigned ll; using u128 = __uint128_t; u128 mod, mh, ml; explicit FastDiv(u64 mod = 1) : mod(mod) { u128 m = u128(-1) / mod; if (m * mod + mod == u128(0)) ++m; mh = m >> 64; ml = m & u64(-1); } u64 umod() const { return mod; } u64 modulo(u128 x) { u128 z = (x & u64(-1)) * ml; z = (x & u64(-1)) * mh + (x >> 64) * ml + (z >> 64); z = (x >> 64) * mh + (z >> 64); x -= z * mod; return x < mod ? x : x - mod; } u64 mul(u64 a, u64 b) { return modulo(u128(a) * b); } }; /** * @brief Fast Division */ #line 2 "library/Math/miller.hpp" struct m64 { using i64 = int64_t; using u64 = uint64_t; using u128 = __uint128_t; static u64 mod; static u64 r; static u64 n2; static u64 get_r() { u64 ret = mod; rep(_,0,5) ret *= 2 - mod * ret; return ret; } static void set_mod(u64 m) { assert(m < (1LL << 62)); assert((m & 1) == 1); mod = m; n2 = -u128(m) % m; r = get_r(); assert(r * mod == 1); } static u64 get_mod() { return mod; } u64 a; m64() : a(0) {} m64(const int64_t &b) : a(reduce((u128(b) + mod) * n2)){}; static u64 reduce(const u128 &b) { return (b + u128(u64(b) * u64(-r)) * mod) >> 64; } u64 get() const { u64 ret = reduce(a); return ret >= mod ? ret - mod : ret; } m64 &operator*=(const m64 &b) { a = reduce(u128(a) * b.a); return *this; } m64 operator*(const m64 &b) const { return m64(*this) *= b; } bool operator==(const m64 &b) const { return (a >= mod ? a - mod : a) == (b.a >= mod ? b.a - mod : b.a); } bool operator!=(const m64 &b) const { return (a >= mod ? a - mod : a) != (b.a >= mod ? b.a - mod : b.a); } m64 pow(u128 n) const { m64 ret(1), mul(*this); while (n > 0) { if (n & 1) ret *= mul; mul *= mul; n >>= 1; } return ret; } }; typename m64::u64 m64::mod, m64::r, m64::n2; bool Miller(ll n){ if(n<2 or (n&1)==0)return (n==2); m64::set_mod(n); ll d=n-1; while((d&1)==0)d>>=1; vector seeds; if(n<(1<<30))seeds={2, 7, 61}; else seeds={2, 325, 9375, 28178, 450775, 9780504}; for(auto& x:seeds){ if(n<=x)break; ll t=d; m64 y=m64(x).pow(t); while(t!=n-1 and y!=1 and y!=n-1){ y*=y; t<<=1; } if(y!=n-1 and (t&1)==0)return 0; } return 1; } /** * @brief Miller-Rabin */ #line 2 "library/Utility/random.hpp" namespace Random { mt19937_64 randgen(chrono::steady_clock::now().time_since_epoch().count()); using u64 = unsigned long long; u64 get() { return randgen(); } template T get(T L) { // [0,L] return get() % (L + 1); } template T get(T L, T R) { // [L,R] return get(R - L) + L; } double uniform() { return double(get(1000000000)) / 1000000000; } string str(int n) { string ret; rep(i, 0, n) ret += get('a', 'z'); return ret; } template void shuffle(Iter first, Iter last) { if (first == last) return; int len = 1; for (auto it = first + 1; it != last; it++) { len++; int j = get(0, len - 1); if (j != len - 1) iter_swap(it, first + j); } } template vector select(int n, T L, T R) { // [L,R] if (n * 2 >= R - L + 1) { vector ret(R - L + 1); iota(ALL(ret), L); shuffle(ALL(ret)); ret.resize(n); return ret; } else { unordered_set used; vector ret; while (SZ(used) < n) { T x = get(L, R); if (!used.count(x)) { used.insert(x); ret.push_back(x); } } return ret; } } void relabel(int n, vector> &es) { shuffle(ALL(es)); vector ord(n); iota(ALL(ord), 0); shuffle(ALL(ord)); for (auto &[u, v] : es) u = ord[u], v = ord[v]; } template vector> genGraph(int n) { vector> cand, es; rep(u, 0, n) rep(v, 0, n) { if (simple and u == v) continue; if (!directed and u > v) continue; cand.push_back({u, v}); } int m = get(SZ(cand)); vector ord; if (simple) ord = select(m, 0, SZ(cand) - 1); else { rep(_, 0, m) ord.push_back(get(SZ(cand) - 1)); } for (auto &i : ord) es.push_back(cand[i]); relabel(n, es); return es; } vector> genTree(int n) { vector> es; rep(i, 1, n) es.push_back({get(i - 1), i}); relabel(n, es); return es; } }; // namespace Random /** * @brief Random */ #line 4 "library/Math/pollard.hpp" vector Pollard(ll n) { if (n <= 1) return {}; if (Miller(n)) return {n}; if ((n & 1) == 0) { vector v = Pollard(n >> 1); v.push_back(2); return v; } for (ll x = 2, y = 2, d;;) { ll c = Random::get(2LL, n - 1); do { x = (__int128_t(x) * x + c) % n; y = (__int128_t(y) * y + c) % n; y = (__int128_t(y) * y + c) % n; d = __gcd(x - y + n, n); } while (d == 1); if (d < n) { vector lb = Pollard(d), rb = Pollard(n / d); lb.insert(lb.end(), ALL(rb)); return lb; } } } /** * @brief Pollard-Rho */ #line 4 "library/Math/primitive.hpp" ll mpow(ll a, i128 t, ll m) { ll res = 1; FastDiv im(m); while (t) { if (t & 1) res = im.mul(res, a); a = im.mul(a, a); t >>= 1; } return res; } ll minv(ll a, ll m) { ll b = m, u = 1, v = 0; while (b) { ll t = a / b; a -= t * b; swap(a, b); u -= t * v; swap(u, v); } u = (u % m + m) % m; return u; } ll getPrimitiveRoot(ll p) { vector ps = Pollard(p - 1); sort(ALL(ps)); rep(x, 1, inf) { for (auto &q : ps) { if (mpow(x, (p - 1) / q, p) == 1) goto fail; } return x; fail:; } assert(0); } ll extgcd(ll a, ll b, ll &p, ll &q) { if (b == 0) { p = 1; q = 0; return a; } ll d = extgcd(b, a % b, q, p); q -= a / b * p; return d; } pair crt(const vector &vs, const vector &ms) { ll V = vs[0], M = ms[0]; rep(i, 1, vs.size()) { ll p, q, v = vs[i], m = ms[i]; if (M < m) swap(M, m), swap(V, v); ll d = extgcd(M, m, p, q); if ((v - V) % d != 0) return {0, -1}; ll md = m / d, tmp = (v - V) / d % md * p % md; V += M * tmp; M *= md; } V = (V % M + M) % M; return {V, M}; } ll ModLog(ll a, ll b, ll p) { ll g = 1; for (ll t = p; t; t >>= 1) g = g * a % p; g = __gcd(g, p); ll t = 1, c = 0; for (; t % g; c++) { if (t == b) return c; t = t * a % p; } if (b % g) return -1; t /= g, b /= g; ll n = p / g, h = 0, gs = 1; for (; h * h < n; h++) gs = gs * a % n; unordered_map bs; for (ll s = 0, e = b; s < h; bs[e] = ++s) e = e * a % n; for (ll s = 0, e = t; s < n;) { e = e * gs % n, s += h; if (bs.count(e)) { return c + s - bs[e]; } } return -1; } ll mod_root(ll k, ll a, ll m) { if (a == 0) return k ? 0 : -1; if (m == 2) return a & 1; k %= m - 1; ll g = gcd(k, m - 1); if (mpow(a, (m - 1) / g, m) != 1) return -1; a = mpow(a, minv(k / g, (m - 1) / g), m); FastDiv im(m); auto _subroot = [&](ll p, int e, ll a) -> ll { // x^(p^e)==a(mod m) ll q = m - 1; int s = 0; while (q % p == 0) { q /= p; s++; } int d = s - e; ll pe = mpow(p, e, m), res = mpow(a, ((pe - 1) * minv(q, pe) % pe * q + 1) / pe, m), c = 1; while (mpow(c, (m - 1) / p, m) == 1) c++; c = mpow(c, q, m); map mp; ll v = 1, block = sqrt(d * p) + 1, bs = mpow(c, mpow(p, s - 1, m - 1) * block % (m - 1), m); rep(i, 0, block + 1) mp[v] = i, v = im.mul(v, bs); ll gs = minv(mpow(c, mpow(p, s - 1, m - 1), m), m); rep(i, 0, d) { ll err = im.mul(a, minv(mpow(res, pe, m), m)); ll pos = mpow(err, mpow(p, d - 1 - i, m - 1), m); rep(j, 0, block + 1) { if (mp.count(pos)) { res = im.mul(res, mpow(c, (block * mp[pos] + j) * mpow(p, i, m - 1) % (m - 1), m)); break; } pos = im.mul(pos, gs); } } return res; }; for (ll d = 2; d * d <= g; d++) if (g % d == 0) { int sz = 0; while (g % d == 0) { g /= d; sz++; } a = _subroot(d, sz, a); } if (g > 1) a = _subroot(g, 1, a); return a; } ull floor_root(ull a, ull k) { if (a <= 1 or k == 1) return a; if (k >= 64) return 1; if (k == 2) return sqrtl(a); constexpr ull LIM = -1; if (a == LIM) a--; auto mul = [&](ull &x, const ull &y) { if (x <= LIM / y) x *= y; else x = LIM; }; auto pw = [&](ull x, ull t) -> ull { ull y = 1; while (t) { if (t & 1) mul(y, x); mul(x, x); t >>= 1; } return y; }; ull ret = (k == 3 ? cbrt(a) - 1 : pow(a, nextafter(1 / double(k), 0))); while (pw(ret + 1, k) <= a) ret++; return ret; } /** * @brief Primitive Function */ #line 6 "sol.cpp" #line 2 "library/Convolution/ntt.hpp" template struct NTT { static constexpr int rank2 = __builtin_ctzll(T::get_mod() - 1); std::array root; // root[i]^(2^i) == 1 std::array iroot; // root[i] * iroot[i] == 1 std::array rate2; std::array irate2; std::array rate3; std::array irate3; NTT() { T g = 2; while (g.pow((T::get_mod() - 1) >> 1) == 1) { g += 1; } root[rank2] = g.pow((T::get_mod() - 1) >> rank2); iroot[rank2] = root[rank2].inv(); for (int i = rank2 - 1; i >= 0; i--) { root[i] = root[i + 1] * root[i + 1]; iroot[i] = iroot[i + 1] * iroot[i + 1]; } { T prod = 1, iprod = 1; for (int i = 0; i <= rank2 - 2; i++) { rate2[i] = root[i + 2] * prod; irate2[i] = iroot[i + 2] * iprod; prod *= iroot[i + 2]; iprod *= root[i + 2]; } } { T prod = 1, iprod = 1; for (int i = 0; i <= rank2 - 3; i++) { rate3[i] = root[i + 3] * prod; irate3[i] = iroot[i + 3] * iprod; prod *= iroot[i + 3]; iprod *= root[i + 3]; } } } void ntt(std::vector &a, bool type = 0) { int n = int(a.size()); int h = __builtin_ctzll((unsigned int)n); a.resize(1 << h); if (type) { int len = h; // a[i, i+(n>>len), i+2*(n>>len), ..] is transformed while (len) { if (len == 1) { int p = 1 << (h - len); T irot = 1; for (int s = 0; s < (1 << (len - 1)); s++) { int offset = s << (h - len + 1); for (int i = 0; i < p; i++) { auto l = a[i + offset]; auto r = a[i + offset + p]; a[i + offset] = l + r; a[i + offset + p] = (unsigned long long)(T::get_mod() + l.v - r.v) * irot.v; ; } if (s + 1 != (1 << (len - 1))) irot *= irate2[__builtin_ctzll(~(unsigned int)(s))]; } len--; } else { // 4-base int p = 1 << (h - len); T irot = 1, iimag = iroot[2]; for (int s = 0; s < (1 << (len - 2)); s++) { T irot2 = irot * irot; T irot3 = irot2 * irot; int offset = s << (h - len + 2); for (int i = 0; i < p; i++) { auto a0 = 1ULL * a[i + offset + 0 * p].v; auto a1 = 1ULL * a[i + offset + 1 * p].v; auto a2 = 1ULL * a[i + offset + 2 * p].v; auto a3 = 1ULL * a[i + offset + 3 * p].v; auto a2na3iimag = 1ULL * T((T::get_mod() + a2 - a3) * iimag.v).v; a[i + offset] = a0 + a1 + a2 + a3; a[i + offset + 1 * p] = (a0 + (T::get_mod() - a1) + a2na3iimag) * irot.v; a[i + offset + 2 * p] = (a0 + a1 + (T::get_mod() - a2) + (T::get_mod() - a3)) * irot2.v; a[i + offset + 3 * p] = (a0 + (T::get_mod() - a1) + (T::get_mod() - a2na3iimag)) * irot3.v; } if (s + 1 != (1 << (len - 2))) irot *= irate3[__builtin_ctzll(~(unsigned int)(s))]; } len -= 2; } } T e = T(n).inv(); for (auto &x : a) x *= e; } else { int len = 0; // a[i, i+(n>>len), i+2*(n>>len), ..] is transformed while (len < h) { if (h - len == 1) { int p = 1 << (h - len - 1); T rot = 1; for (int s = 0; s < (1 << len); s++) { int offset = s << (h - len); for (int i = 0; i < p; i++) { auto l = a[i + offset]; auto r = a[i + offset + p] * rot; a[i + offset] = l + r; a[i + offset + p] = l - r; } if (s + 1 != (1 << len)) rot *= rate2[__builtin_ctzll(~(unsigned int)(s))]; } len++; } else { // 4-base int p = 1 << (h - len - 2); T rot = 1, imag = root[2]; for (int s = 0; s < (1 << len); s++) { T rot2 = rot * rot; T rot3 = rot2 * rot; int offset = s << (h - len); for (int i = 0; i < p; i++) { auto mod2 = 1ULL * T::get_mod() * T::get_mod(); auto a0 = 1ULL * a[i + offset].v; auto a1 = 1ULL * a[i + offset + p].v * rot.v; auto a2 = 1ULL * a[i + offset + 2 * p].v * rot2.v; auto a3 = 1ULL * a[i + offset + 3 * p].v * rot3.v; auto a1na3imag = 1ULL * T(a1 + mod2 - a3).v * imag.v; auto na2 = mod2 - a2; a[i + offset] = a0 + a2 + a1 + a3; a[i + offset + 1 * p] = a0 + a2 + (2 * mod2 - (a1 + a3)); a[i + offset + 2 * p] = a0 + na2 + a1na3imag; a[i + offset + 3 * p] = a0 + na2 + (mod2 - a1na3imag); } if (s + 1 != (1 << len)) rot *= rate3[__builtin_ctzll(~(unsigned int)(s))]; } len += 2; } } } } vector mult(const vector &a, const vector &b) { if (a.empty() or b.empty()) return vector(); int as = a.size(), bs = b.size(); int n = as + bs - 1; if (as <= 30 or bs <= 30) { if (as > 30) return mult(b, a); vector res(n); rep(i, 0, as) rep(j, 0, bs) res[i + j] += a[i] * b[j]; return res; } int m = 1; while (m < n) m <<= 1; vector res(m); rep(i, 0, as) res[i] = a[i]; ntt(res); if (a == b) rep(i, 0, m) res[i] *= res[i]; else { vector c(m); rep(i, 0, bs) c[i] = b[i]; ntt(c); rep(i, 0, m) res[i] *= c[i]; } ntt(res, 1); res.resize(n); return res; } }; /** * @brief Number Theoretic Transform */ #line 2 "library/Math/modint.hpp" template struct fp { unsigned v; static constexpr int get_mod() { return mod; } constexpr unsigned inv() const { assert(v != 0); int x = v, y = mod, p = 1, q = 0, t = 0, tmp = 0; while (y > 0) { t = x / y; x -= t * y, p -= t * q; tmp = x, x = y, y = tmp; tmp = p, p = q, q = tmp; } if (p < 0) p += mod; return p; } constexpr fp(ll x = 0) : v(x >= 0 ? x % mod : (mod - (-x) % mod) % mod) {} fp operator-() const { return fp() - *this; } fp pow(ull t) { fp res = 1, b = *this; while (t) { if (t & 1) res *= b; b *= b; t >>= 1; } return res; } fp &operator+=(const fp &x) { if ((v += x.v) >= mod) v -= mod; return *this; } fp &operator-=(const fp &x) { if ((v += mod - x.v) >= mod) v -= mod; return *this; } fp &operator*=(const fp &x) { v = ull(v) * x.v % mod; return *this; } fp &operator/=(const fp &x) { v = ull(v) * x.inv() % mod; return *this; } fp operator+(const fp &x) const { return fp(*this) += x; } fp operator-(const fp &x) const { return fp(*this) -= x; } fp operator*(const fp &x) const { return fp(*this) *= x; } fp operator/(const fp &x) const { return fp(*this) /= x; } bool operator==(const fp &x) const { return v == x.v; } bool operator!=(const fp &x) const { return v != x.v; } friend istream &operator>>(istream &is, fp &x) { return is >> x.v; } friend ostream &operator<<(ostream &os, const fp &x) { return os << x.v; } }; template void rd(fp &x) { fastio::rd(x.v); } template void wt(fp x) { fastio::wt(x.v); } /** * @brief Modint */ #line 4 "library/Convolution/arbitrary.hpp" using M1 = fp<1045430273>; using M2 = fp<1051721729>; using M3 = fp<1053818881>; NTT N1; NTT N2; NTT N3; constexpr int r_12 = M2(M1::get_mod()).inv(); constexpr int r_13 = M3(M1::get_mod()).inv(); constexpr int r_23 = M3(M2::get_mod()).inv(); constexpr int r_1323 = M3(ll(r_13) * r_23).v; constexpr ll w1 = M1::get_mod(); constexpr ll w2 = ll(w1) * M2::get_mod(); template vector ArbitraryMult(const vector &a, const vector &b) { if (a.empty() or b.empty()) return vector(); int n = a.size() + b.size() - 1; vector res(n); if (min(a.size(), b.size()) <= 60) { rep(i, 0, a.size()) rep(j, 0, b.size()) res[i + j] += T(a[i]) * b[j]; return res; } vector vals[3]; vector a1(ALL(a)), b1(ALL(b)), c1 = N1.mult(a1, b1); vector a2(ALL(a)), b2(ALL(b)), c2 = N2.mult(a2, b2); vector a3(ALL(a)), b3(ALL(b)), c3 = N3.mult(a3, b3); for (M1 x : c1) vals[0].push_back(x.v); for (M2 x : c2) vals[1].push_back(x.v); for (M3 x : c3) vals[2].push_back(x.v); rep(i, 0, n) { ll p = vals[0][i]; ll q = (vals[1][i] + M2::get_mod() - p) * r_12 % M2::get_mod(); ll r = ((vals[2][i] + M3::get_mod() - p) * r_1323 + (M3::get_mod() - q) * r_23) % M3::get_mod(); res[i] = (T(r) * w2 + q * w1 + p); } return res; } template vector ArbitraryMult(const vector &a, const vector &b) { vector A, B; for (auto &x : a) A.push_back(x.v); for (auto &x : b) B.push_back(x.v); return ArbitraryMult(A, B); } /** * @brief Arbitrary Mod Convolution */ #line 3 "library/Math/bigint.hpp" template struct bigint { using u128 = __uint128_t; static const int B = pow(10, D); int sign = 0; vector v; static int get_D() { return D; } static int get_B() { return B; } bigint() {} bigint(const vector &_v, bool _s = false) : sign(_s), v(_v) {} bigint(ll x) { if (x < 0) x *= -1, sign = 1; while (x) { v.push_back(x % B); x /= B; } } bigint(string s) { if (s[0] == '-') s.erase(s.begin()), sign = 1; int add = 0, cnt = 0, base = 1; while (s.size()) { if (cnt == D) { v.push_back(add); cnt = 0; add = 0; base = 1; } add = (s.back() - '0') * base + add; cnt++; base *= 10; s.pop_back(); } if (add) v.push_back(add); } bigint operator-() const { bigint res = *this; res.sign ^= 1; return res; } bigint abs() const { bigint res = *this; res.sign = 0; return res; } int &operator[](const int i) { return v[i]; } int size() const { return v.size(); } void norm() { rep(i, 0, v.size() - 1) { if (v[i] >= 0) { v[i + 1] += v[i] / B; v[i] %= B; } else { int c = (-v[i] + B - 1) / B; v[i] += c * B; v[i + 1] -= c; } } while (!v.empty() and v.back() >= B) { int c = v.back() / B; v.back() %= B; v.push_back(c); } while (!v.empty() and v.back() == 0) v.pop_back(); } string to_str() const { string res; if (v.empty()) return "0"; if (sign) res += '-'; res += to_string(v.back()); for (int i = v.size() - 2; i >= 0; i--) { string add; int w = v[i]; rep(_, 0, D) { add += ('0' + (w % 10)); w /= 10; } reverse(ALL(add)); res += add; } return res; } friend istream &operator>>(istream &is, bigint &x) { string tmp; is >> tmp; x = bigint(tmp); return is; } friend ostream &operator<<(ostream &os, bigint x) { os << x.to_str(); return os; } bigint &operator<<=(const int &x) { if (!v.empty()) { vector add(x, 0); v.insert(v.begin(), ALL(add)); } return *this; } bigint &operator>>=(const int &x) { v = vector(v.begin() + min(x, (int)v.size()), v.end()); return *this; } bigint &operator+=(const bigint &x) { if (sign != x.sign) { *this -= (-x); return *this; } if ((int)v.size() < (int)x.size()) v.resize(x.size(), 0); rep(i, 0, x.size()) { v[i] += x.v[i]; } norm(); return *this; } bigint &operator-=(const bigint &x) { if (sign != x.sign) { *this += (-x); return *this; } if (abs() < x.abs()) { *this = x - (*this); sign ^= 1; return *this; } rep(i, 0, x.size()) { v[i] -= x.v[i]; } norm(); return *this; } bigint &operator*=(const bigint &x) { sign ^= x.sign; auto v1 = ArbitraryMult(v, x.v); u128 add = 0; v.clear(); v.reserve(v1.size() + 3); for (int i = 0;; i++) { if (i >= int(v1.size()) and add == 0) break; if (i < int(v1.size())) add += v1[i]; v.push_back(add % B); add /= B; } norm(); return *this; } bigint div_naive(const bigint &a, const bigint &b) { if (SZ(b) == 1) return a.div(b.v.back()); if (a < b) return bigint(); int norm = B / (b.v.back() + 1); bigint x = a.mul(norm), y = b.mul(norm); int yb = y.v.back(); bigint quo, rem; quo.v.resize(x.size() - y.size() + 1); rem.v = {x.v.end() - y.size(), x.v.end()}; for (int i = x.size() - y.size(); i >= 0; i--) { if (rem.size() == y.size()) { if (rem >= y) { quo[i] = 1; rem -= y; } } else if (rem.size() > y.size()) { ll rb = ll(rem.v.back()) * B + rem[rem.size() - 2]; int q = rb / yb; bigint yq = y.mul(q); while (rem < yq) { q--; yq -= y; } rem -= yq; while (rem >= y) { q++; rem -= y; } quo[i] = q; } if (i) rem.v.insert(rem.v.begin(), x[i - 1]); } return quo; } bigint &operator/=(const bigint &x) { bigint a = abs(), b = x.abs(); sign ^= x.sign; if (a < b) return *this = bigint(); if (b.size() == 1) return *this = a.div(b.v.back()); int deg = a.size() - b.size() + 2, k = deg; while (k > 64) k = (k + 1) >> 1; bigint base(1); base <<= (b.size() + k); bigint inv(div_naive(base, b)); while (k < deg) { bigint y = inv.square(); y.v.insert(y.v.begin(), 0); int l = min(SZ(b), k * 2 + 1); bigint pref; pref.v = {b.v.end() - l, b.v.end()}; y *= pref; y >>= l; inv = inv + inv; inv <<= k + 1; inv -= y; inv.v.erase(inv.v.begin()); k <<= 1; } inv >>= (k - deg); (*this) = a * inv; (*this) >>= int(a.size() + 2); bigint mul = (*this) * b; while (mul + b <= a) { (*this) += bigint(1); mul += b; } while (mul > a) { (*this) -= bigint(1); mul -= b; } return *this; } bigint &operator%=(const bigint &x) { bigint div = (*this) / x; (*this) -= div * x; return *this; } bigint square() const { bigint ret; auto v1 = ArbitraryMult(v, v); u128 add = 0; ret.v.reserve(v1.size() + 3); for (int i = 0;; i++) { if (i >= int(v1.size()) and add == 0) break; if (i < int(v1.size())) add += v1[i]; ret.v.push_back(add % B); add /= B; } return ret; } bigint mul(ll x) const { bigint res; if (x < 0) res.sign ^= 1, x *= -1; u128 add = 0; res.v.reserve(v.size() + 3); for (int i = 0;; i++) { if (i >= int(v.size()) and add == 0) break; if (i < int(v.size())) add += ll(v[i]) * x; res.v.push_back(add % B); add /= B; } return res; } bigint div(ll x) const { bigint res = *this; if (x < 0) res.sign ^= 1, x *= -1; ll add = 0; for (int i = res.v.size() - 1; i >= 0; i--) { add = add * B + res.v[i]; int q = add / x, r = add % x; res.v[i] = q, add = r; } res.norm(); return res; } bigint operator<<(const int &x) const { return bigint(*this) <<= x; } bigint operator>>(const int &x) const { return bigint(*this) >>= x; } bigint operator+(const bigint &x) const { return bigint(*this) += x; } bigint operator-(const bigint &x) const { return bigint(*this) -= x; } bigint operator*(const bigint &x) const { return bigint(*this) *= x; } bigint operator/(const bigint &x) const { return bigint(*this) /= x; } bigint operator%(const bigint &x) const { return bigint(*this) %= x; } bool operator<(const bigint &x) const { if (sign != x.sign) return sign > x.sign; if ((int)v.size() != (int)x.size()) { if (sign) return (int)x.size() < (int)v.size(); else return (int)v.size() < (int)x.size(); } for (int i = v.size() - 1; i >= 0; i--) if (v[i] != x.v[i]) { if (sign) return x.v[i] < v[i]; else return v[i] < x.v[i]; } return false; } bool operator>(const bigint &x) const { return x < *this; } bool operator<=(const bigint &x) const { return !(*this > x); } bool operator>=(const bigint &x) const { return !(*this < x); } bool operator==(const bigint &x) const { return !(*this < x) and !(*this > x); } bool operator!=(const bigint &x) const { return !(*this == x); } }; typedef bigint<9> Bigint; struct Bigfloat { Bigint v; int p = 0; Bigfloat() {} Bigfloat(const ll &_v) { v = Bigint(_v); } Bigfloat(const Bigint &_v, int _p = 0) : v(_v), p(_p) {} void set(int _p) { if (p < _p) { v >>= (_p - p); } else { v <<= (p - _p); } p = _p; } Bigint to_int() const { if (p < 0) return v >> (-p); else return v << p; } Bigfloat &operator+=(const Bigfloat &x) { if (p > x.p) set(x.p), v += x.v; else v += x.v << (x.p - p); return *this; } Bigfloat &operator-=(const Bigfloat &x) { if (p > x.p) set(x.p), v -= x.v; else v -= x.v << (x.p - p); return *this; } Bigfloat square() { Bigfloat res = *this; res.p *= 2; res.v = res.v.square(); return res; } Bigfloat mul(ll x) { Bigfloat res = *this; res.v = v.mul(x); return res; } Bigfloat div(ll x) { Bigfloat res = *this; res.v = v.div(x); return res; } Bigfloat &operator*=(const Bigfloat &x) { p += x.p; v *= x.v; return *this; } Bigfloat &operator/=(const Bigfloat &x) { p -= x.p; v /= x.v; return *this; } Bigfloat operator+(const Bigfloat &x) const { return Bigfloat(*this) += x; } Bigfloat operator-(const Bigfloat &x) const { return Bigfloat(*this) -= x; } Bigfloat operator*(const Bigfloat &x) const { return Bigfloat(*this) *= x; } Bigfloat operator/(const Bigfloat &x) const { return Bigfloat(*this) /= x; } string to_str() { string res = v.abs().to_str(); int d = Bigint::get_D(); if (p * d > 0) res += string(p * d, '0'); else if (-p * d >= 1 and -p * d < (int)res.size()) { res = res.substr(0, (int)res.size() + p * d) + '.' + res.substr((int)res.size() + p * d); } else if (-p * d >= (int)res.size()) res = "0." + string(-p * d - (int)res.size(), '0') + res; if (v.sign) { res.insert(res.begin(), '-'); } return res; } friend ostream &operator<<(ostream &os, Bigfloat x) { os << x.to_str(); return os; } }; Bigfloat sqrt(ll n, int d) { Bigfloat res(Bigint((ll)sqrt(1LL * Bigint::get_B() * Bigint::get_B() / n)), -1), pre; int cur = 1; while (pre.v != res.v) { cur = min(cur << 1, d); pre = res; Bigfloat add = Bigfloat(1) - res.square().mul(n); add.set(-cur); res += (res * add).div(2); res.set(-cur); } return res.mul(n); } /** * @brief Big Integer(Float) */ #line 8 "sol.cpp" Bigint Mpow(Bigint a, Bigint t, Bigint m) { Bigint ret = 1; while (t != 0) { if (t % 2 == 1) ret = (ret * a) % m; a = (a * a) % m; t /= 2; } return ret; } void solve(int _rot) { // print("Case #"+to_string(_rot)+":"); ll B, N, M; read(B, N, M); if (gcd(B, M) != 1) { print(-1); return; } ll M2 = M % B; ll inv = minv(M2, B); ll ord = ModLog(M2, inv, B) + 1; auto ps = Pollard(ord); UNIQUE(ps); bool ch = 0; for (auto &p : ps) if (B % p != 0) { ch = 1; break; } if (ch) { print(-1); return; } if (N == 0) { print(0); return; } { ll n = N * 2 + 10; Bigint Bn = 1, BnN = 1; rep(_, 0, n) Bn *= B; rep(_, 0, n + N) BnN *= B; Bigint mpw = Mpow(M, Bn, BnN); Bigint a = (mpw - 1) / Bn; print(a.to_str()); } } int main() { int t; read(t); rep(rot, 0, t) solve(rot + 1); return 0; }