#define INF 4'000'000'000'000'000'037LL #include using namespace std; namespace { using ll = long long; using uint = unsigned int; using ull = unsigned long long; using pll = pair; #define vc vector template using vvc = vc>; using vpll = vc; #ifdef __SIZEOF_INT128__ using i128 = __int128_t; using u128 = __uint128_t; #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 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 repi(...) overload4(__VA_ARGS__, repi3, repi2, repi1)(__VA_ARGS__) #define fec(...) for (cauto &__VA_ARGS__) #define fem(...) for (auto &__VA_ARGS__) template inline bool chmin(T &a, U b) { return a > b ? a = b, true : false; } template inline constexpr T divfloor(U a, V b) { return T(a) / T(b) - (T(a) % T(b) && (T(a) ^ T(b)) < 0); } template 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); } #define ALL(a) (a).begin(), (a).end() #define eb emplace_back #define LMD(x, fx) ([&](auto x) { return fx; }) 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 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...); } #if __cplusplus < 202002L #else #endif template void unique(V &v) { v.erase(std::unique(ALL(v)), v.end()); } template void rotate(V &v, U k) { const U n = v.size(); k = (k % n + n) % n; std::rotate(v.begin(), v.begin() + k, v.end()); } 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; } }; #else #endif 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_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); } #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 #define dump(...) namespace fastio { 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; } template void rd1_integer(T &x) { if (pil + 100 > pir) load(); char c; do c = ibuf[pil++]; while (c < '-'); bool minus = 0; if constexpr (is_signed::value || is_same_v) { if (c == '-') { minus = 1, c = ibuf[pil++]; } } x = 0; while ('0' <= c) { x = x * 10 + (c & 15), c = ibuf[pil++]; } if constexpr (is_signed::value || is_same_v) { if (minus) x = -x; } } void rd1(ll &x) { rd1_integer(x); } template void rd1(pair &p) { return rd1(p.first), rd1(p.second); } template void rd1_tuple(T &t) { if constexpr (N < std::tuple_size::value) { auto &x = std::get(t); rd1(x); rd1_tuple(t); } } template void rd1(tuple &tpl) { rd1_tuple(tpl); } template void rd1(array &x) { for (auto &d: x) rd1(d); } template void rd1(vc &x) { for (auto &d: x) rd1(d); } void read() {} template void read(H &h, T &... t) { rd1(h), read(t...); } void wt1(const char c) { if (por == SIZ) flush(); obuf[por++] = c; } template void wt1_integer(T x) { if (por > SIZ - 100) flush(); if (x < 0) { obuf[por++] = '-', x = -x; } int outi; for (outi = 96; x >= 10000; outi -= 4) { memcpy(out + outi, pre.num[x % 10000], 4); x /= 10000; } if (x >= 1000) { memcpy(obuf + por, pre.num[x], 4); por += 4; } else if (x >= 100) { memcpy(obuf + por, pre.num[x] + 1, 3); por += 3; } else if (x >= 10) { int q = (x * 103) >> 10; obuf[por] = q | '0'; obuf[por + 1] = (x - q * 10) | '0'; por += 2; } else obuf[por++] = x | '0'; memcpy(obuf + por, out + outi + 4, 96 - outi); por += 96 - outi; } void wt1(int x) { wt1_integer(x); } template , int> = 0> void wt1(T x) { wt1_integer(x); } template void wt1(const pair &val) { wt1(val.first); wt1(' '); wt1(val.second); } template void wt1_tuple(const T &t) { if constexpr (N < std::tuple_size::value) { if constexpr (N > 0) { wt1(' '); } const auto x = std::get(t); wt1(x); wt1_tuple(t); } } template void wt1(const tuple &tpl) { wt1_tuple(tpl); } template void wt1(const array &val) { auto n = val.size(); for (size_t i = 0; i < n; i++) { if (i) wt1(' '); wt1(val[i]); } } template void wt1(const vector &val) { auto n = val.size(); for (size_t i = 0; i < n; i++) { if (i) wt1(' '); wt1(val[i]); } } void print() { wt1('\n'); } template void print(Head &&head, Tail &&... tail) { wt1(head); if (sizeof...(Tail)) wt1(' '); print(forward(tail)...); } } // namespace fastio 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); } 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 LL(...) IN(ll, __VA_ARGS__) #define READVEC(...) internal::READVECnodump(__VA_ARGS__); dump(__VA_ARGS__) #define VEC(T, n, ...) vc __VA_ARGS__; READVEC(n, __VA_ARGS__) #define PRINT fastio::print 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 vc> top(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> top(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> top(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 top(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 top(const tuple...> &tv) { size_t n = get<0>(tv).size(); apply([&](auto &...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; } mt19937_64 mt; template struct Binomial { private: static decltype(T::mod()) mod; static vc fac_, finv_, inv_; public: static void reserve(int n) { if (mod != T::mod()) { mod = T::mod(); fac_ = {1, 1}, finv_ = {1, 1}, inv_ = {0, 1}; } int i = fac_.size(); chmin(n, T::mod() - 1); if (n < i) return; fac_.resize(n + 1), finv_.resize(n + 1), inv_.resize(n + 1); for (; i <= n; i++) { 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 inv(T n) { assert(n != 0); reserve(n.val()); return inv_[n.val()]; } }; template decltype(T::mod()) Binomial::mod{}; template vc Binomial::fac_{}; template vc Binomial::finv_{}; template vc Binomial::inv_{}; namespace internal { constexpr ll powmod32_constexpr(ll x, ll n, int m) { if (m == 1) return 0; uint _m = (uint)m; ull r = 1; ull y = safemod(x, m); while (n) { if (n & 1) r = (r * y) % _m; y = (y * y) % _m; n >>= 1; } return r; } constexpr bool isprime32_constexpr(int 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; constexpr ll bases[3] = {2, 7, 61}; for (ll a : bases) { ll t = d; ll y = powmod32_constexpr(a, t, n); while (t != n - 1 && y != 1 && y != n - 1) { y = y * y % n; t <<= 1; } if (y != n - 1 && t % 2 == 0) return false; } return true; } template constexpr bool isprime32 = isprime32_constexpr(n); 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)); } }; } namespace internal { #define REF static_cast(*this) #define CREF static_cast(*this) #define VAL *static_cast(this) template struct modint_base { mint &operator+=(const mint &rhs) { mint &self = REF; self._v += rhs._v; if (self._v >= self.umod()) self._v -= self.umod(); return self; } mint &operator-=(const mint &rhs) { mint &self = REF; self._v -= rhs._v; if (self._v >= self.umod()) self._v += self.umod(); return self; } mint &operator/=(const mint &rhs) { mint &self = REF; return self = self * rhs.inv(); } mint &operator++() { mint &self = REF; self._v++; if (self._v == self.umod()) self._v = 0; return self; } mint &operator--() { mint &self = REF; if (self._v == 0) self._v = self.umod(); self._v--; return self; } mint operator++(int) { mint res = VAL; ++REF; return res; } mint operator--(int) { mint res = VAL; --REF; return res; } mint operator+() const { return VAL; } mint operator-() const { return mint() - VAL; } mint pow(ll n) const { assert(n >= 0); mint x = VAL, r = 1; while (n) { if (n & 1) r *= x; x *= x; n >>= 1; } return r; } 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 mint(lhs).eq(rhs); } friend bool operator!=(const mint &lhs, const mint &rhs) { return mint(lhs).neq(rhs); } private: bool eq(const mint &rhs) { return REF._v == rhs._v; } bool neq(const mint &rhs) { return REF._v != rhs._v; } }; } template , T>, int> = 0> void rd1(T &x) { ll a; fastio::rd1(a); x = a; } template , T>, int> = 0> void wt1(const T &x) { fastio::wt1(x.val()); } template constexpr tuple extgcd(const T &a, const 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}; } template struct static_modint : internal::modint_base> { using mint = static_modint; private: friend struct internal::modint_base>; uint _v; static constexpr uint umod() { return m; } static constexpr bool prime = internal::isprime32; public: static constexpr int mod() { return m; } static mint raw(int v) { mint x; x._v = v; return x; } static_modint() : _v(0) {} template static_modint(T v) { if constexpr (is_signed_v) { ll x = (ll)(v % (ll)(umod())); if (x < 0) x += umod(); _v = (uint)x; } else if constexpr (is_unsigned_v) { _v = (uint)(v % umod()); } else { static_assert(is_signed_v || is_unsigned_v, "Unsupported Type"); } } int val() const { return (int)_v; } mint& operator*=(const mint &rhs) { ull z = _v; z *= rhs._v; _v = (uint)(z % umod()); return *this; } mint inv() const { if (prime) { assert(_v != 0); return CREF.pow(umod() - 2); } else { auto [g, x, y] = extgcd(_v, m); assert(g == 1); return x; } } }; template struct dynamic_modint : internal::modint_base> { using mint = dynamic_modint; private: friend struct internal::modint_base>; uint _v; static internal::barrett32 bt; static uint umod() { return bt.umod(); } public: static int mod() { return (int)(bt.umod()); } static mint raw(int v) { mint x; x._v = v; return x; } dynamic_modint() : _v(0) {} template dynamic_modint(T v) { if constexpr (is_signed_v) { ll x = (ll)(v % (ll)(umod())); if (x < 0) x += umod(); _v = (uint)x; } else if constexpr (is_unsigned_v) { _v = (uint)(v % umod()); } else { static_assert(is_signed_v || is_unsigned_v, "Unsupported Type"); } } int val() const { return (int)_v; } mint& operator*=(const mint &rhs) { _v = bt.mul(_v, rhs._v); return *this; } mint inv() const { auto [g, x, y] = extgcd(_v, mod()); assert(g == 1); return x; } }; template internal::barrett32 dynamic_modint::bt(998244353); using modint998244353 = static_modint<998244353>; template struct is_static_modint : false_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 inline constexpr bool is_dynamic_modint_v = is_dynamic_modint::value; template constexpr pair crt_mod_constexpr(const array &rs, const array &ms) { 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}; } 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 (powmod32_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 std::array iroot; // root[i] * iroot[i] == 1 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 constexpr (is_static_modint_v) { if constexpr (!internal::isprime32) return false; static constexpr int rank2 = countr_zero(mint::mod() - 1); return n <= (1 << rank2); } else return false; } template void ntt(vc> &) {} template void intt(vc> &) {} template void ntt(vc> &a) { using mint = static_modint; int n = int(a.size()); int h = countr_zero((unsigned int)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 = static_modint; int n = int(a.size()); int h = countr_zero((unsigned int)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_modint; 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_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_modint; auto c = convolution(vc(ALL(a)), vc(ALL(b))); vc c_(c.size()); repi(i, c.size()) c_[i] = c[i].val(); return c_; } namespace internal { constexpr ll powmod64_constexpr(ll x, ll n, ll m) { if (m == 1) return 0; ull _m = (ull)m; ull r = 1; ull y = safemod(x, m); while (n) { u128 y128(y); if (n & 1) r = (y128 * r) % _m; y = (y128 * y) % _m; n >>= 1; } return r; } constexpr bool isprime64_constexpr(ll n) { if (n <= INT_MAX) return isprime32_constexpr(n); if (n % 2 == 0) return false; ll d = n - 1; while (d % 2 == 0) d /= 2; constexpr ll bases[7] = {2, 325, 9375, 28178, 450775, 9780504, 1795265022}; for (ll a : bases) { ll t = d; ll y = powmod64_constexpr(a, t, n); while (t != n - 1 && y != 1 && y != n - 1) { y = (u128(y) * y) % n; t <<= 1; } if (y != n - 1 && t % 2 == 0) return false; } return true; } template constexpr bool isprime64 = isprime64_constexpr(n); 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; // m == 2^b * mx, mx is odd imx = inv64(mx); d = powmod64_constexpr((mx + 1) / 2, b, mx); // 2^{-b} mod mx u128 sq = -u128(mx) % mx; // 2^128 mod mx q = (1 + (((sq - 1) * d) << b)) % m; } ull umod() const { return m; } ull reduce(u128 x) const { ull p = x & MASK(b); // x mod 2^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 static_modint64 : internal::modint_base> { using mint = static_modint64; private: friend struct internal::modint_base>; ull _v; static constexpr ull umod() { return m; } static constexpr bool prime = internal::isprime64; public: static constexpr ll mod() { return m; } static mint raw(ll v) { mint x; x._v = v; return x; } static_modint64() : _v(0) {} template static_modint64(T v) { if constexpr (is_unsigned_v) { _v = (ull)(v % umod()); } else { ll x = (ll)(v % (ll)(umod())); if (x < 0) x += umod(); _v = (ull)x; } } ll val() const { return (ll)_v; } mint& operator*=(const mint &rhs) { u128 z = _v; z *= rhs._v; _v = (ull)(z % umod()); return *this; } mint inv() const { if (prime) { assert(_v != 0); return CREF.pow(umod() - 2); } else { auto [g, x, y] = extgcd(_v, m); assert(g == 1); return x; } } }; template struct dynamic_modint64_odd : internal::modint_base> { using mint = dynamic_modint64_odd; private: friend struct internal::modint_base>; ull _v; // montgomery expression static internal::montgomery64odd mg; static ull umod() { return mg.umod(); } public: static ll mod() { return (ll)(mg.umod()); } dynamic_modint64_odd() : _v(0) {} dynamic_modint64_odd(i128 v) { _v = mg.inv_reduce(v); } ll val() const { return (ll)mg.reduce(_v); } mint& operator*=(const mint &rhs) { _v = mg.reduce(u128(_v) * rhs._v); return *this; } mint inv() const { auto [g, x, y] = extgcd(val(), mod()); assert(g == 1); return x; } }; template internal::montgomery64odd dynamic_modint64_odd::mg((1LL << 61) - 1); template struct dynamic_modint64 : internal::modint_base> { using mint = dynamic_modint64; private: friend struct internal::modint_base>; ull _v; // montgomery expression static internal::montgomery64 mg; static ull umod() { return mg.umod(); } public: static ll mod() { return (ll)(mg.umod()); } dynamic_modint64() : _v(0) {} dynamic_modint64(i128 v) { _v = mg.inv_reduce(v); } ll val() const { return (ll)mg.reduce(_v); } mint& operator*=(const mint &rhs) { _v = mg.reduce(u128(_v) * rhs._v); return *this; } mint inv() const { auto [g, x, y] = extgcd(val(), mod()); assert(g == 1); return x; } }; template internal::montgomery64 dynamic_modint64::mg((1LL << 61) - 1); namespace internal { }; // namespace internal template struct FormalPowerSeries : vc { using F = FormalPowerSeries; using vc::vc; using vc::operator=; using vc::size; using vc::empty; using vc::back; using vc::pop_back; using vc::begin; using vc::resize; using vc::front; FormalPowerSeries(const vc &f) : vc(f) {} int sz() const { return size(); } void shrink() { while (!empty() && back() == 0) pop_back(); } mint get(int i) const { return 0 <= i && i < sz() ? (*this)[i] : 0; } F pre(int len) const { assert(len >= 0); return F(begin(), begin() + min(sz(), len)); } F rev(int d = -1) const { F res(*this); if (d >= 0) res.resize(d); reverse(ALL(res)); return res; } int cnt_nz() const { return count_if(ALL(*this), LMD(x, x != 0)); } tuple nz_front() const { repi(i, sz()) if ((*this)[i] != 0) return {true, i, (*this)[i]}; return {false, -1, 0}; } vc> nz() const { vc> res; repi(i, sz()) if ((*this)[i] != 0) res.eb(i, (*this)[i]); return res; } F operator-() const { F res(*this); fem(a : res) a = -a; return res; } F &operator*=(const mint &k) { fem(a : *this) a *= k; return *this; } F operator*(const mint &k) const { return F(*this) *= k; } friend F operator*(const mint &k, const F &f) { return f * k; } F &operator/=(const mint &k) { *this *= k.inv(); return *this; } F operator/(const mint &k) const { return F(*this) /= k; } F &operator+=(const F &g) { const int n = size(), m = g.size(); resize(max(n, m)); repi(i, m)(*this)[i] += g[i]; return *this; } F operator+(const F &g) const { return F(*this) += g; } F &operator-=(const F &g) { const int n = size(), m = g.size(); resize(max(n, m)); repi(i, m)(*this)[i] -= g[i]; return *this; } F operator-(const F &g) const { return F(*this) -= g; } F &operator*=(const F &g) { return *this = *this * g; } F operator*(const F &g) const { return convolution(*this, g); } F div_sparse_destructive(const F &g, int d = -1) { assert(g.get(0) != 0); if (d < 0) d = sz(); mint iv = g.front().inv(); auto gnz = g.nz(); resize(d); repi(i, d) { fec([j, b] : gnz) { if (j == 0) continue; if (j > i) break; (*this)[i] -= (*this)[i - j] * b; } (*this)[i] *= iv; } return pre(d); } F div_sparse(const F &g, int d = -1) const { return F(*this).div_sparse_destructive(g, d); } F inv(int d = -1) const { assert(get(0) != 0); if (d < 0) d = sz(); if (cnt_nz() <= 200) return F{1}.div_sparse(*this, d); F f, g2, g{front().inv()}; for (int m = 1; m < d; m *= 2) { if (ntt_ok(2 * m)) { f = pre(2 * m), g2 = F(g); f.resize(2 * m), ntt(f); g2.resize(2 * m), ntt(g2); repi(i, 2 * m) f[i] *= g2[i]; intt(f); f >>= m; f.resize(2 * m), ntt(f); repi(i, 2 * m) f[i] *= g2[i]; intt(f); mint iz = mint(2 * m).inv(); iz *= -iz; repi(i, m) f[i] *= iz; g.insert(g.end(), f.begin(), f.begin() + m); } else g = (g * mint(2) - g * g * pre(2 * m)).pre(2 * m); } return g.pre(d); } F &operator/=(const F &g) { if (cnt_nz() <= 200) { div_sparse_destructive(g); return *this; } *this *= g.inv(); return *this; } F operator/(const F &g) const { return F(*this) /= g; } F div_poly(const F &g) const { const int n = sz() - g.sz() + 1; if (n <= 0) return {}; return (rev().pre(n) * g.rev().inv(n)).pre(n).rev(); } pair divmod(const F &g) const { F q = div_poly(g); F r = *this - q * g; r.shrink(); return {q, r}; } F operator%(const F &g) const { return divmod(g).second; } F &operator%=(const F &g) { return *this = *this % g; } F circular_mod(int n) const { F res(n); repi(i, sz()) res[i % n] += (*this)[i]; return res; } F operator<<(int k) const { F res(sz() + k); repi(i, sz()) res[i + k] = (*this)[i]; return res; } F operator>>(int k) const { F res(max(0, sz() - k)); repi(i, sz() - k) res[i] = (*this)[i + k]; return res; } F &operator<<=(int k) { return *this = *this << k; } F &operator>>=(int k) { return *this = *this >> k; } F diff() const { F res(max(0, sz() - 1)); repi(i, 1, size()) res[i - 1] = (*this)[i] * i; return res; } F integ() const { F res(sz() + 1); repi(i, size()) res[i + 1] = (*this)[i] * Binomial::inv(i + 1); return res; } F log(int d = -1) const { assert(get(0) == 1); if (d < 0) d = sz(); F f = pre(d); return (f.diff() / f).pre(d - 1).integ(); } static F diff_eq(const F &a, const F &b, int d) { assert(a.get(0) == 1); assert(d >= 0); if (d == 0) return {}; F f(d); f[0] = 1; auto anz = a.nz(), bnz = b.nz(); repi(k, d - 1) { fec([i, ai] : anz) { if (0 <= k - i + 1) f[k + 1] -= ai * (k - i + 1) * f[k - i + 1]; } fec([j, bj] : bnz) { if (0 <= k - j && k - j < k + 1) f[k + 1] -= bj * f[k - j]; } f[k + 1] *= Binomial::inv(k + 1); } return f; } F exp_sparse(int d = -1) const { assert(get(0) == 0); if (d < 0) d = sz(); return diff_eq(F{1}, -diff(), d); } F pow_sparse(ll k, int d = -1) const { if (d < 0) d = sz(); auto [exi, d0, a0] = nz_front(); if (!exi) { F res(d); if (k == 0 && d > 0) res[0] = 1; return res; } mint ia0 = a0.inv(); F f = ((*this) >> d0) * ia0; if (k >= 0) { F g = diff_eq(f, -k * f.diff(), d - mul_limited(d0, k, d)); F h = (g * a0.pow(k)) << mul_limited(d0, k, d); return h.pre(d); } else { F g = diff_eq(f, -k * f.diff(), d + (d0 * (-k))); F h = (g * ia0.pow(-k)) >> (d0 * (-k)); return h.pre(d); } } F exp(int d = -1) const { assert(get(0) == 0); if (d < 0) d = sz(); if (ntt_ok(2 * d)) { if (cnt_nz() <= 320) return exp_sparse(d); F f{1}, g{1}; F f2, g2, f3, q, s, h, u; g2 = {0}; for (int m = 1; m < d; m *= 2) { mint im = mint(m).inv(), i2m = mint(2 * m).inv(); f2 = f, f2.resize(2 * m), ntt(f2); f3 = f, ntt(f3); repi(i, m) f3.at(i) *= g2.at(i); intt(f3); f3 >>= m / 2; f3.resize(m), ntt(f3); repi(i, m) f3.at(i) *= g2.at(i); intt(f3); repi(i, m / 2) f3.at(i) *= -im * im; g.insert(g.end(), f3.begin(), f3.begin() + m / 2); g2 = g, g2.resize(2 * m), ntt(g2); q = diff(), q.resize(2 * m), fill(q.begin() + m - 1, q.end(), 0); ntt(q); repi(i, 2 * m) q.at(i) *= f2.at(i); intt(q); q = q.circular_mod(m); repi(i, m) q.at(i) *= i2m; q.resize(m + 1); s = ((f.diff() - q) << 1).circular_mod(m); s.resize(2 * m), ntt(s); repi(i, 2 * m) s.at(i) *= g2.at(i); intt(s); repi(i, m) s.at(i) *= i2m; s.resize(m); h = *this, h.resize(2 * m), s.resize(2 * m); u = (h - (s << (m - 1)).integ()) >> m; ntt(u); repi(i, 2 * m) u.at(i) *= f2.at(i); intt(u); repi(i, m) u.at(i) *= i2m; u.resize(m); f.insert(f.end(), u.begin(), u.end()); } return f.pre(d); } else { if (cnt_nz() <= 3000) return exp_sparse(d); F f{1}; for (int m = 1; m < d; m *= 2) { f = (f * (pre(2 * m) + F{1} - f.log(2 * m))).pre(2 * m); } return f.pre(d); } } F pow(ll k, int d = -1) const { if (ntt_ok(2 * d)) { if (cnt_nz() <= 100) return pow_sparse(k, d); } else { if (cnt_nz() <= 1300) return pow_sparse(k, d); } if (d < 0) d = sz(); if (k == 0) { F res(d); res[0] = 1; return res; } repi(i, sz()) { if ((*this)[i] != 0) { mint iv = (*this)[i].inv(); F res = (((*this * iv) >> i).log(d) * mint(k)).exp(d); res *= (*this)[i].pow(k); res = (res << (i * k)).pre(d); if (res.sz() < d) res.resize(d); return res; } if (mul_limited(i + 1, k, d) >= d) return F(d); } return F(d); } }; using mint = modint998244353; using bi = Binomial; using fps = FormalPowerSeries; template pair, FormalPowerSeries> rational_plus ( const pair, FormalPowerSeries> &f, const pair, FormalPowerSeries> &g ) { using F = FormalPowerSeries; cauto &[p_, q_] = f; cauto &[r_, s_] = g; if (min({p_.cnt_nz(), q_.cnt_nz(), r_.cnt_nz(), s_.cnt_nz()}) <= 60) return {p_ * s_ + q_ * r_, q_ * s_}; F p = p_, q = q_, r = r_, s = s_; const int k = p.size(), l = q.size(), m = r.size(), n = s.size(); const int z = bit_ceil(max({k + n - 1, l + m - 1, l + n - 1})); p.resize(z), q.resize(z), r.resize(z), s.resize(z); ntt(p), ntt(q), ntt(r), ntt(s); F a(z), b(z), c(z); repi(i, z) a[i] = p[i] * s[i], b[i] = q[i] * r[i], c[i] = q[i] * s[i]; intt(a), intt(b), intt(c); mint iz = mint(z).inv(); repi(i, z) a[i] *= iz, b[i] *= iz, c[i] *= iz; return {a.pre(k + n - 1) + b.pre(l + m - 1), c.pre(l + n - 1)}; } template pair, FormalPowerSeries> rational_sum (const vc, FormalPowerSeries>> &fs, int d = -1) { using F = FormalPowerSeries; using R = pair; auto dc = [&](auto dc, int l, int r) -> R { if (r - l == 0) return {{1}, {1}}; if (r - l == 1) return fs[l]; const int m = (l + r) / 2; R res = rational_plus(dc(dc, l, m), dc(dc, m, r)); if (d < 0) return res; else return {res.first.pre(d), res.second.pre(d)}; }; return dc(dc, 0, fs.size()); } void init() {} void main2() { LL(N, S); VEC(mint, N, P); fem(p : P) p /= S; vc> fs(N); rep(i, N) fs.at(i) = {{0, P.at(i) * P.at(i)}, {1, P.at(i)}}; auto f = rational_sum(fs); dump(fs, f); mint ans = 0; rep(k, N + 1) { ans += (k + 1) * bi::fac(k) * f.first.at(k); dump(k, f.first.at(k), bi::fac(k) * f.first.at(k)); } 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_IO]\n\n", "32"); cout << fixed << setprecision(20); test(); init(); CERR("\n[SINGLE_TESTCASE]\n\n", "36"); main2(); } }; Main main_dummy; } int main() {}