#define SINGLE_TESTCASE // #define MULTI_TESTCASE // #define AOJ_TESTCASE #define FAST_IO // #define FAST_CIO // #define INTERACTIVE #define INF 4'000'000'000'000'000'037LL #define EPS 1e-11 // https://github.com/miscalculation53/library/tree/wip/template/template_all.hpp // https://github.com/miscalculation53/library/tree/wip/template/template_all_but_modint.hpp // https://github.com/miscalculation53/library/tree/wip/template/template_types.hpp #include namespace { using namespace std; using ll = long long; using uint = unsigned int; using ull = unsigned long long; using pll = pair; using tlll = tuple; #define vc vector template using vvc = vc>; using vstr = vc; 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; } #define cauto const auto // https://github.com/miscalculation53/library/tree/wip/template/template_rep.hpp #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__) // https://github.com/miscalculation53/library/tree/wip/template/template_math.hpp // https://github.com/miscalculation53/library/tree/wip/utils/is_integral_ext.hpp 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; // https://github.com/miscalculation53/library/tree/wip/utils/integer_arithmetic.hpp // https://github.com/miscalculation53/library/tree/wip/utils/make_unsigned_ext.hpp template constexpr bool is_negative(T x) { static_assert(is_integral_ext); if constexpr (is_signed_ext) return x < 0; return false; } // https://github.com/miscalculation53/library/tree/wip/utils/default_infty.hpp namespace default_infty_detail { template struct has_infty : false_type {}; template struct has_infty> : true_type {}; template const T &custom_value() { static const T value = T::infty(); return value; } template inline constexpr bool unsupported = false; } template constexpr decltype(auto) default_infty() { if constexpr (default_infty_detail::has_infty::value) return default_infty_detail::custom_value(); else if constexpr (is_same_v || is_same_v) return T(INF) * T(INF); else if constexpr (is_integral_ext && !is_same_v) { if constexpr (sizeof(T) >= sizeof(ll)) return T(INF); else if constexpr (numeric_limits::digits >= 31) return (T(1) << 30) - 1; else return T(numeric_limits::max() / 2); } else if constexpr (numeric_limits::has_infinity) return numeric_limits::infinity(); else static_assert(default_infty_detail::unsupported, "No default infinity for this type; specify infty explicitly or define T::infty()."); } namespace template_math_detail { } template && is_integral_ext>> inline constexpr T safemod(U a, V b) { const T x = T(a), y = T(b); if constexpr (is_signed_ext) if (y == T(-1)) return T(0); const T r = x % y; return r != T(0) && is_negative(x) != is_negative(y) ? T(r + y) : r; } // https://github.com/miscalculation53/library/tree/wip/template/template_vector.hpp #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) ([&](const auto &x) { return fx; }) #define GEN_VEC(n,i,fi) (gen_vec(n, LMD(i, fi))) template void fill(V &v, const T &value) { if constexpr (is_assignable_v) std::fill(std::begin(v), std::end(v), value); else for (auto &row : v) fill(row, value); } // https://github.com/miscalculation53/library/tree/wip/template/template_algo.hpp // https://github.com/miscalculation53/library/tree/wip/utils/resolved_infty.hpp // https://github.com/miscalculation53/library/tree/wip/utils/resolved_value.hpp template , int> = 0> constexpr T resolved_value() { return T(x); } template , long> = 0> const T &resolved_value() { static const T value = T(x()); return value; } template constexpr decltype(auto) resolved_infty() { if constexpr (is_same_v) return default_infty(); else return resolved_value(); } 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 {}; } }; namespace internal { } 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 cumlsum(const vc &v, int left_index = 0) { return cuml>(v, left_index); } constexpr array DRULgrid = {{{1, 0}, {0, 1}, {-1, 0}, {0, -1}}}; constexpr array DRULplane = {{{0, -1}, {1, 0}, {0, 1}, {-1, 0}}}; // https://github.com/miscalculation53/library/tree/wip/template/template_binsearch.hpp 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; #define DEFAULT_COMP ranges::less namespace internal { }; // https://github.com/miscalculation53/library/tree/wip/template/template_bit.hpp template inline constexpr ull MASK(T k) { return (1ULL << k) - 1ULL; } 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); } 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; } 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; } // https://github.com/miscalculation53/library/tree/wip/template/template_inout.hpp // https://github.com/miscalculation53/library/tree/wip/template/template_dump.hpp // https://github.com/miscalculation53/library/tree/wip/template/template_dump_map.hpp #define CPP_DUMP_DEFINE_DATA(...) #define dump(...) ((void)0) #define local(...) #define oj(...) __VA_ARGS__ #define local_oj(a,b) (b) namespace fastio { template struct unsigned_integer { using type = make_unsigned_t; }; template <> struct unsigned_integer { using type = u128; }; template <> struct unsigned_integer { using type = u128; }; template using unsigned_integer_t = typename unsigned_integer::type; static constexpr uint32_t SIZ = 1 << 17; char ibuf[SIZ]; char obuf[SIZ]; char out[100]; 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() { memcpy(ibuf, ibuf + pil, pir - pil); pir = pir - pil + fread(ibuf + pir - pil, 1, SIZ - pir + pil, stdin); pil = 0; if (pir < SIZ) ibuf[pir++] = '\n'; } inline void flush() { fwrite(obuf, 1, por, stdout); por = 0; } void rd1(char &c) { do { if (pil + 1 > pir) load(); c = ibuf[pil++]; } while (c <= ' '); } void rd1(string &x) { x.clear(); while (true) { if (pil == pir) load(); while (pil < pir && ibuf[pil] <= ' ') ++pil; if (pil < pir) break; } while (true) { uint32_t p = pil; while (pil < pir && ibuf[pil] > ' ') ++pil; x.append(ibuf + p, pil - p); if (pil < pir) { ++pil; return; } load(); } } template void rd1_real(T &x) { string s; rd1(s); if constexpr (!is_same_v) { auto [p, ec] = from_chars(s.data(), s.data() + s.size(), x); if (ec == errc{} && p == s.data() + s.size()) return; } if constexpr (is_same_v) x = stold(s); else x = stod(s); } template void rd1_integer(T &x) { using U = unsigned_integer_t; bool minus = false; U val = 0; if constexpr (check_buffer) if (pil + 100 > pir) load(); uint32_t p = pil; while (ibuf[p] < '-') ++p; if constexpr (is_signed::value || is_same_v) { if (ibuf[p] == '-') minus = true, ++p; } while ('0' <= ibuf[p]) val = val * 10 + (ibuf[p++] & 15); pil = p; if constexpr (is_signed::value || is_same_v) { if (minus) { const U min_abs = U(numeric_limits::max()) + 1; x = val == min_abs ? numeric_limits::lowest() : -T(val); } else x = T(val); } else x = T(val); } void rd1(int &x) { rd1_integer(x); } void rd1(ll &x) { rd1_integer(x); } void rd1(i128 &x) { rd1_integer(x); } void rd1(uint &x) { rd1_integer(x); } void rd1(ull &x) { rd1_integer(x); } void rd1(u128 &x) { rd1_integer(x); } void rd1(double &x) { rd1_real(x); } void rd1(long double &x) { rd1_real(x); } template void rd1(tuple &tpl) { apply([](auto &...x) { (rd1(x), ...); }, tpl); } template void rd1(vc &x) { for (auto &d: x) rd1(d); } template void read(T &...x) { if constexpr (sizeof...(T) <= SIZ / 100 && ((!is_same_v && (is_integral_v || is_same_v || is_same_v)) && ...)) { if (pil + 100 * sizeof...(T) > pir) load(); (rd1_integer(x), ...); } else (rd1(x), ...); } void wt1(const char c) { if (por == SIZ) flush(); obuf[por++] = c; } void wt1(string_view s) { while (!s.empty()) { if (por == SIZ) flush(); size_t n = min(s.size(), SIZ - por); memcpy(obuf + por, s.data(), n); por += n; s.remove_prefix(n); } } template void wt1_integer(T x) { if (por > SIZ - 100) flush(); using U = unsigned_integer_t; U ux; if constexpr (is_signed::value || is_same_v) { if (x < 0) obuf[por++] = '-', ux = U(0) - U(x); else ux = U(x); } else ux = x; int outi; for (outi = 96; ux >= 10000; outi -= 4) { memcpy(out + outi, pre.num[ux % 10000], 4); ux /= 10000; } if (ux >= 1000) { memcpy(obuf + por, pre.num[ux], 4); por += 4; } else if (ux >= 100) { memcpy(obuf + por, pre.num[ux] + 1, 3); por += 3; } else if (ux >= 10) { int q = (ux * 103) >> 10; obuf[por] = q | '0'; obuf[por + 1] = (ux - q * 10) | '0'; por += 2; } else obuf[por++] = ux | '0'; memcpy(obuf + por, out + outi + 4, 96 - outi); por += 96 - outi; } template void wt1_real(T x) { if constexpr (!is_same_v) { auto [p, ec] = to_chars(out, out + sizeof(out), x, chars_format::fixed, 15); if (ec == errc{}) { wt1(string_view(out, p)); return; } } ostringstream oss; oss << fixed << setprecision(15) << x; wt1(oss.str()); } void wt1(int x) { wt1_integer(x); } template , int> = 0> void wt1(T x) { wt1_integer(x); } void wt1(i128 x) { wt1_integer(x); } void wt1(u128 x) { wt1_integer(x); } void wt1(double x) { wt1_real(x); } void wt1(long double x) { wt1_real(x); } template void wt1(const pair &val); template void wt1(const tuple &tpl); template void wt1(const array &val); template void wt1(const vector &val); template void write(T &&...x) { (wt1(std::forward(x)), ...); } template void print(T &&...x) { if constexpr (sizeof...(T)) { int i = 0; ((i++ ? wt1(' ') : void(), wt1(std::forward(x))), ...); } wt1('\n'); } } struct Dummy { Dummy() { atexit(fastio::flush); } } dummy; namespace internal { template void READnodump(Ts &...a) { fastio::read(a...); } template void READVECnodump(int n, vc &...v) { (v.resize(n), ...); READnodump(v...); } }; #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' #define WRITE fastio::write #define PRINT fastio::print #define PRINTEXIT(...) do { PRINT(__VA_ARGS__); exit(0); } while (false) #define PRINTRETURN(...) do { PRINT(__VA_ARGS__); return; } while (false) #define PRINTVEXIT(...) do { PRINTV(__VA_ARGS__); exit(0); } while (false) #define PRINTVRETURN(...) do { PRINTV(__VA_ARGS__); return; } while (false) namespace internal { template auto &tuple_add_impl(A &a, const B &b, const index_sequence) { ((get(a) += get(b)), ...); return a; } }; template tuple &operator+=(tuple &a, const Tp &b) { return internal::tuple_add_impl(a, b, make_index_sequence>>{}); } template void offset(vc &v, const Add &add) { for (auto &vi : v) vi += 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 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; } 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 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; } namespace internal { }; #define UNZIP(vt,...) auto [__VA_ARGS__] = unzip(vt) #define ZIP(vt,...) auto vt = internal::zip_vectors(__VA_ARGS__) // https://github.com/miscalculation53/library/tree/wip/template/template_random.hpp mt19937_64 mt; bool randbool(double p) { assert(0 <= p && p <= 1); return bernoulli_distribution(p)(mt); } namespace internal { }; // https://github.com/miscalculation53/library/tree/wip/math/modint/template_modint.hpp // https://github.com/miscalculation53/library/tree/wip/math/modint/modint.hpp // https://github.com/miscalculation53/library/tree/wip/math/modint/modint_internal_static.hpp // https://github.com/miscalculation53/library/tree/wip/utils/larger_int.hpp namespace larger_int_detail { } template struct larger_int { private: static constexpr bool check(); static_assert(check()); public: 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 struct Rational; template struct larger_int> { using type = Rational::type>; }; template using larger_int_t = typename larger_int::type; // https://github.com/miscalculation53/library/tree/wip/math/modint/modint_internal_isprime.hpp 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 = isprime_constexpr(n); }; namespace internal { template struct policy_static { using mod_type = decltype(M); using value_type = make_unsigned_t; using calc_type = conditional_t<(M <= UINT_MAX), ull, larger_int_t>; static constexpr bool is_prime = isprime_constexpr(M); static constexpr value_type init(value_type v) { return v; } static constexpr mod_type val(value_type v); }; }; // https://github.com/miscalculation53/library/tree/wip/math/modint/modint_internal_barrett32.hpp namespace internal { struct barrett32 { uint m = 0; ull im = 0; barrett32() = default; 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{}; static value_type init(value_type v) { return v; } static mod_type val(value_type v); }; }; // https://github.com/miscalculation53/library/tree/wip/math/modint/modint_internal_montgomery64.hpp 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 = 0, im = 0, sq = 0; montgomery64odd() = default; 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 = 0, mx = 0, imx = 0, d = 0, q = 0; uint b = 0; montgomery64() = default; 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{}; static value_type init(value_type v) { return reducer.inv_reduce(v); } static mod_type val(value_type v); }; 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{}; static value_type init(value_type v) { return reducer.inv_reduce(v); } static mod_type val(value_type v); }; }; // https://github.com/miscalculation53/library/tree/wip/math/extgcd.hpp 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: modint_impl(); template >> modint_impl(T v); M val() const; 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 M safe_hash_key(const mint &x) { return x.val(); } friend void rd1(mint &x) { long long a; fastio::rd1(a); x = a; } friend void wt1(const mint &x) { fastio::wt1(x.val()); } }; }; template using static_modint32 = internal::modint_impl>; template using dynamic_modint32 = 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, PowInv>> : true_type {}; template struct is_static_modint, PowInv>> : 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; // https://github.com/miscalculation53/library/tree/wip/math/modint/power_table.hpp // https://github.com/miscalculation53/library/tree/wip/math/modint/binomial.hpp template struct Binomial { private: inline static decltype(T::mod()) mod; public: inline static vc fac_, finv_, inv_; }; // https://github.com/miscalculation53/library/tree/wip/math/modint/stom.hpp // https://github.com/miscalculation53/library/tree/wip/math/modint/to_rational.hpp // https://github.com/miscalculation53/library/tree/wip/math/svp2d.hpp namespace cpp_dump { template struct rat_value { T p, q; friend ostream &operator<<(ostream &os, const rat_value &x) { os << x.p; if (x.q != 1) os << '/' << x.q; return os; } }; template inline constexpr bool rat_is_map = false; template inline constexpr bool rat_is_map> = true; template inline constexpr bool rat_is_map> = true; template inline constexpr bool rat_is_map> = true; template inline constexpr bool rat_is_map> = true; struct rat_closure : std::ranges::range_adaptor_closure { }; constexpr rat_closure rat() { return rat_closure{}; } } using mint = modint998244353; // using mint = modint1000000007; // using mint = static_modint<1000000000>; // using mint = modint; void init() { oj(mt.seed(random_device()())); } // https://github.com/miscalculation53/library/tree/wip/graph/maxflow.hpp // https://github.com/miscalculation53/library/tree/wip/ds/csr.hpp template struct CSR { protected: int n, m; vc start; vc elist; vc len; inline int get_last(int i) const { if constexpr (is_erasable) return start[i] + len[i]; else return start[i + 1]; } template struct RowBase { using iterator = Iter; using reference = typename iterator_traits::reference; private: iterator begi, endi; public: RowBase(const iterator &begi, const iterator &endi) : begi(begi), endi(endi) {} inline iterator begin() const { return begi; } inline iterator end() const { return endi; } template inline I size() const { return endi - begi; } inline bool empty() const { return size() == 0; } inline reference operator[](int i) const { return *(begi + i); } inline reference back() const { assert(!empty()); return *prev(endi); } }; using Row = RowBase::iterator>; using ConstRow = RowBase::const_iterator>; public: CSR() : n(0), m(0), start{0} {} CSR(const vc &row_sizes) : n(row_sizes.size()) { fec(s : row_sizes) assert(s >= 0); start = cumlsum(row_sizes); m = start.back(); elist.resize(m); if constexpr (is_erasable) len = row_sizes; } template CSR(int n, const vc> &ies) : n(n), m(ies.size()), start(n, 0), elist(m) { if constexpr (is_erasable) len.resize(n); fec([ i, e ] : ies) { assert(0 <= i && i < n); start[i]++; } start = cumlsum(start); if constexpr (is_erasable) repi(i, n) len[i] = start[i + 1] - start[i]; auto cnt = start; repi(j, m) { cauto & [ i, e ] = ies[j]; int &k = cnt[i]; elist[k] = e; k++; } } CSR(const vvc &vv) : n(vv.size()), start(n + 1, 0) { m = 0; fec(row : vv) m += row.size(); elist.resize(m); if constexpr (is_erasable) len.resize(n); for (int i = 0, j = 0; i < n; i++) { start[i] = j; if constexpr (is_erasable) len[i] = vv[i].size(); fec(e : vv[i]) { elist[j] = e; j++; } } start.back() = m; } Row operator[](int i) { return Row(elist.begin() + start[i], elist.begin() + get_last(i)); } int offset(int i) const { assert(0 <= i && i <= n); return start[i]; } template I size() const { return n; } vc &get_elist() { return elist; } }; template struct MaxFlow { static_assert(is_integral_ext && !is_same_v); struct edge { int from, to; Cap cap, flow; CPP_DUMP_DEFINE_DATA(from, to, cap, flow); }; private: struct ResidualEdge { int to, rev; Cap cap; }; int n = 0; mutable CSR g; mutable vc pos; mutable vc pending; mutable CSR originals; mutable vc merged_pos; mutable vc merged_flow; mutable int built = 0; mutable bool dirty = false; bool merge_parallel_edges() const { array seen{}; const int m = num_of_edges(), samples = min(m, 1024); int repeated = 0, omitted = 0; for (int j = 0; j < samples; j++) { const int i = m <= 1024 ? j : int(ll(j / 16) * m / 64) + j % 16; const auto e = get_edge(i); if (e.from == e.to || e.cap == 0) { omitted++; continue; } const uint64_t key = (uint64_t(e.from) << 32) | uint32_t(e.to); const int slot = (key * 0x9e3779b97f4a7c15ULL) >> 54; if (seen[slot] == key) repeated++; seen[slot] = key; } return 2 * omitted <= samples && repeated >= max(1, samples / 64); } Cap finish_flow(Cap result) { const auto &data = g.get_elist(); for (int k = 0; k < int(merged_pos.size()); k++) { const Cap current = data[data[merged_pos[k]].rev].cap; const Cap previous = merged_flow[k]; if (current == previous) continue; Cap remaining = current > previous ? current - previous : previous - current; for (auto &e : originals[k]) { const Cap f = min(remaining, current > previous ? e.cap - e.flow : e.flow); if (current > previous) e.flow += f; else e.flow -= f; remaining -= f; if (remaining == 0) break; } assert(remaining == 0); merged_flow[k] = current; } return result; } void check_flow(int s, int t, Cap limit) const { assert(0 <= s && s < n && 0 <= t && t < n && s != t); assert(Cap(0) <= limit); } Cap flow_bound(int s, int t, Cap limit) const { auto cut_bound = [&](int v, bool reverse) { Cap sum = 0; for (const auto &e : g[v]) if (e.to != v) { const Cap c = reverse ? g.get_elist()[e.rev].cap : e.cap; if (c >= limit - sum) return limit; sum += c; } return sum; }; limit = min(limit, cut_bound(s, false)); return min(limit, cut_bound(t, true)); } int dinic_budget_factor(int s, int t, Cap progress) const { Cap unit = 0; auto check = [&](const vc &use) { auto selected = [&](int v) { return use.empty() || use[v]; }; unit = [&] { Cap result = 0; for (int v = 0; v < n; v++) if (selected(v)) for (const auto &e : g[v]) if (e.to != v && selected(e.to)) { Cap a = e.cap; while (a != 0) { const Cap b = result % a; result = a, a = b; } if (result == 1) return result; } return result; }(); if (unit == 0) return true; for (int v = 0; v < n; v++) if (v != s && v != t && selected(v)) { int in = 0, out = 0; auto add = [&](int &sum, Cap cap) { if (sum > 8) return; const Cap value = unit == 1 ? cap : cap / unit; sum = value > Cap(8 - sum) ? 9 : sum + int(value); }; for (const auto &e : g[v]) if (e.to != v && selected(e.to)) { add(out, e.cap), add(in, g.get_elist()[e.rev].cap); if (in > 8 && out > 8) return false; } } return true; }; if (check({})) return 8; auto use = min_cut(s); vc reaches_t(n, false); vc que(n); int front = 0, back = 0; reaches_t[t] = true, que[back++] = t; while (front < back) { const int v = que[front++]; for (const auto &e : g[v]) if (!reaches_t[e.to] && g.get_elist()[e.rev].cap > 0) reaches_t[e.to] = true, que[back++] = e.to; } for (int v = 0; v < n; v++) use[v] = use[v] && reaches_t[v]; auto small_remaining_cut = [&](int v, bool reverse) { Cap remaining = progress; for (const auto &e : g[v]) if (e.to != v && use[e.to]) { const Cap cap = reverse ? g.get_elist()[e.rev].cap : e.cap; if (cap > remaining) return false; remaining -= cap; } return true; }; if (small_remaining_cut(s, false) || small_remaining_cut(t, true)) return 8; if (check(use)) return 8; vc large_reachable(n, false); front = back = 0; large_reachable[s] = true, que[back++] = s; while (front < back) { const int v = que[front++]; for (const auto &e : g[v]) if (use[e.to] && !large_reachable[e.to]) { const Cap value = unit == 1 ? e.cap : e.cap / unit; if (value <= 8) continue; if (e.to == t) return 1; large_reachable[e.to] = true, que[back++] = e.to; } } int cut = 0; for (int v = 0; v < n; v++) if (large_reachable[v]) for (const auto &e : g[v]) if (use[e.to] && !large_reachable[e.to]) { const Cap value = unit == 1 ? e.cap : e.cap / unit; if (value > Cap(32 - cut)) return 1; cut += int(value); } return 64; } template Cap run_dinic(int s, int t, Cap limit, ll &budget, bool &finished, int *structure = nullptr) { vc level(n), iter(n), que(n), path; vc bound, sent; Cap total = 0; ll work = 0; finished = true; while (total < limit) { fill(ALL(level), -1); if constexpr (adaptive) work += n; int front = 0, back = 0; level[s] = 0; que[back++] = s; while (front < back && level[t] == -1) { const int v = que[front++]; for (const auto &e : g[v]) { if constexpr (adaptive) work++; if (e.cap == 0 || level[e.to] != -1) continue; level[e.to] = level[v] + 1; if (e.to == t) break; que[back++] = e.to; } } if (level[t] == -1) return total; ll phase_budget = budget * (level[t] <= 16 ? 8 : 2); if constexpr (adaptive) if (work >= phase_budget) { if (structure && *structure == 0) { *structure = total >= limit - total ? 8 : dinic_budget_factor(s, t, total); budget *= *structure, phase_budget *= *structure; } if (work >= phase_budget) { finished = false; return total; } } fill(ALL(iter), 0); if constexpr (adaptive) work += n; bool stop = false; if (level[t] <= 64) { auto dfs = [&](auto &self, int v, Cap up) -> Cap { if (v == s) return up; Cap res = 0; const int level_v = level[v]; const auto row = g[v]; for (int &i = iter[v]; i < int(row.size()); i++) { if constexpr (adaptive) if (++work >= phase_budget) { stop = true; return res; } auto &e = row[i]; if (level_v <= level[e.to] || g.get_elist()[e.rev].cap == 0) continue; const Cap f = self(self, e.to, min(up - res, g.get_elist()[e.rev].cap)); e.cap += f, g.get_elist()[e.rev].cap -= f; res += f; if (res == up || stop) return res; } level[v] = n; return res; }; total += dfs(dfs, t, limit - total); if (stop) { finished = false; return total; } continue; } if (path.empty()) path.resize(n), bound.resize(n), sent.resize(n); int depth = 0; path[0] = t, bound[0] = limit - total, sent[0] = 0; while (true) { const int v = path[depth]; if (v != s && sent[depth] < bound[depth] && !stop) { int &i = iter[v]; const auto row = g[v]; while (i < int(row.size())) { if constexpr (adaptive) if (++work >= phase_budget) { stop = true; break; } const auto &e = row[i]; if (level[e.to] < level[v] && g.get_elist()[e.rev].cap > 0) break; i++; } if (!stop && i < int(row.size())) { const auto &e = row[i]; const Cap next_bound = min(bound[depth] - sent[depth], g.get_elist()[e.rev].cap); path[++depth] = e.to; bound[depth] = next_bound, sent[depth] = 0; continue; } if (!stop) level[v] = n; } const Cap f = v == s ? bound[depth] : sent[depth]; if (depth == 0) { total += f; break; } const int p = path[--depth]; const auto row = g[p]; auto &e = row[iter[p]]; e.cap += f; g.get_elist()[e.rev].cap -= f; sent[depth] += f; if (sent[depth] < bound[depth]) iter[p]++; } if (stop) { finished = false; return total; } } return total; } public: MaxFlow() = default; explicit MaxFlow(int n) : n(n) { assert(0 <= n); } template I size() const { return n; } template I num_of_edges() const { return built + pending.size(); } int add_edge(int from, int to, Cap cap) { assert(0 <= from && from < n && 0 <= to && to < n); assert(Cap(0) <= cap); const int id = num_of_edges(); pending.push_back({from, to, cap, 0}); return id; } void build() const { if (pending.empty() && !dirty && g.template size() == n) return; const int m = num_of_edges(); if (!merge_parallel_edges()) { vc degree(n, 0); int omitted = 0; for (int i = 0; i < m; i++) { const auto e = get_edge(i); if (e.from == e.to || e.cap == 0) omitted++; else degree[e.from]++, degree[e.to]++; } CSR next_graph(degree); CSR next_originals(vc{omitted}); vc next_pos(m), count(n, 0); auto &data = next_graph.get_elist(); int special = 0; for (int i = 0; i < m; i++) { const auto e = get_edge(i); if (e.from == e.to || e.cap == 0) { next_originals.get_elist()[special] = e; next_pos[i] = ~special++; continue; } const int a = next_graph.offset(e.from) + count[e.from]++; const int b = next_graph.offset(e.to) + count[e.to]++; data[a] = {e.to, b, Cap(e.cap - e.flow)}; data[b] = {e.from, a, e.flow}; next_pos[i] = a; } swap(g, next_graph), swap(originals, next_originals); pos.swap(next_pos); vc().swap(merged_pos), vc().swap(merged_flow); built = pos.size(), dirty = false; vc().swap(pending); return; } vc degree(n, 0); for (int i = 0; i < m; i++) { const auto e = get_edge(i); if (e.from != e.to && e.cap > 0) degree[e.from]++; } CSR outgoing(degree); vc count(n, 0); for (int i = 0; i < m; i++) { const auto e = get_edge(i); if (e.from != e.to && e.cap > 0) outgoing.get_elist()[outgoing.offset(e.from) + count[e.from]++] = i; } vc combined; vc group(m, -1), last(n, -1), multiplicity; combined.reserve(outgoing.get_elist().size()); multiplicity.reserve(outgoing.get_elist().size()); for (int v = 0; v < n; v++) for (int i : outgoing[v]) { const auto e = get_edge(i); int k = last[e.to]; if (k == -1 || combined[k].from != v || e.cap > numeric_limits::max() - combined[k].cap) { k = combined.size(); last[e.to] = k; combined.eb(e), multiplicity.eb(1); } else { combined[k].cap += e.cap, combined[k].flow += e.flow; multiplicity[k]++; } group[i] = k; } fill(ALL(degree), 0); vc special_row(combined.size(), -1), row_sizes; vc next_merged_pos; vc next_merged_flow; for (int k = 0; k < int(combined.size()); k++) { const auto &e = combined[k]; degree[e.from]++, degree[e.to]++; if (multiplicity[k] > 1) { special_row[k] = row_sizes.size(); row_sizes.eb(multiplicity[k]); next_merged_flow.eb(e.flow); } } row_sizes.eb(m - outgoing.get_elist().size()); CSR next_originals(row_sizes); CSR next_graph(degree); vc next_pos(m), combined_pos(combined.size(), -1); next_merged_pos.resize(next_merged_flow.size()); fill(ALL(count), 0); auto &data = next_graph.get_elist(); for (int i = 0; i < m; i++) { const int k = group[i]; if (k == -1 || combined_pos[k] != -1) continue; const auto &e = combined[k]; const int a = next_graph.offset(e.from) + count[e.from]++; const int b = next_graph.offset(e.to) + count[e.to]++; data[a] = {e.to, b, Cap(e.cap - e.flow)}; data[b] = {e.from, a, e.flow}; combined_pos[k] = a; if (special_row[k] != -1) next_merged_pos[special_row[k]] = a; } count.assign(row_sizes.size(), 0); for (int i = 0; i < m; i++) { const int k = group[i]; if (k != -1 && special_row[k] == -1) next_pos[i] = combined_pos[k]; else { const int row = k == -1 ? int(row_sizes.size()) - 1 : special_row[k]; const int index = next_originals.offset(row) + count[row]++; next_originals.get_elist()[index] = get_edge(i); next_pos[i] = ~index; } } swap(g, next_graph); swap(originals, next_originals); pos.swap(next_pos); merged_pos.swap(next_merged_pos), merged_flow.swap(next_merged_flow); built = pos.size(); dirty = false; vc().swap(pending); } edge get_edge(int i) const { assert(0 <= i && i < num_of_edges()); if (i >= built) return pending[i - built]; if (pos[i] < 0) return originals.get_elist()[~pos[i]]; const auto &data = g.get_elist(); const auto &e = data[pos[i]], &r = data[e.rev]; return {r.to, e.to, Cap(e.cap + r.cap), r.cap}; } private: template Cap run_push_relabel(int s, int t, Cap flow_limit, ll budget, bool &finished, F resume) { const int source = n, N = n + 1, max_h = 2 * N; Cap return_capacity = flow_limit; vc excess(N, 0); vc height(N), iter(N), count(max_h + 1), head(max_h + 1, -1); vc next(N), que(N); vc active(N, false); int highest = -1; ll work = 0; ll spent = 0; bool interrupted = false; if constexpr (bounded) spent = 8LL * N; auto exhausted = [&] { if constexpr (bounded) if (spent >= budget) { const ll extension = resume(); if (extension == 0) { interrupted = true; return true; } budget = spent + extension; } return false; }; const ll global_budget = 2LL * g.get_elist().size() + 6LL * N; auto activate = [&](int v) { if (v == source || v == t || excess[v] == 0 || active[v]) return; assert(height[v] < max_h); next[v] = head[height[v]], head[height[v]] = v; active[v] = true; highest = max(highest, height[v]); }; auto rebuild_active = [&](int current) { if constexpr (bounded) spent += 4LL * N; fill(ALL(head), -1); fill(ALL(active), false); highest = -1; for (int v = 0; v < n; v++) if (v != current) activate(v); }; auto global_relabel = [&] { if constexpr (bounded) spent += 4LL * N; fill(ALL(height), max_h); fill(ALL(iter), 0); int front = 0, back = 0; height[t] = 0, que[back++] = t; while (front < back) { const int v = que[front++]; for (const auto &e : g[v]) { if constexpr (bounded) spent++; if (height[e.to] == max_h && g.get_elist()[e.rev].cap > 0) height[e.to] = height[v] + 1, que[back++] = e.to; } } front = back = 0; height[source] = N; if (height[s] == max_h && return_capacity > 0) height[s] = N + 1, que[back++] = s; while (front < back) { const int v = que[front++]; for (const auto &e : g[v]) { if constexpr (bounded) spent++; if (height[e.to] == max_h && g.get_elist()[e.rev].cap > 0) height[e.to] = height[v] + 1, que[back++] = e.to; } } fill(ALL(count), 0); for (int h : height) count[h]++; rebuild_active(-1); work = 0; }; excess[s] = flow_limit; global_relabel(); while (highest >= 0 && !exhausted()) { if constexpr (bounded) spent++; if (head[highest] == -1) { highest--; continue; } const int v = head[highest]; head[highest] = next[v], active[v] = false; const auto row = g[v]; while (excess[v] > 0 && !exhausted()) { int &i = iter[v]; const int degree = row.size(); if (i == degree + (v == s)) { const int old = height[v]; int h = max_h; if (v == s) { work++; if constexpr (bounded) spent++; if (return_capacity > 0) h = N + 1; } for (const auto &e : g[v]) { work++; if constexpr (bounded) spent++; if (e.cap > 0) h = min(h, height[e.to] + 1); } assert(h < max_h); height[v] = h, i = 0; count[old]--, count[h]++; if (old < N && count[old] == 0) { if constexpr (bounded) spent += n; for (int u = 0; u < n; u++) if (old < height[u] && height[u] < N) { count[height[u]]--; height[u] = N + 1, count[N + 1]++, iter[u] = 0; } rebuild_active(v); } continue; } work++; if constexpr (bounded) spent++; if (i == degree) { if (return_capacity > 0 && height[v] == N + 1) { const Cap f = min(excess[v], return_capacity); return_capacity -= f; excess[v] -= f, excess[source] += f; } else i++; continue; } auto &e = row[i]; if (e.cap > 0 && height[v] == height[e.to] + 1) { const Cap f = min(excess[v], e.cap); e.cap -= f, g.get_elist()[e.rev].cap += f; excess[v] -= f, excess[e.to] += f; activate(e.to); } else i++; } if (excess[v] > 0) break; if (work >= global_budget) global_relabel(); } finished = !interrupted; const Cap result = excess[t]; return result; } public: Cap flow(int s, int t, Cap flow_limit = numeric_limits::max()) { check_flow(s, t, flow_limit); if (flow_limit == 0) return 0; build(); flow_limit = flow_bound(s, t, flow_limit); if (flow_limit == 0) return 0; bool finished; ll budget = n + ll(g.get_elist().size()); if (n < 128 || g.get_elist().size() < 1024) return finish_flow(run_dinic(s, t, flow_limit, budget, finished)); int structure = 0; Cap total = run_dinic(s, t, flow_limit, budget, finished, &structure); if (finished || total == flow_limit) return finish_flow(total); if (structure == 0) { structure = total >= flow_limit - total ? 8 : dinic_budget_factor(s, t, total); if (structure > 1) { budget *= structure; total += run_dinic(s, t, flow_limit - total, budget, finished); if (finished || total == flow_limit) return finish_flow(total); } } vc saved; saved.reserve(g.get_elist().size()); for (const auto &e : g.get_elist()) saved.eb(e.cap); CSR dinic_graph; Cap extra = 0; auto resume = [&]() -> ll { if (dinic_graph.template size() == 0) { dinic_graph = g; int index = 0; for (auto &e : dinic_graph.get_elist()) e.cap = saved[index++]; vc().swap(saved); } swap(g, dinic_graph); budget *= 2; bool done; extra += run_dinic(s, t, flow_limit - total - extra, budget, done); swap(g, dinic_graph); return done ? 0 : 8 * budget; }; const Cap f = run_push_relabel(s, t, flow_limit - total, 8 * budget, finished, resume); if (finished) return finish_flow(total + f); swap(g, dinic_graph); return finish_flow(total + extra); } vc min_cut(int s) const { assert(0 <= s && s < n); build(); vc seen(n, false); vc que(n); int front = 0, back = 0; seen[s] = true, que[back++] = s; while (front < back) { const int v = que[front++]; for (const auto &e : g[v]) if (e.cap > 0 && !seen[e.to]) seen[e.to] = true, que[back++] = e.to; } return seen; } }; void main2() { LL(N, M, S, T); S--, T--; VEC(tlll, M, UVW); offset(UVW, tlll{-1, -1, 0}); UNZIP(UVW, U, V, W); MaxFlow mf(N); fec([ u, v, w ] : UVW) mf.add_edge(u, v, w); PRINT(mf.flow(S, T)); } void test() { } // https://github.com/miscalculation53/library/tree/wip/template/template_main.hpp template struct Main { Main() { cauto CERR = [](string val, string color) { string s = "\033[" + color + "m" + val + "\033[m"; }; CERR("\n[FAST_IO]\n\n", "32"); cout << fixed << setprecision(20); init(); CERR("\n[SINGLE_TESTCASE]\n\n", "36"); main2(); } }; Main main_dummy; } int main() {}