#define SINGLE_TESTCASE #define FAST_CIO #define INF 4'000'000'000'000'000'037LL #define EPS 1e-11 #include using namespace std; namespace { using ld = decltype(EPS); using ll = long long; using uint = unsigned int; using ull = unsigned long long; using pll = pair; using tlll = tuple; using tllll = tuple; #define vc vector template using vvc = vc>; using vpll = vc; using vstr = vc; #ifdef __SIZEOF_INT128__ using i128 = __int128_t; using u128 = __uint128_t; i128 stoi128(const string &s) { const bool neg = s.front() == '-'; u128 res = 0; for (int i = neg; i < (int)s.size(); i++) res = 10 * res + s[i] - '0'; if (neg) return -i128(res - 1) - 1; return i128(res); } string i128tos(i128 x) { if (x == 0) return "0"; string sign = "", res = ""; u128 ux; if (x < 0) ux = u128(-(x + 1)) + 1, sign = "-"; else ux = x; while (ux > 0) { res += '0' + ux % 10; ux /= 10; } reverse(res.begin(), res.end()); return sign + res; } istream &operator>>(istream &is, i128 &a) { string s; is >> s; a = stoi128(s); return is; } ostream &operator<<(ostream &os, const i128 &a) { os << i128tos(a); return os; } #endif #define cauto const auto #define overload4(_1, _2, _3, _4, name, ...) name #define rep1(i, n) for (ll i = 0, nnnnn = ll(n); i < nnnnn; i++) #define rep2(i, l, r) for (ll i = ll(l), rrrrr = ll(r); i < rrrrr; i++) #define rep3(i, l, r, d) for (ll i = ll(l), rrrrr = ll(r), ddddd = ll(d); ddddd > 0 ? i < rrrrr : i > rrrrr; i += d) #define rep(...) overload4(__VA_ARGS__, rep3, rep2, rep1)(__VA_ARGS__) #define repi1(i, n) for (int i = 0, nnnnn = int(n); i < nnnnn; i++) #define repi2(i, l, r) for (int i = int(l), rrrrr = int(r); i < rrrrr; i++) #define repi3(i, l, r, d) for (int i = int(l), rrrrr = int(r), ddddd = int(d); ddddd > 0 ? i < rrrrr : i > rrrrr; i += d) #define repi(...) overload4(__VA_ARGS__, repi3, repi2, repi1)(__VA_ARGS__) #define fe(...) for (auto __VA_ARGS__) #define fec(...) for (cauto &__VA_ARGS__) #define fem(...) for (auto &__VA_ARGS__) template constexpr bool is_integral_ext = is_integral_v || is_same_v || is_same_v; template constexpr bool is_signed_ext = is_signed_v || is_same_v; template constexpr bool is_unsigned_ext = is_unsigned_v || is_same_v; template constexpr T monoid_default_infty() { if constexpr (numeric_limits::is_specialized && numeric_limits::is_integer) return numeric_limits::max() / T(2); else return T(INF); } template constexpr T monoid_default_identity() { if constexpr (numeric_limits::is_specialized && numeric_limits::is_integer) return numeric_limits::max(); else return T(INF) + T(1); } template inline bool chmin(T &a, U b) { return a > b ? a = b, true : false; } template inline bool chmax(T &a, U b) { return a < b ? a = b, true : false; } template && is_integral_ext>> inline constexpr T divfloor(U a, V b) { return T(a) / T(b) - (T(a) % T(b) && (T(a) ^ T(b)) < 0); } template && is_integral_ext>> inline constexpr T divceil(U a, V b) { return T(a) / T(b) + (T(a) % T(b) && (T(a) ^ T(b)) >= 0); } template && is_integral_ext>> inline constexpr T divround(U a, V b) { return divfloor(2 * T(a) + T(b), 2 * T(b)); } template && is_integral_ext>> inline constexpr T safemod(U a, V b) { return T(a) - T(b) * divfloor(a, b); } template constexpr T ipow(U a, V b) { assert(b >= 0); if (b == 0) return 1; if (a == 0 || a == 1) return a; if (a < 0 && a == -1) return b & 1 ? -1 : 1; T res = 1, tmp = a; while (true) { if (b & 1) res *= tmp; b >>= 1; if (b == 0) break; tmp *= tmp; } return res; } template T mul_limited(A a, B b, M m) { assert(a >= 0 && b >= 0 && m >= 0); if (b == 0) return 0; return T(a) > T(m) / T(b) ? T(m) : T(a) * T(b); } template T mul_limited(A a, B b) { return mul_limited(a, b, INF); } template T pow_limited(A a, B b, M m) { assert(a >= 0 && b >= 0 && m >= 0); if (a <= 1 || b == 0) return min(ipow(a, b), T(m)); T res = 1, tmp = a; while (true) { if (b & 1) { if (res > T(m) / tmp) return m; res *= tmp; } b >>= 1; if (b == 0) break; if (tmp > T(m) / tmp) return m; tmp *= tmp; } return res; } template T pow_limited(A a, B b) { return pow_limited(a, b, INF); } template constexpr T iroot(A a, K k) { assert(a >= 0 && k >= 1); if (a <= 1 || k == 1) return a; if (k == 2) { const T aa = T(a); T x = T(sqrtl((long double)a)); while (x > aa / x) x--; while (x < numeric_limits::max()) { const T y = x + 1; if (y > aa / y) break; x = y; } return x; } auto isok = [&](T x) -> bool { if (x == 0) return true; T res = 1, k2 = k; while (true) { if (k2 & 1) { if (res > T(a) / x) return false; res *= x; } k2 >>= 1; if (k2 == 0) break; if (x > T(a) / x) return false; x *= x; } return res <= T(a); }; T x = pow(a, 1.0 / k); bool up = true; while (!isok(x)) up = false, x--; if (up) { while (x < numeric_limits::max() && isok(x + 1)) x++; } return x; } template constexpr T iroot_ceil(A a, K k) { T x = iroot(a, k); return ipow(x, k) == a ? x : x + 1; } template int SGN(A a, D eps = EPS) { return int(a > eps) - int(a < -eps); } template vc base_repr(U val, V base) { assert(val >= 0); assert(base >= 2); if (val == 0) return {0}; vc a; while (val > 0) { a.emplace_back(val % base); val /= base; } reverse(a.begin(), a.end()); return a; } template vc base_repr(U val, V base, int n) { assert(val >= 0); assert(base >= 2); assert(n >= 0); vc a(n); repi(i, n) { a[i] = val % base; val /= base; } reverse(a.begin(), a.end()); return a; } template string base_repr_str(U val, int base) { assert(val >= 0); assert(2 <= base && base <= 36); auto a = base_repr(val, base); string s = ""; for (cauto &ai : a) s += (ai < 10 ? '0' + ai : (use_upper ? 'A' : 'a') + (ai - 10)); return s; } template string base_repr_str(U val, int base, int n) { assert(val >= 0); assert(2 <= base && base <= 36); assert(n >= 0); auto a = base_repr(val, base, n); string s = ""; for (cauto &ai : a) s += (ai < 10 ? '0' + ai : (use_upper ? 'A' : 'a') + (ai - 10)); return s; } #define ALL(a) (a).begin(), (a).end() template inline T SZ(const V &x) { return x.size(); } #define eb emplace_back #define LMD(x, fx) ([&](auto x) { return fx; }) template auto gen_vec(int n, const F &f) { vc res(n); repi(i, n) res[i] = f(i); return res; } #define GEN_VEC(n, i, fi) (gen_vec(n, LMD(i, fi))) template auto dvec(const V (&sz)[d], const T &init) { if constexpr (i < d) return vc(sz[i], dvec(sz, init)); else return init; } template T ctol(const char &c, const string &s) { repi(i, SZ(s)) if (s[i] == c) return i; return -1; } template vc stov(const string &s, char first) { return gen_vec(SZ(s), [&](int i) -> T { return s[i] - first; }); } template vc stov(const string &s, const string &t) { return gen_vec(SZ(s), [&](int i) -> T { return ctol(s[i], t); }); } template string vtos(const vc &v, char first) { string res = ""; fe(vi : v) res += vi + first; return res; } template string vtos(const vc &v, const string &t) { string res = ""; fe(vi : v) res += t[vi]; return res; } template vc concat(const vvc &vs) { vc res; for (cauto &v : vs) res.insert(res.end(), ALL(v)); return res; } template vc concat(const vc &v) { return v; } template vc concat(vc v, const vc &...vs) { (v.insert(v.end(), ALL(vs)), ...); return v; } template vc merged(const vc &a, const vc &b) { vc res; merge(ALL(a), ALL(b), back_inserter(res)); return res; } template T vecget(const vc &v, I i, const T &dflt_negative = -monoid_default_infty(), const T &dflt_positive = monoid_default_infty()) { if (i < 0) return dflt_negative; if (i >= SZ(v)) return dflt_positive; return v[i]; } template auto SUM(const V &v) { typename V::value_type s{}; fec(vi : v) s += vi; return s; } template T SUM(const V &v) { T s{}; fec(vi : v) s += vi; return s; } template auto MAX(const V &v) { return *max_element(ALL(v)); } template auto MIN(const V &v) { return *min_element(ALL(v)); } template I ARGMAX(const V &v) { return max_element(ALL(v)) - v.begin(); } template I ARGMIN(const V &v) { return min_element(ALL(v)) - v.begin(); } template T mex(const V &a) { int n = a.size(); vector exists(n, false); repi(i, n) if (0 <= a[i] && a[i] < n) exists[a[i]] = true; repi(x, n) if (!exists[x]) return x; return n; } template bool is_permutation(const vc &p) { const int n = p.size(); vc b(n, false); repi(i, n) { if (!(0 <= p[i] && p[i] < n)) return false; b[p[i]] = true; } return all_of(ALL(b), [](bool bi) { return bi; }); } template vc permid(const int &n, const int &base_index = 0) { vc p(n); repi(i, n) p[i] = i + base_index; return p; } template vc perminv(const vc &p) { if (p.empty()) return {}; const int n = p.size(); vc q(MAX(p) + 1); repi(i, n) if (p[i] >= 0) q[p[i]] = i; return q; } template vc permuted(const vc &a, const vc &p) { const int n = p.size(); vc res(n); repi(i, n) { assert(0 <= p[i] && p[i] < U(a.size())); res[i] = a[p[i]]; } return res; } template vc permuted(const vc &p, const vc &q, const vc &...rs) { return permuted(permuted(p, q), rs...); } template V reversed(const V &v) { return V(v.rbegin(), v.rend()); } #if __cplusplus < 202002L template V sorted(V v, Args&&... args) { sort(ALL(v), forward(args)...); return v; } #else template V sorted(V v, Args&&... args) { ranges::sort(v, forward(args)...); return v; } #endif template void unique(V &v) { v.erase(std::unique(ALL(v)), v.end()); } template V uniqued(V v) { unique(v); return v; } template void sortunique(V &v) { sort(ALL(v)); unique(v); } template V sortuniqued(V v) { sortunique(v); return v; } template void rotate(V &v, U k) { const U n = v.size(); if (n == 0) return; k = (k % n + n) % n; std::rotate(v.begin(), v.begin() + k, v.end()); } template V rotated(V v, U k) { rotate(v, k); return v; } template vvc top(const vvc &a) { if (a.empty()) return {}; const int n = a.size(), m = a[0].size(); vvc b(m, vc(n)); repi(i, n) { assert(SZ(a[i]) == m); repi(j, m) b[j][i] = a[i][j]; } return b; } vstr top(const vstr &a) { vvc a_(a.size()); repi(i, SZ(a)) a_[i] = {ALL(a[i])}; vvc b_ = top(a_); vstr b(b_.size()); repi(i, SZ(b)) b[i] = {ALL(b_[i])}; return b; } template struct has_e0 : false_type {}; template struct has_e0> : true_type {}; template inline constexpr bool has_e0_v = has_e0::value; template struct MonoidAdd { using S = T; static constexpr S op(S a, S b) { return a + b; } static constexpr S e() { if constexpr (has_e0_v) return S::e0(); else return {}; } }; template ()> struct MonoidMin { using S = T; static constexpr S op(S a, S b) { return min(a, b); } static constexpr S e() { return infty; } }; template ()> struct MonoidMax { using S = T; static constexpr S op(S a, S b) { return max(a, b); } static constexpr S e() { return -infty; } }; template typename M::S pow_monoid(typename M::S a, ll k) { typename M::S c = M::e(); for (; k; k >>= 1) { if (k & 1) c = M::op(c, a); a = M::op(a, a); } return c; } template vc cuml(const vc &v, int left_index = 0) { const int n = v.size(); vc res(n + 1); res[0] = M::e(); repi(i, n) res[i + 1] = M::op(res[i], v[i]); res.erase(res.begin(), res.begin() + left_index); return res; } template vc cumr(const vc &v, int right_index = 0) { return reversed(cuml(reversed(v), right_index)); } template vc cumlsum(const vc &v, int left_index = 0) { return cuml>(v, left_index); } template vc cumrsum(const vc &v, int right_index = 0) { return cumr>(v, right_index); } template vc cumlmin(const vc &v, int left_index = 0) { return cuml>(v, left_index); } template vc cumrmin(const vc &v, int right_index = 0) { return cumr>(v, right_index); } template vc cumlmax(const vc &v, int left_index = 0) { return cuml>(v, left_index); } template vc cumrmax(const vc &v, int right_index = 0) { return cumr>(v, right_index); } template vc adjd(const vc &v, int left_index = 0, int right_index = 0) { int n = v.size(); assert(0 <= left_index && 0 <= right_index && left_index + right_index <= n + 1); vc res(n + 1); if (n == 0) { res[0] = T{}; res.erase(res.end() - right_index, res.end()); res.erase(res.begin(), res.begin() + left_index); return res; } res[0] = v[0]; repi(i, 1, n) res[i] = v[i] - v[i - 1]; res[n] = -v[n - 1]; res.erase(res.end() - right_index, res.end()); res.erase(res.begin(), res.begin() + left_index); return res; } const vpll DRULgrid = {{1, 0}, {0, 1}, {-1, 0}, {0, -1}}; const vpll DRULplane = {{0, -1}, {1, 0}, {0, 1}, {-1, 0}}; template struct is_random_access_iterator { static constexpr bool value = is_same_v< typename iterator_traits::iterator_category, random_access_iterator_tag >; }; template constexpr bool is_random_access_iterator_v = is_random_access_iterator::value; #if __cplusplus < 202002L struct identity { template constexpr T &&operator()(T &&t) const noexcept { return forward(t); } }; namespace internal { template inline T bound_helper(const V &v, Judge judge) { int l = -1, r = v.size(); while (r - l > 1) { int m = (l + r) / 2; if (judge(m)) l = m; else r = m; } return r; } }; template , class Proj = identity> inline T LB(const V &v, const Value &val, Comp comp = {}, Proj proj = {}) { return internal::bound_helper(v, [&](int i) -> bool { return comp(proj(*(v.begin() + i)), val); }); } template , class Proj = identity> inline T UB(const V &v, const Value &val, Comp comp = {}, Proj proj = {}) { return internal::bound_helper(v, [&](int i) -> bool { return !comp(val, proj(*(v.begin() + i))); }); } #define DEFAULT_COMP less<> #else template inline T LB(const V &v, const Value &val, Comp comp = {}, Proj proj = {}) { return ranges::lower_bound(v, val, comp, proj) - v.begin(); } template inline T UB(const V &v, const Value &val, Comp comp = {}, Proj proj = {}) { return ranges::upper_bound(v, val, comp, proj) - v.begin(); } #define DEFAULT_COMP ranges::less #endif template inline auto lt_max(const V &v, const Value &val, Comp comp = {}, Proj proj = {}) -> enable_if_t, T> { return LB(v, val, comp, proj) - 1; } template inline auto leq_max(const V &v, const Value &val, Comp comp = {}, Proj proj = {}) -> enable_if_t, T> { return UB(v, val, comp, proj) - 1; } template inline auto gt_min(const V &v, const Value &val, Comp comp = {}, Proj proj = {}) -> enable_if_t, T> { return UB(v, val, comp, proj); } template inline auto geq_min(const V &v, const Value &val, Comp comp = {}, Proj proj = {}) -> enable_if_t, T> { return LB(v, val, comp, proj); } template inline auto lt_cnt(const V &v, const Value &val, Comp comp = {}, Proj proj = {}) -> enable_if_t, T> { return LB(v, val, comp, proj); } template inline auto leq_cnt(const V &v, const Value &val, Comp comp = {}, Proj proj = {}) -> enable_if_t, T> { return UB(v, val, comp, proj); } template inline auto gt_cnt(const V &v, const Value &val, Comp comp = {}, Proj proj = {}) -> enable_if_t, T> { return SZ(v) - UB(v, val, comp, proj); } template inline auto geq_cnt(const V &v, const Value &val, Comp comp = {}, Proj proj = {}) -> enable_if_t, T> { return SZ(v) - LB(v, val, comp, proj); } template inline auto in_cnt(const V &v, L l, R r, Comp comp = {}, Proj proj = {}) -> enable_if_t, T> { if (l > r) return 0; return lt_cnt(v, r, comp, proj) - lt_cnt(v, l, comp, proj); } template inline auto lt_max(const V &v, const Value &val) -> enable_if_t, typename V::const_iterator> { auto it = v.lower_bound(val); return it == v.begin() ? v.end() : prev(it); } template inline auto leq_max(const V &v, const Value &val) -> enable_if_t, typename V::const_iterator> { auto it = v.upper_bound(val); return it == v.begin() ? v.end() : prev(it); } template inline auto gt_min(const V &v, const Value &val) -> enable_if_t, typename V::const_iterator> { return v.upper_bound(val); } template inline auto geq_min(const V &v, const Value &val) -> enable_if_t, typename V::const_iterator> { return v.lower_bound(val); } namespace internal { template bool binsearch_adjacent(T a, T b) { if (a < b) return a + 1 == b; if (b < a) return b + 1 == a; return false; } }; template pair binsearch(const Judge &judge, InitOk init_ok, InitNg init_ng, bool check_ok = true, bool check_ng = true) { T ok(init_ok), ng(init_ng); if (check_ok) assert(judge(ok)); if (check_ng) assert(!judge(ng)); while (!internal::binsearch_adjacent(ok, ng)) { T mid = (ok & ng) + ((ok ^ ng) >> 1); (judge(mid) ? ok : ng) = mid; } return {ok, ng}; } template T binsearch_real(const Judge &judge, InitOk init_ok, InitNg init_ng, int iteration_count = 100, bool check_ok = true, bool check_ng = true) { T ok(init_ok), ng(init_ng); if (check_ok) assert(judge(ok)); if (check_ng) assert(!judge(ng)); repi(_, iteration_count) { T mid = (ok + ng) / 2; (judge(mid) ? ok : ng) = mid; } return ok; } template pair expsearch(const Judge &judge, InitVal init_val, bool positive = true) { T cur(init_val), step = 1; const bool cur_ok = judge(cur); while (true) { T nxt; if (positive) nxt = cur > numeric_limits::max() - step ? numeric_limits::max() : cur + step; else nxt = cur < numeric_limits::lowest() + step ? numeric_limits::lowest() : cur - step; assert(nxt != cur && "the boundary must exist in the searched direction"); if (nxt == cur || judge(nxt) != cur_ok) { T ok = cur_ok ? cur : nxt; T ng = cur_ok ? nxt : cur; return binsearch(judge, ok, ng, false, false); } cur = nxt; if (step > numeric_limits::max() / 2) step = numeric_limits::max(); else step *= 2; } } template inline constexpr ull pow2(T k) { return 1ULL << k; } template inline constexpr ull MASK(T k) { return (1ULL << k) - 1ULL; } #if __cplusplus < 202002L inline constexpr ull bit_width(ull x) { return x == 0 ? 0 : 64 - __builtin_clzll(x); } inline constexpr ull bit_floor(ull x) { return x == 0 ? 0ULL : 1ULL << (bit_width(x) - 1); } inline constexpr ull bit_ceil(ull x) { return x == 0 ? 1ULL : 1ULL << bit_width(x - 1); } inline constexpr ull countr_zero(ull x) { assert(x != 0); return __builtin_ctzll(x); } inline constexpr ull popcount(ull x) { return __builtin_popcountll(x); } inline constexpr bool has_single_bit(ull x) { return popcount(x) == 1; } #else inline constexpr ll bit_width(ll x) { return std::bit_width((ull)x); } inline constexpr ll bit_floor(ll x) { return std::bit_floor((ull)x); } inline constexpr ll bit_ceil(ll x) { return std::bit_ceil((ull)x); } inline constexpr ll countr_zero(ll x) { assert(x != 0); return std::countr_zero((ull)x); } inline constexpr ll popcount(ll x) { return std::popcount((ull)x); } inline constexpr bool has_single_bit(ll x) { return std::has_single_bit((ull)x); } #endif inline constexpr ull lsb_pos(ull x) { assert(x != 0); return countr_zero(x); } inline constexpr ull msb_pos(ull x) { assert(x != 0); return bit_width(x) - 1; } inline constexpr ull lsb_mask(ull x) { assert(x != 0); return x & -x; } inline constexpr ull msb_mask(ull x) { assert(x != 0); return bit_floor(x); } inline constexpr bool btest(ull x, uint k) { return (x >> k) & 1; } template inline void bset(T &x, uint k, bool b = 1) { b ? x |= (1ULL << k) : x &= ~(1ULL << k); } template inline void bflip(T &x, uint k) { x ^= (1ULL << k); } inline constexpr bool bsubset(ull x, ull y) { return (x & y) == x; } inline constexpr bool bsupset(ull x, ull y) { return (x & y) == y; } inline constexpr ull bsetminus(ull x, ull y) { return x & ~y; } #define dump(...) #define local(...) #define oj(...) __VA_ARGS__ #define local_oj(a, b) (b) template vc content(queue que) { vc res; while (!que.empty()) { res.eb(que.front()); que.pop(); } return res; } template vc content(priority_queue pque) { vc res; while (!pque.empty()) { res.eb(pque.top()); pque.pop(); } return res; } template auto content(const T &obj) { return obj.content(); } template istream &operator>>(istream &is, pair &p) { is >> p.first >> p.second; return is; } template istream &operator>>(istream &is, tuple &t) { apply([&](auto &...a) { (is >> ... >> a); }, t); return is; } template istream &operator>>(istream &is, array &a) { for (size_t i = 0; i < n; i++) is >> a[i]; return is; } template istream &operator>>(istream &is, vc &a) { const size_t n = a.size(); for (size_t i = 0; i < n; i++) is >> a[i]; return is; } namespace internal { template void CIN(Ts &...a) { (cin >> ... >> a); } template void READnodump(Ts &...a) { CIN(a...); } template void READVECnodump(int n, vc &v) { v.resize(n); READnodump(v); } template void READVECnodump(int n, vc &v, vc &...vs) { READVECnodump(n, v), READVECnodump(n, vs...); } template void READVEC2nodump(int n, int m, vvc &v) { v.assign(n, vc(m)); READnodump(v); } template void READVEC2nodump(int n, int m, vvc &v, vvc &...vs) { READVEC2nodump(n, m, v), READVEC2nodump(n, m, vs...); } template void READJAGnodump(int n, vvc &v) { v.resize(n); repi(i, n) { int k; READnodump(k); READVECnodump(k, v[i]); } } template void READJAGnodump(int n, vvc &v, vvc &...vs) { READJAGnodump(n, v), READJAGnodump(n, vs...); } }; // namespace internal #define READ(...) internal::READnodump(__VA_ARGS__); dump(__VA_ARGS__) #define IN(T, ...) T __VA_ARGS__; READ(__VA_ARGS__) #define CHAR(...) IN(char, __VA_ARGS__) #define INT(...) IN(int, __VA_ARGS__) #define LL(...) IN(ll, __VA_ARGS__) #define STR(...) IN(string, __VA_ARGS__) #define ARR(T, n, ...) array __VA_ARGS__; READ(__VA_ARGS__) #define READVEC(...) internal::READVECnodump(__VA_ARGS__); dump(__VA_ARGS__) #define READVEC2(...) internal::READVEC2nodump(__VA_ARGS__); dump(__VA_ARGS__) #define VEC(T, n, ...) vc __VA_ARGS__; READVEC(n, __VA_ARGS__) #define VEC2(T, n, m, ...) vvc __VA_ARGS__; READVEC2(n, m, __VA_ARGS__) #define READJAG(...) internal::READJAGnodump(__VA_ARGS__); dump(__VA_ARGS__) #define JAG(T, n, ...) vvc __VA_ARGS__; READJAG(n, __VA_ARGS__) #define ENDL '\n' template ostream &operator<<(ostream &os, const pair &p) { os << p.first << ' ' << p.second; return os; } namespace internal { template void cout_tuple(ostream &os, const T &t) { if constexpr (N < std::tuple_size::value) { if constexpr (N > 0) { os << ' '; } const auto x = std::get(t); os << x; cout_tuple(os, t); } } }; // namespace internal template ostream &operator<<(ostream &os, const tuple &t) { internal::cout_tuple(os, t); return os; } template ostream &operator<<(ostream &os, const array &a) { for (size_t i = 0; i < n; i++) { if (i) os << ' '; os << a[i]; } return os; } template ostream &operator<<(ostream &os, const vc &v) { const size_t n = v.size(); for (size_t i = 0; i < n; i++) { if (i) os << ' '; os << v[i]; } return os; } namespace internal { template void COUTW() {} template void COUTW(const Ts &...a) { (cout << ... << a); } inline void COUTP() { cout << ENDL; } template void COUTP(const T &a) { cout << a << ENDL; } template void COUTP(const T &a, const Ts &...b) { cout << a; (cout << ... << (cout << ' ', b)); cout << ENDL; } }; // namespace internal #define WRITE internal::COUTW #define PRINT internal::COUTP #define PRINTEXIT(...) do { PRINT(__VA_ARGS__); exit(0); } while (false) #define PRINTRETURN(...) do { PRINT(__VA_ARGS__); return; } while (false) template void PRINTV(const vc &v) { for (auto &vi : v) PRINT(vi); } #define PRINTVEXIT(...) do { PRINTV(__VA_ARGS__); exit(0); } while (false) #define PRINTVRETURN(...) do { PRINTV(__VA_ARGS__); return; } while (false) template pair &operator+=(pair &a, const P &b) { a.first += b.first; a.second += b.second; return a; } template pair operator+(pair a, const P &b) { return a += b; } template array &operator+=(array &a, const A &b) { for (size_t i = 0; i < n; i++) a[i] += b[i]; return a; } template array operator+(array a, const A &b) { return a += b; } namespace internal { template auto &tuple_add_impl(A &a, const B &b, const index_sequence) { ((get(a) += get(b)), ...); return a; } }; // namespace internal template tuple &operator+=(tuple &a, const Tp &b) { return internal::tuple_add_impl(a, b, make_index_sequence>>{}); } template tuple operator+(tuple a, const Tp &b) { return a += b; } template void offset(vc &v, const Add &add) { for (auto &vi : v) vi += add; } template void offset(vvc &v, const Add &add) { for (auto &vi : v) for (auto &vij : vi) vij += add; } template array, m> unzip(const vc> &vt) { const size_t n = vt.size(); array, m> tv; tv.fill(vc(n)); for (size_t i = 0; i < n; i++) for (size_t j = 0; j < m; j++) tv[j][i] = vt[i][j]; return tv; } template vc> zip(const array, m> &tv) { if (tv.empty()) return {}; const size_t n = tv[0].size(); vc> vt(n); for (size_t j = 0; j < m; j++) { assert(tv[j].size() == n); for (size_t i = 0; i < n; i++) vt[i][j] = tv[j][i]; } return vt; } template pair, vc> unzip(const vc> &vt) { const size_t n = vt.size(); pair, vc> tv; tv.first.resize(n), tv.second.resize(n); for (size_t i = 0; i < n; i++) tie(tv.first[i], tv.second[i]) = vt[i]; return tv; } template vc> zip(const pair, vc> &tv) { const size_t n = tv.first.size(); assert(n == tv.second.size()); vc> vt(n); for (size_t i = 0; i < n; i++) vt[i] = make_pair(tv.first[i], tv.second[i]); return vt; } namespace internal { template auto vt_to_tv_impl(V &tv, const Tp &t, index_sequence, size_t index) { ((get(tv)[index] = get(t)), ...); } template auto tv_to_vt_impl(const Tp &tv, index_sequence, size_t index) { return make_tuple(get(tv)[index]...); } }; template auto unzip(const vc> &vt) { const size_t n = vt.size(); tuple...> tv; apply([&](auto &...v) { ((v.resize(n)), ...); }, tv); for (size_t i = 0; i < n; i++) internal::vt_to_tv_impl(tv, vt[i], make_index_sequence>{}, i); return tv; } template auto zip(const tuple...> &tv) { size_t n = get<0>(tv).size(); apply([&](auto &...v) { ((void(v), assert(v.size() == n)), ...); }, tv); vc> vt(n); for (size_t i = 0; i < n; i++) vt[i] = internal::tv_to_vt_impl(tv, index_sequence_for{}, i); return vt; } #define UNZIP(vt, ...) auto [__VA_ARGS__] = unzip(vt) #define ZIP(vt, ...) auto vt = zip(tuple{__VA_ARGS__}) mt19937_64 mt; template T randint(U1 l, U2 r) { assert(T(l) <= T(r)); return T(l) + mt() % (T(r) - T(l) + 1); } template T randrange(U1 l, U2 r) { assert(T(l) < T(r)); return T(l) + mt() % (T(r) - T(l)); } template array random_sample_range_array(U1 l, U2 r) { assert(T(r) - T(l) >= T(k)); array res; repi(i, k) res[i] = randint(T(l), T(r) - T(k)); sort(ALL(res)); repi(i, k) res[i] += i; if (!does_sort) shuffle(ALL(res), mt); return res; } template vc random_sample_range_vector(U1 l, U2 r, int k) { assert(T(r) - T(l) >= T(k)); vc res(k); repi(i, k) res[i] = randint(T(l), T(r) - T(k)); sort(ALL(res)); repi(i, k) res[i] += i; if (!does_sort) shuffle(ALL(res), mt); return res; } template struct larger_int { using type = T; }; #define LARGER_INT(T, U) \ template <> \ struct larger_int \ { \ using type = U; \ }; LARGER_INT(signed char, short) LARGER_INT(short, int) LARGER_INT(int, long long) LARGER_INT(long, __int128_t) LARGER_INT(long long, __int128_t) LARGER_INT(unsigned char, unsigned short) LARGER_INT(unsigned short, unsigned int) LARGER_INT(unsigned int, unsigned long long) LARGER_INT(unsigned long, __uint128_t) LARGER_INT(unsigned long long, __uint128_t) #undef LARGER_INT template using larger_int_t = typename larger_int::type; namespace internal { template constexpr ll powmod_constexpr(ll x, ll n, T m) { if (m == 1) return 0; using U = make_unsigned_t; using L = larger_int_t; U r = 1, y = safemod(x, m); while (n) { if (n & 1) r = L(r) * y % m; y = L(y) * y % m; n >>= 1; } return r; } template constexpr bool isprime_constexpr(T n) { if constexpr (sizeof(T) > 4) { if (n <= INT_MAX) return isprime_constexpr(n); } if (n <= 1) return false; if (n == 2 || n == 7 || n == 61) return true; if (n % 2 == 0) return false; ll d = n - 1; while (d % 2 == 0) d /= 2; using U = make_unsigned_t; using L = larger_int_t; auto miller_rabin = [&](const auto &bases) constexpr { for (ll a : bases) { ll t = d, y = powmod_constexpr(a, t, n); while (t != n - 1 && y != 1 && y != n - 1) { y = L(y) * y % n; t <<= 1; } if (y != n - 1 && t % 2 == 0) return false; } return true; }; if constexpr(sizeof(T) <= 4) { constexpr ll bases[3] = {2, 7, 61}; return miller_rabin(bases); } else { constexpr ll bases[7] = {2, 325, 9375, 28178, 450775, 9780504, 1795265022}; return miller_rabin(bases); } } template constexpr bool isprime = isprime_constexpr(n); }; namespace internal { template struct policy_static { using mod_type = decltype(M); using value_type = make_unsigned_t; using calc_type = larger_int_t; static constexpr bool is_prime = isprime_constexpr(M); static constexpr mod_type mod() { return M; } static constexpr value_type umod() { return M; } static constexpr value_type init(value_type v) { return v; } static constexpr mod_type val(value_type v) { return v; } static constexpr value_type mul(value_type a, value_type b) { return (value_type)((calc_type(a) * b) % M); } }; }; namespace internal { struct barrett32 { uint m; ull im; explicit barrett32(uint m) : m(m), im((ull)(-1) / m + 1) {} uint umod() const { return m; } uint mul(uint a, uint b) const { ull z = a; z *= b; ull x = ull((u128(z) * im) >> 64); ull y = x * m; return uint(z - y + (z < y ? m : 0)); } }; template struct policy_barrett32 { using value_type = uint; using calc_type = ull; using mod_type = int; static constexpr bool is_prime = false; static inline barrett32 reducer{998244353}; static void set_mod(mod_type m) { reducer = barrett32(m); } static mod_type mod() { return reducer.umod(); } static value_type umod() { return reducer.umod(); } static value_type init(value_type v) { return v; } static mod_type val(value_type v) { return v; } static value_type mul(value_type a, value_type b) { return reducer.mul(a, b); } }; }; namespace internal { inline constexpr ull inv64(ull a) { ull x = a; while (a * x != 1) x *= 2 - a * x; return x; } struct montgomery64odd { ull m, im, sq; explicit montgomery64odd(ull m) : m(m), im(inv64(m)), sq(-u128(m) % m) {} ull umod() const { return m; } ull reduce(u128 x) const { auto t = (x + u128(m) * (-im * ull(x))) >> 64; if (t >= m) t -= m; return (ull)t; } ull inv_reduce(i128 v) const { return reduce(u128(v % m + m) * sq); } }; struct montgomery64 { ull m, mx, imx, d, q; uint b; explicit montgomery64(ull m) : m(m) { b = countr_zero(m), mx = m >> b; imx = inv64(mx); d = powmod_constexpr((mx + 1) / 2, b, mx); u128 sq = -u128(mx) % mx; q = (1 + (((sq - 1) * d) << b)) % m; } ull umod() const { return m; } ull reduce(u128 x) const { if (b == 0) { auto t = (x + u128(mx) * (-imx * ull(x))) >> 64; if (t >= m) t -= m; return (ull)t; } ull p = x & MASK(b); x = (x >> b) + p * d; ull y = p << (64 - b); auto t = (x + u128(mx) * (imx * (y - ull(x)))) >> (64 - b); if (t >= m) { t -= m; if (t >= m) t -= m; } return (ull)t; } ull inv_reduce(i128 v) const { return reduce(u128(v % m + m) * q); } }; template struct policy_montgomery64_odd { using value_type = ull; using calc_type = u128; using mod_type = ll; static constexpr bool is_prime = false; static inline montgomery64odd reducer{(1LL << 61) - 1}; static void set_mod(mod_type m) { reducer = montgomery64odd(m); } static mod_type mod() { return reducer.umod(); } static value_type umod() { return reducer.umod(); } static value_type init(value_type v) { return reducer.inv_reduce(v); } static mod_type val(value_type v) { return reducer.reduce(v); } static value_type mul(value_type a, value_type b) { return reducer.reduce((calc_type)a * b); } }; template struct policy_montgomery64 { using value_type = ull; using calc_type = u128; using mod_type = ll; static constexpr bool is_prime = false; static inline montgomery64 reducer{(1LL << 61) - 1}; static void set_mod(mod_type m) { reducer = montgomery64(m); } static mod_type mod() { return reducer.umod(); } static value_type umod() { return reducer.umod(); } static value_type init(value_type v) { return reducer.inv_reduce(v); } static mod_type val(value_type v) { return reducer.reduce(v); } static value_type mul(value_type a, value_type b) { return reducer.reduce((calc_type)a * b); } }; }; template constexpr tuple extgcd(T a, T b) { if (a == 0 && b == 0) return {0, 0, 0}; T x1 = 1, y1 = 0, z1 = a; T x2 = 0, y2 = 1, z2 = b; while (z2 != 0) { T q = z1 / z2; tie(x1, x2) = make_pair(x2, x1 - q * x2); tie(y1, y2) = make_pair(y2, y1 - q * y2); tie(z1, z2) = make_pair(z2, z1 - q * z2); } if (z1 < 0) x1 = -x1, y1 = -y1, z1 = -z1; return {z1, x1, y1}; } namespace internal { template struct modint_impl { using V = typename Policy::value_type; using M = typename Policy::mod_type; using mint = modint_impl; private: V _v; public: static constexpr M mod() { return Policy::mod(); } template static auto set_mod(M m) -> decltype(T::set_mod(m)) { return T::set_mod(m); } static mint raw(V v) { mint x; x._v = v; return x; } modint_impl() : _v(0) {} template >> modint_impl(T v) { V rem; if constexpr (is_signed_ext) { using S = make_signed_t; S x = v % S(Policy::umod()); if (x < 0) x += Policy::umod(); rem = x; } else rem = V(v % Policy::umod()); _v = Policy::init(rem); }; M val() const { return Policy::val(_v); } mint &operator+=(const mint &rhs) { _v += rhs._v; if (_v >= Policy::umod()) _v -= Policy::umod(); return *this; } mint &operator-=(const mint &rhs) { _v -= rhs._v; if (_v >= Policy::umod()) _v += Policy::umod(); return *this; } mint &operator*=(const mint &rhs) { _v = Policy::mul(_v, rhs._v); return *this; } mint &operator/=(const mint &rhs) { return *this *= rhs.inv(); } mint &operator++() { _v++; if (_v == Policy::umod()) _v = 0; return *this; } mint &operator--() { if (_v == 0) _v = Policy::umod(); _v--; return *this; } mint operator++(int) { mint res = *this; ++(*this); return res; } mint operator--(int) { mint res = *this; --(*this); return res; } mint operator+() const { return *this; } mint operator-() const { return mint() - *this; } template mint pow(T n) const { assert(n >= 0); mint x = *this, r = 1; while (n) { if (n & 1) r *= x; x *= x; n >>= 1; } return r; } mint inv() const { if constexpr (Policy::is_prime) { return pow(mod() - 2); } else { auto [g, x, y] = extgcd(val(), mod()); assert(g == 1); return mint(x); } } friend mint operator+(const mint &lhs, const mint &rhs) { return mint(lhs) += rhs; } friend mint operator-(const mint &lhs, const mint &rhs) { return mint(lhs) -= rhs; } friend mint operator*(const mint &lhs, const mint &rhs) { return mint(lhs) *= rhs; } friend mint operator/(const mint &lhs, const mint &rhs) { return mint(lhs) /= rhs; } friend bool operator==(const mint &lhs, const mint &rhs) { return lhs._v == rhs._v; } friend bool operator!=(const mint &lhs, const mint &rhs) { return lhs._v != rhs._v; } friend std::istream &operator>>(std::istream &is, mint &x) { long long a; is >> a; x = a; return is; } friend std::ostream &operator<<(std::ostream &os, const mint &x) { os << x.val(); return os; } }; }; template using static_modint32 = internal::modint_impl>; template using dynamic_modint32 = internal::modint_impl>; template using static_modint64 = internal::modint_impl>; template using dynamic_modint64_odd = internal::modint_impl>; template using dynamic_modint64 = internal::modint_impl>; using modint998244353 = static_modint32<998244353>; template struct is_modint : std::false_type { }; template struct is_modint> : std::true_type { }; template inline constexpr bool is_modint_v = is_modint::value; template struct is_static_modint : false_type {}; template struct is_static_modint> : true_type {}; template struct is_static_modint> : true_type {}; template inline constexpr bool is_static_modint_v = is_static_modint::value; template struct is_dynamic_modint : false_type {}; template struct is_dynamic_modint> : true_type {}; template struct is_dynamic_modint> : true_type {}; template struct is_dynamic_modint> : true_type {}; template inline constexpr bool is_dynamic_modint_v = is_dynamic_modint::value; template struct has_mod : std::false_type { }; template struct has_mod> : std::true_type { }; template struct PowerTable { private: decltype(mint::mod()) mod; mint base; vc pw; public: PowerTable() {} PowerTable(const mint &base) : mod(mint::mod()), base(base), pw(1, 1) {} void reserve(int n) { if (mod != mint::mod()) { mod = mint::mod(); pw = {1}; } int i = pw.size(); if (n < i) return; pw.resize(n + 1); for (; i <= n; i++) pw[i] = pw[i - 1] * base; } mint pow(int n) { reserve(n); return pw[n]; } }; template struct Binomial { private: inline static decltype(T::mod()) mod; public: inline static vc fac_, finv_, inv_; static void reserve(int n) { if constexpr (is_dynamic_modint_v) { if (mod != T::mod()) { mod = T::mod(); fac_ = {1, 1}, finv_ = {1, 1}, inv_ = {0, 1}; } } else { if (fac_.empty()) fac_ = {1, 1}, finv_ = {1, 1}, inv_ = {0, 1}; } if (n < SZ(fac_)) return; chmin(n, T::mod() - 1); int si = fac_.size(); fac_.resize(n + 1), finv_.resize(n + 1), inv_.resize(n + 1); repi(i, si, n + 1) { fac_[i] = fac_[i - 1] * T::raw(i); inv_[i] = -inv_[T::mod() % i] * T::raw(T::mod() / i); finv_[i] = finv_[i - 1] * inv_[i]; } } static T fac(int n) { assert(n >= 0); if (n >= T::mod()) return 0; reserve(n); return fac_[n]; } static T finv(int n) { assert(n < T::mod()); if (n < 0) return 0; reserve(n); return finv_[n]; } static T inv(int n) { n %= T::mod(); if (n < 0) n += T::mod(); assert(n != 0); reserve(n); return inv_[n]; } static T P(int n, int k) { if (n < k) return 0; if (n < 0 || k < 0) return 0; if (n >= T::mod()) return 0; reserve(n); return fac_[n] * finv_[n - k]; } static T C(int n, int k) { if (n < k) return 0; if (n < 0 || k < 0) return 0; if (n >= T::mod()) return 0; reserve(n); return fac_[n] * finv_[k] * finv_[n - k]; } static T H(int n, int k) { if (n == 0 && k == 0) return 1; return C(n + k - 1, k); } }; template mint stom(string s) { mint res = 0; fec(c : s) { res *= 10; res += c - '0'; } return res; } template > pair svp2d(const pair &a, const pair &b) { assert((a != pair{0, 0} && b != pair{0, 0})); auto [a1, a2] = a; auto [b1, b2] = b; if ((U)a1 * a1 + (U)a2 * a2 < (U)b1 * b1 + (U)b2 * b2) swap(a1, b1), swap(a2, b2); while ((U)a1 * a1 + (U)a2 * a2 > (U)b1 * b1 + (U)b2 * b2) { swap(a1, b1), swap(a2, b2); T k = divround((U)a1 * b1 + (U)a2 * b2, (U)a1 * a1 + (U)a2 * a2); b1 -= k * a1, b2 -= k * a2; if (b1 == 0 && b2 == 0) return {a1, a2}; } return {a1, a2}; } template pair mint_to_rat(const mint &x) { auto [p, q] = svp2d({x.val(), 1}, {mint::mod(), 0}); if (q < 0) p = -p, q = -q; return {p, q}; } #if __cplusplus >= 202302L namespace cpp_dump { struct mint_to_rat_fn { template constexpr auto operator()(const T &x) const -> decltype(mint_to_rat(x)) { return mint_to_rat(x); } }; struct rat_closure : std::ranges::range_adaptor_closure { template constexpr auto operator()(T &&t) const { if constexpr (!std::ranges::range && std::invocable(t))>) { return mint_to_rat_fn{}(std::forward(t)); } else if constexpr (std::ranges::range) { using Ref = std::ranges::range_reference_t; if constexpr (std::invocable(std::declval()))>) { return std::forward(t) | std::views::transform(mint_to_rat_fn{}); } else if constexpr (std::ranges::range) { return std::forward(t) | std::views::transform([this](auto &&inner) { return (*this)(std::forward(inner)); }); } else { static_assert(false); } } else { static_assert(false); } } }; template requires(!std::ranges::range && std::invocable) constexpr auto operator|(T &&t, const rat_closure &c) { return c(std::forward(t)); } constexpr rat_closure rat() { return rat_closure{}; } } #endif using mint = modint998244353; using bi = Binomial; void init() { oj(mt.seed(random_device()())); } template constexpr tuple crt2(R0 r0_, R1 r1_, M0 m0_, M1 m1_) { T m0 = m0_, m1 = m1_; assert(m0 >= 1 && m1 >= 1); T r0 = safemod(r0_, m0), r1 = safemod(r1_, m1); if (m0 < m1) swap(r0, r1), swap(m0, m1); if (m0 % m1 == 0) { if (r0 % m1 != r1) return {false, 0, 0}; return {true, r0, m0}; } auto [g, im, _] = extgcd(m0, m1); T u1 = m1 / g; if ((r1 - r0) % g) return {false, 0, 0}; T x = (r1 - r0) / g % u1 * im % u1; r0 += x * m0; m0 *= u1; if (r0 < 0) r0 += m0; return {true, r0, m0}; } template constexpr tuple crt(const V1 &rs, const V2 &ms) { assert(rs.size() == ms.size()); const int n = rs.size(); T r = 0, m = 1; repi(i, n) { auto [ok, nr, nm] = crt2(r, rs[i], m, ms[i]); if (!ok) return {false, 0, 0}; r = nr, m = nm; } return {true, r, m}; } template pair crt_mod(const V1 &rs, const V2 &ms) { using T = decay_t; assert(rs.size() == ms.size()); const int n = rs.size(); mint r = 0, m = 1; if constexpr (sizeof(T) <= 4) { vc ba; ba.reserve(n); repi(i, n) ba.eb(ms[i]); vc rr(n, 0), mm(n, 1); repi(i, n) { assert(ms[i] >= 1); auto [g, im, _] = extgcd(mm[i], ms[i]); assert(g == 1); if (im < 0) im += ms[i]; ll diff = safemod(ll(rs[i]) - rr[i], ms[i]); uint t = ba[i].mul(diff, im); r += t * m, m *= ms[i]; repi(j, i + 1, n) { rr[j] += ba[j].mul(t, mm[j]); if (rr[j] >= (uint)ms[j]) rr[j] -= ms[j]; mm[j] = ba[j].mul(mm[j], ms[i]); } } } else { vc rr(n, 0), mm(n, 1); repi(i, n) { assert(ms[i] >= 1); auto [g, im, _] = extgcd(mm[i], ms[i]); assert(g == 1); if (im < 0) im += ms[i]; i128 diff = safemod((i128)rs[i] - rr[i], ms[i]); ull t = (ull)((u128)diff * im % ms[i]); r += t * m, m *= ms[i]; repi(j, i + 1, n) { rr[j] += (ull)((u128)t * mm[j] % ms[j]); if (rr[j] >= (ull)ms[j]) rr[j] -= ms[j]; mm[j] = (ull)((u128)mm[j] * ms[i] % ms[j]); } } } return {r, m}; } template constexpr pair crt_mod_constexpr(const array &rs, const array &ms) { using T = larger_int_t; assert(rs.size() == ms.size()); mint r = 0, m = 1; array rr{}, mm; fill(ALL(mm), 1); repi(i, n) { assert(ms[i] >= U2(1)); assert(U1(0) <= rs[i] && U2(rs[i]) < ms[i]); auto [g, im, _] = extgcd(mm[i], ms[i]); assert(g == 1); T t = safemod((rs[i] - rr[i]) * im, ms[i]); r += t * m, m *= ms[i]; repi(j, i + 1, n) { rr[j] += t * mm[j] % ms[j]; if (rr[j] >= ms[j]) rr[j] -= ms[j]; mm[j] *= ms[i], mm[j] %= ms[j]; } } return {r, m}; } template bool pre_crt(const V1 &rs, V2 &ms) { using T = typename V2::value_type; assert(rs.size() == ms.size()); const int n = rs.size(); repi(i, n) repi(j, i + 1, n) { T g = gcd(ms[i], ms[j]); if ((rs[i] - rs[j]) % g) return false; ms[i] /= g, ms[j] /= g; T gi = gcd(ms[i], g), gj = g / gi; do { g = gcd(gi, gj); gi *= g, gj /= g; } while (g > 1); ms[i] *= gi, ms[j] *= gj; } return true; } template T convolution_point_get(const vc &a, const vc &b, int p) { const int n = a.size(), m = b.size(); T res = 0; repi(i, max(0, p - m + 1), min(n, p + 1)) res += a[i] * b[p - i]; return res; } namespace internal { constexpr int primitive_root_constexpr(int m) { if (m == 2) return 1; if (m == 167772161) return 3; if (m == 469762049) return 3; if (m == 754974721) return 11; if (m == 998244353) return 3; if (m == 1107296257) return 10; if (m == 1711276033) return 29; if (m == 1811939329) return 13; if (m == 2013265921) return 31; if (m == 2113929217) return 5; int divs[20] = {}; divs[0] = 2; int cnt = 1; int x = (m - 1) / 2; while (x % 2 == 0) x /= 2; for (int i = 3; (long long)(i)*i <= x; i += 2) { if (x % i == 0) { divs[cnt++] = i; while (x % i == 0) { x /= i; } } } if (x > 1) { divs[cnt++] = x; } for (int g = 2;; g++) { bool ok = true; for (int i = 0; i < cnt; i++) { if (powmod_constexpr(g, (m - 1) / divs[i], m) == 1) { ok = false; break; } } if (ok) return g; } } template constexpr int primitive_root_for_convolution = primitive_root_constexpr(m); template > struct fft_info { static constexpr int rank2 = countr_zero(mint::mod() - 1); std::array root; // root[i]^(2^i) == 1 (i <= rank2) std::array iroot; // root[i] * iroot[i] == 1 (i <= rank2) std::array rate2; std::array irate2; std::array rate3; std::array irate3; fft_info() { root[rank2] = mint(g).pow((mint::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]; } { mint 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]; } } { mint 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]; } } } }; } // namespace internal template bool ntt_ok(int n) { if (n <= 0) return false; if constexpr (is_static_modint_v) { if constexpr (!internal::isprime) return false; static constexpr int rank2 = countr_zero(mint::mod() - 1); return n <= (1 << rank2); } else return false; } template void ntt(vc> &) { assert(false); } template void intt(vc> &) { assert(false); } template void ntt(vc>> &a) { using mint = internal::modint_impl>; int n = int(a.size()); assert(n > 0); int h = countr_zero((unsigned int)n); assert(n == (1 << h)); assert(ntt_ok(n)); static const internal::fft_info info; 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); mint 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 *= info.rate2[countr_zero(~(unsigned int)(s))]; } len++; } else { int p = 1 << (h - len - 2); mint rot = 1, imag = info.root[2]; for (int s = 0; s < (1 << len); s++) { mint rot2 = rot * rot; mint rot3 = rot2 * rot; int offset = s << (h - len); for (int i = 0; i < p; i++) { auto mod2 = 1ULL * mint::mod() * mint::mod(); auto a0 = 1ULL * a[i + offset].val(); auto a1 = 1ULL * a[i + offset + p].val() * rot.val(); auto a2 = 1ULL * a[i + offset + 2 * p].val() * rot2.val(); auto a3 = 1ULL * a[i + offset + 3 * p].val() * rot3.val(); auto a1na3imag = 1ULL * mint(a1 + mod2 - a3).val() * imag.val(); 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 *= info.rate3[countr_zero(~(unsigned int)(s))]; } len += 2; } } } template void intt(vc>> &a) { using mint = internal::modint_impl>; int n = int(a.size()); assert(n > 0); int h = countr_zero((unsigned int)n); assert(n == (1 << h)); assert(ntt_ok(n)); static const internal::fft_info info; int len = h; // a[i, i+(n>>len), i+2*(n>>len), ..] is transformed while (len) { if (len == 1) { int p = 1 << (h - len); mint 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)mint::mod() + l.val() - (uint)r.val()) * irot.val(); ; } if (s + 1 != (1 << (len - 1))) irot *= info.irate2[countr_zero(~(unsigned int)(s))]; } len--; } else { int p = 1 << (h - len); mint irot = 1, iimag = info.iroot[2]; for (int s = 0; s < (1 << (len - 2)); s++) { mint irot2 = irot * irot; mint irot3 = irot2 * irot; int offset = s << (h - len + 2); for (int i = 0; i < p; i++) { auto a0 = 1ULL * a[i + offset + 0 * p].val(); auto a1 = 1ULL * a[i + offset + 1 * p].val(); auto a2 = 1ULL * a[i + offset + 2 * p].val(); auto a3 = 1ULL * a[i + offset + 3 * p].val(); auto a2na3iimag = 1ULL * mint((mint::mod() + a2 - a3) * iimag.val()).val(); a[i + offset] = a0 + a1 + a2 + a3; a[i + offset + 1 * p] = (a0 + (mint::mod() - a1) + a2na3iimag) * irot.val(); a[i + offset + 2 * p] = (a0 + a1 + (mint::mod() - a2) + (mint::mod() - a3)) * irot2.val(); a[i + offset + 3 * p] = (a0 + (mint::mod() - a1) + (mint::mod() - a2na3iimag)) * irot3.val(); } if (s + 1 != (1 << (len - 2))) irot *= info.irate3[countr_zero(~(unsigned int)(s))]; } len -= 2; } } } namespace internal { template vc convolution_naive(const vc &a, const vc &b) { const int n = a.size(), m = b.size(); const int cnta = n - count(ALL(a), 0), cntb = m - count(ALL(b), 0); vc c(n + m - 1); if ((ll)m * cnta > (ll)n * cntb) { repi(j, m) { if (b[j] == 0) continue; repi(i, n) c[i + j] += a[i] * b[j]; } } else { repi(i, n) { if (a[i] == 0) continue; repi(j, m) c[i + j] += a[i] * b[j]; } } return c; } template vc convolution_ntt(vc a, vc b) { const int n = a.size(), m = b.size(); const int z = bit_ceil(n + m - 1); if (a == b) { a.resize(z); ntt(a); repi(i, z) a[i] *= a[i]; } else { a.resize(z), b.resize(z); ntt(a), ntt(b); repi(i, z) a[i] *= b[i]; } intt(a); mint iz = mint(z).inv(); fem(ai : a) ai *= iz; a.resize(n + m - 1); return a; } template void convolution_crt_helper(const vc &a, const vc &b, vc> &cs) { using mint = static_modint32; const int n = a.size(), m = b.size(); auto c = convolution_ntt(vc(ALL(a)), vc(ALL(b))); repi(i, n + m - 1) cs[i][j] = c[i].val(); } template vc convolution_crt(const vc &a, const vc &b) { const int n = a.size(), m = b.size(); constexpr size_t k = sizeof...(ms); vc> cs(n + m - 1); constexpr array ms_arr = {ms...}; [&](index_sequence) { (convolution_crt_helper(a, b, cs), ...); }(make_index_sequence{}); vc c(n + m - 1); repi(i, n + m - 1) c[i] = get<1>(crt(cs[i], ms_arr)); return c; } template vc convolution_crt_mod(const vc &a, const vc &b) { const int n = a.size(), m = b.size(); constexpr size_t k = sizeof...(ms); vc> cs(n + m - 1); constexpr array ms_arr = {ms...}; [&](index_sequence) { (convolution_crt_helper(a, b, cs), ...); }(make_index_sequence{}); vc c(n + m - 1); repi(i, n + m - 1) c[i] = crt_mod_constexpr(cs[i], ms_arr).first; return c; } } // namespace internal template ::value>> vc convolution(const vc &a, const vc &b) { const int n = a.size(), m = b.size(); const int cnta = n - count(ALL(a), 0), cntb = m - count(ALL(b), 0); if (n == 0 || m == 0) return {}; if (ntt_ok(n + m - 1)) { if (min(cnta, cntb) <= 60) return internal::convolution_naive(a, b); return internal::convolution_ntt(a, b); } else { if (min(cnta, cntb) <= 300) return internal::convolution_naive(a, b); assert(ntt_ok>(n + m - 1) && "|a| + |b| - 1 <= 2^26"); vc a_(n), b_(m); repi(i, n) a_[i] = a[i].val(); repi(j, m) b_[j] = b[j].val(); return internal::convolution_crt_mod(a_, b_); } } template ::value>> vc convolution(const vc &a, const vc &b) { using mint = static_modint32; auto c = convolution(vc(ALL(a)), vc(ALL(b))); vc c_(c.size()); repi(i, c.size()) c_[i] = c[i].val(); return c_; } void main2() { LL(N, M, K); bi::reserve(1000); vc f(N + 1); rep(i, 1, N + 1) f[i] = bi::finv_[i]; vc g(N + 1); auto dp = dvec({M + 2, M + 2, N + 1}, mint(0)); dp[M + 1][M + 1][0] = bi::fac(N); rep(k, M + 1, K + 1) dp[M + 1][M + 1][0] *= k; rep(k, M, 0, -1) { rep(x, M + 2) { dp[k][x] = convolution(f, dp[k + 1][x]); dp[k][x].resize(N + 1); } rep(n, N + 1) rep(x, M + 2) { if (dp[k + 1][x][n] == 0) continue; { ll nn = n; ll ck = nn - n; ll nx = ck == 0 ? k : x; mint coef = bi::finv(ck); dp[k][nx][nn] += coef * dp[k + 1][x][n]; } } rep(nx, M + 2) rep(nn, N + 1) dp[k][nx][nn] *= (k <= K ? nx : 1); } mint ans = 0; rep(x, M + 2) ans += dp[1][x][N]; PRINT(ans); } void test() { } template struct Main { Main() { cauto CERR = [](string val, string color) { string s = "\033[" + color + "m" + val + "\033[m"; /* コードテストで確認する際にコメントアウトを外す cerr << val; //*/ }; CERR("\n[FAST_CIO]\n\n", "32"); cin.tie(0); ios::sync_with_stdio(false); cout << fixed << setprecision(20); init(); CERR("\n[SINGLE_TESTCASE]\n\n", "36"); main2(); } }; Main main_dummy; } int main() {}