// #pragma GCC target("avx2") // #pragma GCC optimize("O3") // #pragma GCC optimize("unroll-loops") // #define INTERACTIVE #include using namespace std; namespace templates { // type using ll = long long; using ull = unsigned long long; using Pii = pair; using Pil = pair; using Pli = pair; using Pll = pair; template using pq = priority_queue; template using qp = priority_queue, greater>; // clang-format off #define vec(T, A, ...) vector A(__VA_ARGS__); #define vvec(T, A, h, ...) vector> A(h, vector(__VA_ARGS__)); #define vvvec(T, A, h1, h2, ...) vector>> A(h1, vector>(h2, vector(__VA_ARGS__))); // clang-format on // for loop #define fori1(a) for (ll _ = 0; _ < (a); _++) #define fori2(i, a) for (ll i = 0; i < (a); i++) #define fori3(i, a, b) for (ll i = (a); i < (b); i++) #define fori4(i, a, b, c) for (ll i = (a); ((c) > 0 || i > (b)) && ((c) < 0 || i < (b)); i += (c)) #define overload4(a, b, c, d, e, ...) e #define fori(...) overload4(__VA_ARGS__, fori4, fori3, fori2, fori1)(__VA_ARGS__) // declare and input // clang-format off #define INT(...) int __VA_ARGS__; inp(__VA_ARGS__); #define LL(...) ll __VA_ARGS__; inp(__VA_ARGS__); #define STRING(...) string __VA_ARGS__; inp(__VA_ARGS__); #define CHAR(...) char __VA_ARGS__; inp(__VA_ARGS__); #define DOUBLE(...) double __VA_ARGS__; STRING(str___); __VA_ARGS__ = stod(str___); #define VEC(T, A, n) vector A(n); inp(A); #define VVEC(T, A, n, m) vector> A(n, vector(m)); inp(A); // clang-format on // const value const ll MOD1 = 1000000007; const ll MOD9 = 998244353; const double PI = acos(-1); // other macro #if !defined(RIN__LOCAL) && !defined(INTERACTIVE) #define endl "\n" #endif #define spa ' ' #define len(A) ll(A.size()) #define all(A) begin(A), end(A) // function vector stoc(string &S) { int n = S.size(); vector ret(n); for (int i = 0; i < n; i++) ret[i] = S[i]; return ret; } string ctos(vector &S) { int n = S.size(); string ret = ""; for (int i = 0; i < n; i++) ret += S[i]; return ret; } template auto min(const T &a) { return *min_element(all(a)); } template auto max(const T &a) { return *max_element(all(a)); } template auto clamp(T &a, const S &l, const S &r) { return (a > r ? r : a < l ? l : a); } template inline bool chmax(T &a, const S &b) { return (a < b ? a = b, 1 : 0); } template inline bool chmin(T &a, const S &b) { return (a > b ? a = b, 1 : 0); } template inline bool chclamp(T &a, const S &l, const S &r) { auto b = clamp(a, l, r); return (a != b ? a = b, 1 : 0); } template T sum(vector &A) { T tot = 0; for (auto a : A) tot += a; return tot; } template vector compression(vector X) { sort(all(X)); X.erase(unique(all(X)), X.end()); return X; } // input and output namespace io { // __int128_t std::ostream &operator<<(std::ostream &dest, __int128_t value) { std::ostream::sentry s(dest); if (s) { __uint128_t tmp = value < 0 ? -value : value; char buffer[128]; char *d = std::end(buffer); do { --d; *d = "0123456789"[tmp % 10]; tmp /= 10; } while (tmp != 0); if (value < 0) { --d; *d = '-'; } int len = std::end(buffer) - d; if (dest.rdbuf()->sputn(d, len) != len) { dest.setstate(std::ios_base::badbit); } } return dest; } // vector template istream &operator>>(istream &is, vector &A) { for (auto &a : A) is >> a; return is; } template ostream &operator<<(ostream &os, vector &A) { for (size_t i = 0; i < A.size(); i++) { os << A[i]; if (i != A.size() - 1) os << ' '; } return os; } // vector> template istream &operator>>(istream &is, vector> &A) { for (auto &a : A) is >> a; return is; } template ostream &operator<<(ostream &os, vector> &A) { for (size_t i = 0; i < A.size(); i++) { os << A[i]; if (i != A.size() - 1) os << endl; } return os; } // pair template istream &operator>>(istream &is, pair &A) { is >> A.first >> A.second; return is; } template ostream &operator<<(ostream &os, pair &A) { os << A.first << ' ' << A.second; return os; } // vector> template istream &operator>>(istream &is, vector> &A) { for (size_t i = 0; i < A.size(); i++) { is >> A[i]; } return is; } template ostream &operator<<(ostream &os, vector> &A) { for (size_t i = 0; i < A.size(); i++) { os << A[i]; if (i != A.size() - 1) os << endl; } return os; } // tuple template struct TuplePrint { static ostream &print(ostream &os, const T &t) { TuplePrint::print(os, t); os << ' ' << get(t); return os; } }; template struct TuplePrint { static ostream &print(ostream &os, const T &t) { os << get<0>(t); return os; } }; template ostream &operator<<(ostream &os, const tuple &t) { TuplePrint::print(os, t); return os; } // io functions void FLUSH() { cout << flush; } void print() { cout << endl; } template void print(Head &&head, Tail &&...tail) { cout << head; if (sizeof...(Tail)) cout << spa; print(std::forward(tail)...); } template void prisep(vector &A, S sep) { int n = A.size(); for (int i = 0; i < n; i++) { cout << A[i]; if (i != n - 1) cout << sep; } cout << endl; } template void priend(T A, S end) { cout << A << end; } template void prispa(T A) { priend(A, spa); } template bool printif(bool f, T A, S B) { if (f) print(A); else print(B); return f; } template void inp(T &...a) { (cin >> ... >> a); } } // namespace io using namespace io; // read graph vector> read_edges(int n, int m, bool direct = false, int indexed = 1) { vector> edges(n, vector()); for (int i = 0; i < m; i++) { INT(u, v); u -= indexed; v -= indexed; edges[u].push_back(v); if (!direct) edges[v].push_back(u); } return edges; } vector> read_tree(int n, int indexed = 1) { return read_edges(n, n - 1, false, indexed); } template vector>> read_wedges(int n, int m, bool direct = false, int indexed = 1) { vector>> edges(n, vector>()); for (int i = 0; i < m; i++) { INT(u, v); T w; inp(w); u -= indexed; v -= indexed; edges[u].push_back({v, w}); if (!direct) edges[v].push_back({u, w}); } return edges; } template vector>> read_wtree(int n, int indexed = 1) { return read_wedges(n, n - 1, false, indexed); } // yes / no namespace yesno { // yes inline bool yes(bool f = true) { cout << (f ? "yes" : "no") << endl; return f; } inline bool Yes(bool f = true) { cout << (f ? "Yes" : "No") << endl; return f; } inline bool YES(bool f = true) { cout << (f ? "YES" : "NO") << endl; return f; } // no inline bool no(bool f = true) { cout << (!f ? "yes" : "no") << endl; return f; } inline bool No(bool f = true) { cout << (!f ? "Yes" : "No") << endl; return f; } inline bool NO(bool f = true) { cout << (!f ? "YES" : "NO") << endl; return f; } // possible inline bool possible(bool f = true) { cout << (f ? "possible" : "impossible") << endl; return f; } inline bool Possible(bool f = true) { cout << (f ? "Possible" : "Impossible") << endl; return f; } inline bool POSSIBLE(bool f = true) { cout << (f ? "POSSIBLE" : "IMPOSSIBLE") << endl; return f; } // impossible inline bool impossible(bool f = true) { cout << (!f ? "possible" : "impossible") << endl; return f; } inline bool Impossible(bool f = true) { cout << (!f ? "Possible" : "Impossible") << endl; return f; } inline bool IMPOSSIBLE(bool f = true) { cout << (!f ? "POSSIBLE" : "IMPOSSIBLE") << endl; return f; } // Alice Bob inline bool Alice(bool f = true) { cout << (f ? "Alice" : "Bob") << endl; return f; } inline bool Bob(bool f = true) { cout << (f ? "Bob" : "Alice") << endl; return f; } // Takahashi Aoki inline bool Takahashi(bool f = true) { cout << (f ? "Takahashi" : "Aoki") << endl; return f; } inline bool Aoki(bool f = true) { cout << (f ? "Aoki" : "Takahashi") << endl; return f; } } // namespace yesno using namespace yesno; } // namespace templates using namespace templates; template struct Modint { int x; Modint() : x(0) {} Modint(int64_t y) { if (y >= 0) x = y % MOD; else x = (y % MOD + MOD) % MOD; } Modint &operator+=(const Modint &p) { x += p.x; if (x >= MOD) x -= MOD; return *this; } Modint &operator-=(const Modint &p) { x -= p.x; if (x < 0) x += MOD; return *this; } Modint &operator*=(const Modint &p) { x = int(1LL * x * p.x % MOD); return *this; } Modint &operator/=(const Modint &p) { *this *= p.inverse(); return *this; } Modint &operator%=(const Modint &p) { assert(p.x == 0); return *this; } Modint operator-() const { return Modint(-x); } Modint &operator++() { x++; if (x == MOD) x = 0; return *this; } Modint &operator--() { if (x == 0) x = MOD; x--; return *this; } Modint operator++(int) { Modint result = *this; ++*this; return result; } Modint operator--(int) { Modint result = *this; --*this; return result; } friend Modint operator+(const Modint &lhs, const Modint &rhs) { return Modint(lhs) += rhs; } friend Modint operator-(const Modint &lhs, const Modint &rhs) { return Modint(lhs) -= rhs; } friend Modint operator*(const Modint &lhs, const Modint &rhs) { return Modint(lhs) *= rhs; } friend Modint operator/(const Modint &lhs, const Modint &rhs) { return Modint(lhs) /= rhs; } friend Modint operator%(const Modint &lhs, const Modint &rhs) { assert(rhs.x == 0); return Modint(lhs); } bool operator==(const Modint &p) const { return x == p.x; } bool operator!=(const Modint &p) const { return x != p.x; } bool operator<(const Modint &rhs) const { return x < rhs.x; } bool operator<=(const Modint &rhs) const { return x <= rhs.x; } bool operator>(const Modint &rhs) const { return x > rhs.x; } bool operator>=(const Modint &rhs) const { return x >= rhs.x; } Modint inverse() const { int a = x, b = MOD, u = 1, v = 0, t; while (b > 0) { t = a / b; a -= t * b; u -= t * v; std::swap(a, b); std::swap(u, v); } return Modint(u); } Modint pow(int64_t k) const { Modint ret(1); Modint y(x); while (k > 0) { if (k & 1) ret *= y; y *= y; k >>= 1; } return ret; } std::pair to_frac(int max_n = 1000) const { int y = x; for (int i = 1; i <= max_n; i++) { if (y <= max_n) { return {y, i}; } else if (MOD - y <= max_n) { return {-(MOD - y), i}; } y = (y + x) % MOD; } return {-1, -1}; } friend std::ostream &operator<<(std::ostream &os, const Modint &p) { return os << p.x; } friend std::istream &operator>>(std::istream &is, Modint &p) { int64_t y; is >> y; p = Modint(y); return (is); } static int get_mod() { return MOD; } }; struct Arbitrary_Modint { int x; static int MOD; static void set_mod(int mod) { MOD = mod; } Arbitrary_Modint() : x(0) {} Arbitrary_Modint(int64_t y) { if (y >= 0) x = y % MOD; else x = (y % MOD + MOD) % MOD; } Arbitrary_Modint &operator+=(const Arbitrary_Modint &p) { x += p.x; if (x >= MOD) x -= MOD; return *this; } Arbitrary_Modint &operator-=(const Arbitrary_Modint &p) { x -= p.x; if (x < 0) x += MOD; return *this; } Arbitrary_Modint &operator*=(const Arbitrary_Modint &p) { x = int(1LL * x * p.x % MOD); return *this; } Arbitrary_Modint &operator/=(const Arbitrary_Modint &p) { *this *= p.inverse(); return *this; } Arbitrary_Modint &operator%=(const Arbitrary_Modint &p) { assert(p.x == 0); return *this; } Arbitrary_Modint operator-() const { return Arbitrary_Modint(-x); } Arbitrary_Modint &operator++() { x++; if (x == MOD) x = 0; return *this; } Arbitrary_Modint &operator--() { if (x == 0) x = MOD; x--; return *this; } Arbitrary_Modint operator++(int) { Arbitrary_Modint result = *this; ++*this; return result; } Arbitrary_Modint operator--(int) { Arbitrary_Modint result = *this; --*this; return result; } friend Arbitrary_Modint operator+(const Arbitrary_Modint &lhs, const Arbitrary_Modint &rhs) { return Arbitrary_Modint(lhs) += rhs; } friend Arbitrary_Modint operator-(const Arbitrary_Modint &lhs, const Arbitrary_Modint &rhs) { return Arbitrary_Modint(lhs) -= rhs; } friend Arbitrary_Modint operator*(const Arbitrary_Modint &lhs, const Arbitrary_Modint &rhs) { return Arbitrary_Modint(lhs) *= rhs; } friend Arbitrary_Modint operator/(const Arbitrary_Modint &lhs, const Arbitrary_Modint &rhs) { return Arbitrary_Modint(lhs) /= rhs; } friend Arbitrary_Modint operator%(const Arbitrary_Modint &lhs, const Arbitrary_Modint &rhs) { assert(rhs.x == 0); return Arbitrary_Modint(lhs); } bool operator==(const Arbitrary_Modint &p) const { return x == p.x; } bool operator!=(const Arbitrary_Modint &p) const { return x != p.x; } bool operator<(const Arbitrary_Modint &rhs) { return x < rhs.x; } bool operator<=(const Arbitrary_Modint &rhs) { return x <= rhs.x; } bool operator>(const Arbitrary_Modint &rhs) { return x > rhs.x; } bool operator>=(const Arbitrary_Modint &rhs) { return x >= rhs.x; } Arbitrary_Modint inverse() const { int a = x, b = MOD, u = 1, v = 0, t; while (b > 0) { t = a / b; a -= t * b; u -= t * v; std::swap(a, b); std::swap(u, v); } return Arbitrary_Modint(u); } Arbitrary_Modint pow(int64_t k) const { Arbitrary_Modint ret(1); Arbitrary_Modint y(x); while (k > 0) { if (k & 1) ret *= y; y *= y; k >>= 1; } return ret; } friend std::ostream &operator<<(std::ostream &os, const Arbitrary_Modint &p) { return os << p.x; } friend std::istream &operator>>(std::istream &is, Arbitrary_Modint &p) { int64_t y; is >> y; p = Arbitrary_Modint(y); return (is); } static int get_mod() { return MOD; } }; int Arbitrary_Modint::MOD = 998244353; using modint9 = Modint<998244353>; using modint1 = Modint<1000000007>; using modint = Arbitrary_Modint; using mint = modint9; template struct NumberTheoreticTransform { static std::vector roots, iroots, rate3, irate3; static int max_base; NumberTheoreticTransform() = default; static void init() { if (!roots.empty()) return; const unsigned mod = mint::get_mod(); auto tmp = mod - 1; max_base = 0; while (tmp % 2 == 0) { tmp >>= 1; max_base++; } mint root = 2; while (root.pow((mod - 1) >> 1) == 1) root++; roots.resize(max_base + 1); iroots.resize(max_base + 1); rate3.resize(max_base + 1); irate3.resize(max_base + 1); roots[max_base] = root.pow((mod - 1) >> max_base); iroots[max_base] = mint(1) / roots[max_base]; for (int i = max_base - 1; i >= 0; i--) { roots[i] = roots[i + 1] * roots[i + 1]; iroots[i] = iroots[i + 1] * iroots[i + 1]; } mint prod = 1, iprod = 1; for (int i = 0; i <= max_base - 3; i++) { rate3[i] = roots[i + 3] * prod; irate3[i] = iroots[i + 3] * iprod; prod *= iroots[i + 3]; iprod *= roots[i + 3]; } } static void ntt(std::vector &A) { init(); int n = int(A.size()); int h = __builtin_ctz(n); int le = 0; mint imag = roots[2]; if (h & 1) { int p = 1 << (h - 1); for (int i = 0; i < p; i++) { auto r = A[i + p]; A[i + p] = A[i] - r; A[i] += r; } le++; } for (; le + 1 < h; le += 2) { int p = 1 << (h - le - 2); for (int i = 0; i < p; i++) { auto a0 = A[i]; auto a1 = A[i + p]; auto a2 = A[i + 2 * p]; auto a3 = A[i + 3 * p]; auto a1na3imag = (a1 - a3) * imag; A[i] = a0 + a2 + a1 + a3; A[i + p] = a0 + a2 - (a1 + a3); A[i + 2 * p] = a0 - a2 + a1na3imag; A[i + 3 * p] = a0 - a2 - a1na3imag; } mint rot = rate3[0]; for (int s = 1; s < (1 << le); s++) { int offset = s << (h - le); mint rot2 = rot * rot; mint rot3 = rot2 * rot; for (int i = 0; i < p; i++) { auto a0 = A[i + offset]; auto a1 = A[i + offset + p] * rot; auto a2 = A[i + offset + 2 * p] * rot2; auto a3 = A[i + offset + 3 * p] * rot3; auto a1na3imag = (a1 - a3) * imag; A[i + offset] = a0 + a2 + a1 + a3; A[i + offset + p] = a0 + a2 - (a1 + a3); A[i + offset + 2 * p] = a0 - a2 + a1na3imag; A[i + offset + 3 * p] = a0 - a2 - a1na3imag; } rot *= rate3[__builtin_ctz(~s)]; } } } static void intt(std::vector &A, bool f = true) { init(); int n = int(A.size()); int h = __builtin_ctz(n); int le = h; mint iimag = iroots[2]; for (; le > 1; le -= 2) { int p = 1 << (h - le); for (int i = 0; i < p; i++) { auto a0 = A[i]; auto a1 = A[i + p]; auto a2 = A[i + 2 * p]; auto a3 = A[i + 3 * p]; auto a2na3iimag = (a2 - a3) * iimag; A[i] = a0 + a1 + a2 + a3; A[i + p] = a0 - a1 + a2na3iimag; A[i + 2 * p] = a0 + a1 - (a2 + a3); A[i + 3 * p] = a0 - a1 - a2na3iimag; } mint irot = irate3[0]; for (int s = 1; s < (1 << (le - 2)); s++) { int offset = s << (h - le + 2); mint irot2 = irot * irot; mint irot3 = irot2 * irot; for (int i = 0; i < p; i++) { auto a0 = A[i + offset]; auto a1 = A[i + offset + p]; auto a2 = A[i + offset + 2 * p]; auto a3 = A[i + offset + 3 * p]; auto a2na3iimag = (a2 - a3) * iimag; A[i + offset] = a0 + a1 + a2 + a3; A[i + offset + p] = (a0 - a1 + a2na3iimag) * irot; A[i + offset + 2 * p] = (a0 + a1 - (a2 + a3)) * irot2; A[i + offset + 3 * p] = (a0 - a1 - a2na3iimag) * irot3; } irot *= irate3[__builtin_ctz(~s)]; } } if (le >= 1) { int p = 1 << (h - 1); for (int i = 0; i < p; i++) { auto ajp = A[i] - A[i + p]; A[i] += A[i + p]; A[i + p] = ajp; } } if (f) { mint inv = mint(1) / n; for (int i = 0; i < n; i++) { A[i] *= inv; } } } static std::vector multiply(std::vector A, std::vector B) { int need = int(A.size() + B.size()) - 1; if (std::min(A.size(), B.size()) < 60u) { std::vector C(need, 0); for (size_t i = 0; i < A.size(); i++) for (size_t j = 0; j < B.size(); j++) { C[i + j] += A[i] * B[j]; } return C; } int sz = 1; while (sz < need) sz <<= 1; A.resize(sz, 0); B.resize(sz, 0); ntt(A); ntt(B); mint inv = mint(1) / sz; for (int i = 0; i < sz; i++) A[i] *= B[i] * inv; intt(A, false); A.resize(need); return A; } }; template std::vector NumberTheoreticTransform::roots = std::vector(); template std::vector NumberTheoreticTransform::iroots = std::vector(); template std::vector NumberTheoreticTransform::rate3 = std::vector(); template std::vector NumberTheoreticTransform::irate3 = std::vector(); template int NumberTheoreticTransform::max_base = 0; template T modpow(T a, S b, T MOD) { T ret = 1; while (b > 0) { if (b & 1) { ret *= a; ret %= MOD; } a *= a; a %= MOD; b >>= 1; } return ret; } template T cipolla(T a, T MOD) { if (MOD == 2) return a; else if (a == 0) return 0; else if (modpow(a, (MOD - 1) / 2, MOD) != 1) return -1; T b = 1; while (modpow((b * b + MOD - a) % MOD, (MOD - 1) / 2, MOD) == 1) { b++; } T base = (b * b + MOD - a) % MOD; auto multi = [&](T a0, T b0, T a1, T b1) -> std::pair { return {(a0 * a1 + (b0 * b1 % MOD) * base) % MOD, (a0 * b1 + b0 * a1) % MOD}; }; auto pow_ = [&](auto self, T a, T b, T n) -> std::pair { if (n == 0) return {1, 0}; auto tmp = multi(a, b, a, b); auto ret = self(self, tmp.first, tmp.second, n / 2); if (n & 1) { ret = multi(ret.first, ret.second, a, b); } return ret; }; return pow_(pow_, b, 1LL, (MOD + 1) / 2).first; } template struct FormalPowerSeries : std::vector { using std::vector::vector; using FPS = FormalPowerSeries; static std::vector inv_x; void shrink() { while (this->size() && this->back() == mint(0)) { this->pop_back(); } } FPS &operator+=(const FPS &A) { if (A.size() > this->size()) this->resize(A.size()); for (size_t i = 0; i < A.size(); i++) (*this)[i] += A[i]; return *this; } FPS &operator+=(const mint &x) { if (this->empty()) this->resize(1); (*this)[0] += x; return *this; } FPS &operator-=(const FPS &A) { if (A.size() > this->size()) this->resize(A.size()); for (size_t i = 0; i < A.size(); i++) (*this)[i] -= A[i]; return *this; } FPS &operator-=(const mint &x) { if (this->empty()) this->resize(1); (*this)[0] -= x; return *this; } FPS &operator*=(const FPS &A) { if (this->empty() || A.empty()) { this->clear(); return *this; } auto res = NumberTheoreticTransform::multiply(*this, A); return *this = {begin(res), end(res)}; } FPS &operator*=(const mint &x) { for (size_t i = 0; i < this->size(); i++) (*this)[i] *= x; return *this; } FPS operator+(const FPS &A) const { return FPS(*this) += A; } FPS operator+(const mint &x) const { return FPS(*this) += x; } FPS operator-(const FPS &A) const { return FPS(*this) -= A; } FPS operator-(const mint &x) const { return FPS(*this) -= x; } FPS operator*(const FPS &A) const { return FPS(*this) *= A; } FPS operator*(const mint &x) const { return FPS(*this) *= x; } FPS operator-() const { FPS ret(this->size); for (size_t i = 0; i < this->size(); i++) ret[i] = -(*this)[i]; return ret; } FPS inv(int deg = -1) { assert((*this)[0] != mint(0)); if (deg == -1) deg = this->size(); FPS g = {mint(1) / (*this)[0]}; int l = 1; while (l < deg) { FPS tmp = g * 2; l <<= 1; FPS tmp2; g *= g; if (int(this->size()) >= l) tmp2 = FPS({this->begin(), this->begin() + l}) * g; else tmp2 = (*this) * g; g = tmp - tmp2; g.resize(l); } g.resize(deg); return g; } void iinv(int deg = -1) { *this = inv(deg); } FPS differential() { FPS ret(this->size() - 1); for (size_t i = 0; i < this->size() - 1; i++) ret[i] = (*this)[i + 1] * (i + 1); return ret; } void idifferential() { *this = this->differential(); } void extend_inv(int n) { int bn = inv_x.size(); if (n >= bn) { inv_x.resize(n + 1, 0); if (bn == 0) { inv_x[0] = 0; inv_x[1] = 1; bn = 2; } long long mod = mint::get_mod(); for (int i = bn; i <= n; i++) { inv_x[i] = mod - inv_x[mod % i].x * (mod / i) % mod; } } } FPS integral() { extend_inv(this->size()); FPS ret(this->size() + 1); for (size_t i = 0; i < this->size(); i++) ret[i + 1] = (*this)[i] * inv_x[i + 1]; return ret; } void iintegral() { *this = this->integral(); } FPS log(int deg = -1) { assert((*this)[0] == mint(1)); if (deg == -1) deg = this->size(); FPS B = (this->differential()) * (this->inv()); B.resize(deg - 1); return B.integral(); } void ilog(int deg = -1) { *this = this->log(deg); } FPS exp(int deg = -1) { assert((*this)[0] == mint(0)); if (deg == -1) deg = this->size(); FPS g = {1}; int l = 1; while (l < deg * 2) { l *= 2; FPS tmp = {1}; tmp -= g.log(l); if (int(this->size()) >= l) tmp += FPS({this->begin(), this->begin() + l}); else tmp += (*this); g *= tmp; g.resize(l); } g.resize(deg); return g; } void iexp(int deg = -1) { *this = this->exp(deg); } FPS pow(long long k, int deg = -1) { if (deg == -1) deg = this->size(); if (k == 0) { FPS ret(deg, 0); ret[0] = 1; return ret; } int p = -1; for (int i = 0; i < deg; i++) { if ((*this)[i] != 0) { p = i; break; } } if (p == -1 || p > deg / k) { FPS ret(deg, 0); return ret; } mint inv = mint(1) / (*this)[p]; FPS A = FPS({(*this).begin() + p, (*this).end()}); A *= inv; A.ilog(deg); A *= k % mint::get_mod(); A.iexp(deg); FPS B(p * k, 0); B.insert(B.end(), A.begin(), A.begin() + (deg - p * k)); B *= (*this)[p].pow(k); return B; } void ipow(long long k, int deg = -1) { *this = this->pow(k, deg); } FPS sqrt(int deg = -1) { if (deg == -1) deg = this->size(); if (this->size() == 0u) { FPS ret(deg, 0); return ret; } if ((*this)[0] == mint(0)) { for (size_t i = 1; i < this->size(); i++) { if ((*this)[i] != 0) { if (i & 1) { FPS ret; return ret; } if (deg <= int(i / 2)) break; FPS ret = FPS({this->begin() + i, this->end()}).sqrt(deg - i / 2); if (ret.size() == 0u) return ret; FPS ret2(i / 2, 0); ret2.insert(ret2.end(), ret.begin(), ret.end()); std::swap(ret, ret2); if (int(ret.size()) < deg) ret.resize(deg); return ret; } } FPS ret(deg, 0); return ret; } long long sq = cipolla((*this)[0].x, mint::get_mod()); if (sq == -1) { FPS ret; return ret; } mint inv2 = mint(1) / 2; FPS g = {sq}; int l = 1; while (l < deg) { l *= 2; if (int(this->size()) >= l) g += FPS({this->begin(), this->begin() + l}) * g.inv(l); else g += (*this) * g.inv(l); g *= inv2; } g.resize(deg); return g; } void isqrt(int deg = -1) { *this = this->sqrt(deg); } FPS taylorshift(mint a) { auto A = (*this); int deg = A.size(); extend_inv(deg); mint fac = 1; for (int i = 0; i < deg; i++) { A[i] *= fac; fac *= (i + 1); } reverse(A.begin(), A.end()); FPS g(deg, 0); g[0] = 1; for (int i = 1; i < deg; i++) g[i] = g[i - 1] * a * inv_x[i]; A *= g; if (int(A.size()) > deg) A.resize(deg); reverse(A.begin(), A.end()); mint invfac = 1; for (int i = 0; i < deg; i++) { A[i] *= invfac; invfac *= inv_x[i + 1]; } return A; } void itaylorshift(mint a) { int deg = this->size(); extend_inv(deg); mint fac = 1; for (int i = 0; i < deg; i++) { (*this)[i] *= fac; fac *= (i + 1); } reverse(this->begin(), this->end()); FPS g(deg, 0); g[0] = 1; for (int i = 1; i < deg; i++) g[i] = g[i - 1] * a * inv_x[i]; (*this) *= g; if (int(this->size()) > deg) this->resize(deg); reverse(this->begin(), this->end()); mint invfac = 1; for (int i = 0; i < deg; i++) { (*this)[i] *= invfac; invfac *= inv_x[i + 1]; } } std::pair division_of_polynomial(FPS G) { FPS F = *this; if (F.size() < G.size()) { return {{}, F}; } reverse(F.begin(), F.end()); reverse(G.begin(), G.end()); int deg = F.size() - G.size() + 1u; auto Q = F * G.inv(deg); if (int(Q.size()) > deg) Q.resize(deg); reverse(Q.begin(), Q.end()); reverse(F.begin(), F.end()); reverse(G.begin(), G.end()); auto R = F - G * Q; R.shrink(); return {Q, R}; } std::vector multipoint_evaluation(std::vector &X) { int m = X.size(); int m2 = 1; while (m2 <= m - 1) m2 *= 2; std::vector G(m2 << 1, FPS(1, 1)); for (int i = 0; i < m; i++) G[m2 + i] = {-X[i], 1}; for (int i = m2 - 1; i >= 0; i--) G[i] = G[i << 1] * G[(i << 1) | 1]; G[1] = this->division_of_polynomial(G[1]).second; for (int i = 2; i < m2 + m; i++) G[i] = G[i >> 1].division_of_polynomial(G[i]).second; std::vector Y(m); for (int i = 0; i < m; i++) { if (G[m2 + i].empty()) Y[i] = 0; else Y[i] = G[m2 + i][0]; } return Y; } std::vector multipoint_evaluation(std::vector &X) { int m = X.size(); int m2 = 1; while (m2 <= m - 1) m2 *= 2; std::vector G(m2 << 1, FPS(1, 1)); for (int i = 0; i < m; i++) G[m2 + i] = {-X[i], 1}; for (int i = m2 - 1; i >= 0; i--) G[i] = G[i << 1] * G[(i << 1) | 1]; G[1] = this->division_of_polynomial(G[1]).second; for (int i = 2; i < m2 + m; i++) G[i] = G[i >> 1].division_of_polynomial(G[i]).second; std::vector Y(m); for (int i = 0; i < m; i++) { if (G[m2 + i].empty()) Y[i] = 0; else Y[i] = G[m2 + i][0].x; } return Y; } friend std::ostream &operator<<(std::ostream &os, const FPS &A) { for (size_t i = 0; i < A.size(); i++) { os << A[i]; if (i != A.size() - 1) os << ' '; } return os; } friend std::istream &operator>>(std::istream &is, FPS &A) { for (size_t i = 0; i < A.size(); i++) { is >> A[i]; } return (is); } }; template std::vector FormalPowerSeries::inv_x = std::vector(); using FPS = FormalPowerSeries; bool isPrime(long long n) { if (n <= 1) return false; else if (n == 2) return true; else if (n % 2 == 0) return false; long long A[7] = {2, 325, 9375, 28178, 450775, 9780504, 1795265022}; long long s = 0; long long d = n - 1; while (d % 2 == 0) { d /= 2; s++; } for (auto a : A) { if (a % n == 0) return true; long long x = modpow<__int128_t>(a, d, n); if (x != 1) { bool ng = true; for (int i = 0; i < s; i++) { if (x == n - 1) { ng = false; break; }; x = __int128_t(x) * x % n; } if (ng) return false; } } return true; } long long pollard(long long N) { if (N % 2 == 0) return 2; if (isPrime(N)) return N; long long step = 0; auto f = [&](long long x) -> long long { return (__int128_t(x) * x + step) % N; }; while (true) { ++step; long long x = step, y = f(x); while (true) { long long p = std::gcd(y - x + N, N); if (p == 0 || p == N) break; if (p != 1) return p; x = f(x); y = f(f(y)); } } } std::vector primefact(long long N) { if (N == 1) return {}; long long p = pollard(N); if (p == N) return {p}; std::vector left = primefact(p); std::vector right = primefact(N / p); left.insert(left.end(), right.begin(), right.end()); std::sort(left.begin(), left.end()); return left; } template mint BostanMori(std::vector P, std::vector Q, unsigned long long n) { while (n > 0) { std::vector R(Q.size()); for (size_t i = 0; i < Q.size(); i++) { if (i & 1) R[i] = -Q[i]; else R[i] = Q[i]; } Q = NumberTheoreticTransform::multiply(Q, R); int lq = Q.size(); for (int i = 0; i < lq; i += 2) { Q[i / 2] = Q[i]; } Q.resize((lq + 1) / 2); P = NumberTheoreticTransform::multiply(P, R); if (n & 1) { int lp = P.size(); for (int i = 1; i < lp; i += 2) { P[i / 2] = P[i]; } P.resize(lp / 2); } else { int lp = P.size(); for (int i = 0; i < lp; i += 2) { P[i / 2] = P[i]; } P.resize((lp + 1) / 2); } n >>= 1; } return P[0] / Q[0]; } vector primes = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97}; void solve() { LL(n, s); if (n == 1) { print(mint(s)); return; } vec(ll, P, 0); ll ss = n; for (auto p : primes) { while (ss % p == 0) { P.push_back(p); ss /= p; } } vec(Pll, rle, 0); int m = 0; { ll b = P[0]; ll row = 0; for (auto p : P) { if (p == b) { row++; } else { rle.emplace_back(b, row); chmax(m, row); row = 1; b = p; } } rle.emplace_back(b, row); chmax(m, row); } // FPS X(m + 1, 1); // X = X.pow(s - 1); vec(mint, X, m + 1); mint t = 1; fori(i, m + 1) { X[i] = t; t *= (s - 1 + i); t /= i + 1; } mint ans = 1; for (auto [p, e] : rle) { mint np = 1; mint tot = 0; fori(i, e + 1) { tot += np * X[e - i]; np *= p; } ans *= tot; } ans *= s; print(ans); } int main() { #ifndef INTERACTIVE std::cin.tie(0)->sync_with_stdio(0); #endif // std::cout << std::fixed << std::setprecision(12); int t; t = 1; std::cin >> t; while (t--) solve(); return 0; } // // #pragma GCC target("avx2") // // #pragma GCC optimize("O3") // // #pragma GCC optimize("unroll-loops") // // #define INTERACTIVE // // #include "kyopro-cpp/template.hpp" // // #include "misc/Modint.hpp" // using mint = modint9; // #include "polynomial/FormalPowerSeries.hpp" // using FPS = FormalPowerSeries; // #include "math/pollard_rho.hpp" // #include "polynomial/BostanMori.hpp" // vector primes = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, // 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97}; // // void solve() { // LL(n, s); // if (n == 1) { // print(mint(s)); // return; // } // vec(ll, P, 0); // ll ss = n; // for (auto p : primes) { // while (ss % p == 0) { // P.push_back(p); // ss /= p; // } // } // // vec(Pll, rle, 0); // int m = 0; // { // ll b = P[0]; // ll row = 0; // for (auto p : P) { // if (p == b) { // row++; // } else { // rle.emplace_back(b, row); // chmax(m, row); // row = 1; // b = p; // } // } // rle.emplace_back(b, row); // chmax(m, row); // } // // // FPS X(m + 1, 1); // // X = X.pow(s - 1); // vec(mint, X, m + 1); // mint t = 1; // fori(i, m + 1) { // X[i] = t; // t *= (s - 1 + i); // t /= i + 1; // } // // mint ans = 1; // for (auto [p, e] : rle) { // mint np = 1; // mint tot = 0; // fori(i, e + 1) { // tot += np * X[e - i]; // np *= p; // } // ans *= tot; // } // // ans *= s; // print(ans); // } // // int main() { // #ifndef INTERACTIVE // std::cin.tie(0)->sync_with_stdio(0); // #endif // // std::cout << std::fixed << std::setprecision(12); // int t; // t = 1; // std::cin >> t; // while (t--) solve(); // return 0; // }