// #pragma GCC target("avx2") // #pragma GCC optimize("O3") // #pragma GCC optimize("unroll-loops") #include using namespace std; // #define INTERACTIVE namespace templates { // type using ll = long long; using ull = unsigned long long; using Pii = pair; using Pil = pair; using Pli = pair; using Pll = pair; template using pq = priority_queue; template using qp = priority_queue, greater>; // clang-format off #define vec(T, A, ...) vector A(__VA_ARGS__); #define vvec(T, A, h, ...) vector> A(h, vector(__VA_ARGS__)); #define vvvec(T, A, h1, h2, ...) vector>> A(h1, vector>(h2, vector(__VA_ARGS__))); // clang-format on // for loop #define fori1(a) for (ll _ = 0; _ < (a); _++) #define fori2(i, a) for (ll i = 0; i < (a); i++) #define fori3(i, a, b) for (ll i = (a); i < (b); i++) #define fori4(i, a, b, c) for (ll i = (a); ((c) > 0 || i > (b)) && ((c) < 0 || i < (b)); i += (c)) #define overload4(a, b, c, d, e, ...) e #define fori(...) overload4(__VA_ARGS__, fori4, fori3, fori2, fori1)(__VA_ARGS__) // declare and input // clang-format off #define INT(...) int __VA_ARGS__; inp(__VA_ARGS__); #define LL(...) ll __VA_ARGS__; inp(__VA_ARGS__); #define STRING(...) string __VA_ARGS__; inp(__VA_ARGS__); #define CHAR(...) char __VA_ARGS__; inp(__VA_ARGS__); #define DOUBLE(...) double __VA_ARGS__; STRING(str___); __VA_ARGS__ = stod(str___); #define VEC(T, A, n) vector A(n); inp(A); #define VVEC(T, A, n, m) vector> A(n, vector(m)); inp(A); // clang-format on // const value const ll MOD1 = 1000000007; const ll MOD9 = 998244353; const double PI = acos(-1); // other macro #if !defined(RIN__LOCAL) && !defined(INTERACTIVE) #define endl "\n" #endif #define spa ' ' #define len(A) ll(A.size()) #define all(A) begin(A), end(A) // function vector stoc(string &S) { int n = S.size(); vector ret(n); for (int i = 0; i < n; i++) ret[i] = S[i]; return ret; } string ctos(vector &S) { int n = S.size(); string ret = ""; for (int i = 0; i < n; i++) ret += S[i]; return ret; } template auto min(const T &a) { return *min_element(all(a)); } template auto max(const T &a) { return *max_element(all(a)); } template auto clamp(T &a, const S &l, const S &r) { return (a > r ? r : a < l ? l : a); } template inline bool chmax(T &a, const S &b) { return (a < b ? a = b, 1 : 0); } template inline bool chmin(T &a, const S &b) { return (a > b ? a = b, 1 : 0); } template inline bool chclamp(T &a, const S &l, const S &r) { auto b = clamp(a, l, r); return (a != b ? a = b, 1 : 0); } template T sum(vector &A) { T tot = 0; for (auto a : A) tot += a; return tot; } template vector compression(vector X) { sort(all(X)); X.erase(unique(all(X)), X.end()); return X; } // input and output namespace io { // vector template istream &operator>>(istream &is, vector &A) { for (auto &a : A) is >> a; return is; } template ostream &operator<<(ostream &os, vector &A) { for (size_t i = 0; i < A.size(); i++) { os << A[i]; if (i != A.size() - 1) os << ' '; } return os; } // vector> template istream &operator>>(istream &is, vector> &A) { for (auto &a : A) is >> a; return is; } template ostream &operator<<(ostream &os, vector> &A) { for (size_t i = 0; i < A.size(); i++) { os << A[i]; if (i != A.size() - 1) os << endl; } return os; } // pair template istream &operator>>(istream &is, pair &A) { is >> A.first >> A.second; return is; } template ostream &operator<<(ostream &os, pair &A) { os << A.first << ' ' << A.second; return os; } // vector> template istream &operator>>(istream &is, vector> &A) { for (size_t i = 0; i < A.size(); i++) { is >> A[i]; } return is; } template ostream &operator<<(ostream &os, vector> &A) { for (size_t i = 0; i < A.size(); i++) { os << A[i]; if (i != A.size() - 1) os << endl; } return os; } // tuple template struct TuplePrint { static ostream &print(ostream &os, const T &t) { TuplePrint::print(os, t); os << ' ' << get(t); return os; } }; template struct TuplePrint { static ostream &print(ostream &os, const T &t) { os << get<0>(t); return os; } }; template ostream &operator<<(ostream &os, const tuple &t) { TuplePrint::print(os, t); return os; } // io functions void FLUSH() { cout << flush; } void print() { cout << endl; } template void print(Head &&head, Tail &&...tail) { cout << head; if (sizeof...(Tail)) cout << spa; print(std::forward(tail)...); } template void prisep(vector &A, S sep) { int n = A.size(); for (int i = 0; i < n; i++) { cout << A[i]; if (i != n - 1) cout << sep; } cout << endl; } template void priend(T A, S end) { cout << A << end; } template void prispa(T A) { priend(A, spa); } template bool printif(bool f, T A, S B) { if (f) print(A); else print(B); return f; } template void inp(T &...a) { (cin >> ... >> a); } } // namespace io using namespace io; // read graph vector> read_edges(int n, int m, bool direct = false, int indexed = 1) { vector> edges(n, vector()); for (int i = 0; i < m; i++) { INT(u, v); u -= indexed; v -= indexed; edges[u].push_back(v); if (!direct) edges[v].push_back(u); } return edges; } vector> read_tree(int n, int indexed = 1) { return read_edges(n, n - 1, false, indexed); } template vector>> read_wedges(int n, int m, bool direct = false, int indexed = 1) { vector>> edges(n, vector>()); for (int i = 0; i < m; i++) { INT(u, v); T w; inp(w); u -= indexed; v -= indexed; edges[u].push_back({v, w}); if (!direct) edges[v].push_back({u, w}); } return edges; } template vector>> read_wtree(int n, int indexed = 1) { return read_wedges(n, n - 1, false, indexed); } // yes / no namespace yesno { // yes inline bool yes(bool f = true) { cout << (f ? "yes" : "no") << endl; return f; } inline bool Yes(bool f = true) { cout << (f ? "Yes" : "No") << endl; return f; } inline bool YES(bool f = true) { cout << (f ? "YES" : "NO") << endl; return f; } // no inline bool no(bool f = true) { cout << (!f ? "yes" : "no") << endl; return f; } inline bool No(bool f = true) { cout << (!f ? "Yes" : "No") << endl; return f; } inline bool NO(bool f = true) { cout << (!f ? "YES" : "NO") << endl; return f; } // possible inline bool possible(bool f = true) { cout << (f ? "possible" : "impossible") << endl; return f; } inline bool Possible(bool f = true) { cout << (f ? "Possible" : "Impossible") << endl; return f; } inline bool POSSIBLE(bool f = true) { cout << (f ? "POSSIBLE" : "IMPOSSIBLE") << endl; return f; } // impossible inline bool impossible(bool f = true) { cout << (!f ? "possible" : "impossible") << endl; return f; } inline bool Impossible(bool f = true) { cout << (!f ? "Possible" : "Impossible") << endl; return f; } inline bool IMPOSSIBLE(bool f = true) { cout << (!f ? "POSSIBLE" : "IMPOSSIBLE") << endl; return f; } // Alice Bob inline bool Alice(bool f = true) { cout << (f ? "Alice" : "Bob") << endl; return f; } inline bool Bob(bool f = true) { cout << (f ? "Bob" : "Alice") << endl; return f; } // Takahashi Aoki inline bool Takahashi(bool f = true) { cout << (f ? "Takahashi" : "Aoki") << endl; return f; } inline bool Aoki(bool f = true) { cout << (f ? "Aoki" : "Takahashi") << endl; return f; } } // namespace yesno using namespace yesno; } // namespace templates using namespace templates; template struct Modint { int x; Modint() : x(0) {} Modint(int64_t y) { if (y >= 0) x = y % MOD; else x = (y % MOD + MOD) % MOD; } Modint &operator+=(const Modint &p) { x += p.x; if (x >= MOD) x -= MOD; return *this; } Modint &operator-=(const Modint &p) { x -= p.x; if (x < 0) x += MOD; return *this; } Modint &operator*=(const Modint &p) { x = int(1LL * x * p.x % MOD); return *this; } Modint &operator/=(const Modint &p) { *this *= p.inverse(); return *this; } Modint &operator%=(const Modint &p) { assert(p.x == 0); return *this; } Modint operator-() const { return Modint(-x); } Modint &operator++() { x++; if (x == MOD) x = 0; return *this; } Modint &operator--() { if (x == 0) x = MOD; x--; return *this; } Modint operator++(int) { Modint result = *this; ++*this; return result; } Modint operator--(int) { Modint result = *this; --*this; return result; } friend Modint operator+(const Modint &lhs, const Modint &rhs) { return Modint(lhs) += rhs; } friend Modint operator-(const Modint &lhs, const Modint &rhs) { return Modint(lhs) -= rhs; } friend Modint operator*(const Modint &lhs, const Modint &rhs) { return Modint(lhs) *= rhs; } friend Modint operator/(const Modint &lhs, const Modint &rhs) { return Modint(lhs) /= rhs; } friend Modint operator%(const Modint &lhs, const Modint &rhs) { assert(rhs.x == 0); return Modint(lhs); } bool operator==(const Modint &p) const { return x == p.x; } bool operator!=(const Modint &p) const { return x != p.x; } bool operator<(const Modint &rhs) const { return x < rhs.x; } bool operator<=(const Modint &rhs) const { return x <= rhs.x; } bool operator>(const Modint &rhs) const { return x > rhs.x; } bool operator>=(const Modint &rhs) const { return x >= rhs.x; } Modint inverse() const { int a = x, b = MOD, u = 1, v = 0, t; while (b > 0) { t = a / b; a -= t * b; u -= t * v; std::swap(a, b); std::swap(u, v); } return Modint(u); } Modint pow(int64_t k) const { Modint ret(1); Modint y(x); while (k > 0) { if (k & 1) ret *= y; y *= y; k >>= 1; } return ret; } friend std::ostream &operator<<(std::ostream &os, const Modint &p) { return os << p.x; } friend std::istream &operator>>(std::istream &is, Modint &p) { int64_t y; is >> y; p = Modint(y); return (is); } static int get_mod() { return MOD; } }; struct Arbitrary_Modint { int x; static int MOD; static void set_mod(int mod) { MOD = mod; } Arbitrary_Modint() : x(0) {} Arbitrary_Modint(int64_t y) { if (y >= 0) x = y % MOD; else x = (y % MOD + MOD) % MOD; } Arbitrary_Modint &operator+=(const Arbitrary_Modint &p) { x += p.x; if (x >= MOD) x -= MOD; return *this; } Arbitrary_Modint &operator-=(const Arbitrary_Modint &p) { x -= p.x; if (x < 0) x += MOD; return *this; } Arbitrary_Modint &operator*=(const Arbitrary_Modint &p) { x = int(1LL * x * p.x % MOD); return *this; } Arbitrary_Modint &operator/=(const Arbitrary_Modint &p) { *this *= p.inverse(); return *this; } Arbitrary_Modint &operator%=(const Arbitrary_Modint &p) { assert(p.x == 0); return *this; } Arbitrary_Modint operator-() const { return Arbitrary_Modint(-x); } Arbitrary_Modint &operator++() { x++; if (x == MOD) x = 0; return *this; } Arbitrary_Modint &operator--() { if (x == 0) x = MOD; x--; return *this; } Arbitrary_Modint operator++(int) { Arbitrary_Modint result = *this; ++*this; return result; } Arbitrary_Modint operator--(int) { Arbitrary_Modint result = *this; --*this; return result; } friend Arbitrary_Modint operator+(const Arbitrary_Modint &lhs, const Arbitrary_Modint &rhs) { return Arbitrary_Modint(lhs) += rhs; } friend Arbitrary_Modint operator-(const Arbitrary_Modint &lhs, const Arbitrary_Modint &rhs) { return Arbitrary_Modint(lhs) -= rhs; } friend Arbitrary_Modint operator*(const Arbitrary_Modint &lhs, const Arbitrary_Modint &rhs) { return Arbitrary_Modint(lhs) *= rhs; } friend Arbitrary_Modint operator/(const Arbitrary_Modint &lhs, const Arbitrary_Modint &rhs) { return Arbitrary_Modint(lhs) /= rhs; } friend Arbitrary_Modint operator%(const Arbitrary_Modint &lhs, const Arbitrary_Modint &rhs) { assert(rhs.x == 0); return Arbitrary_Modint(lhs); } bool operator==(const Arbitrary_Modint &p) const { return x == p.x; } bool operator!=(const Arbitrary_Modint &p) const { return x != p.x; } bool operator<(const Arbitrary_Modint &rhs) { return x < rhs.x; } bool operator<=(const Arbitrary_Modint &rhs) { return x <= rhs.x; } bool operator>(const Arbitrary_Modint &rhs) { return x > rhs.x; } bool operator>=(const Arbitrary_Modint &rhs) { return x >= rhs.x; } Arbitrary_Modint inverse() const { int a = x, b = MOD, u = 1, v = 0, t; while (b > 0) { t = a / b; a -= t * b; u -= t * v; std::swap(a, b); std::swap(u, v); } return Arbitrary_Modint(u); } Arbitrary_Modint pow(int64_t k) const { Arbitrary_Modint ret(1); Arbitrary_Modint y(x); while (k > 0) { if (k & 1) ret *= y; y *= y; k >>= 1; } return ret; } friend std::ostream &operator<<(std::ostream &os, const Arbitrary_Modint &p) { return os << p.x; } friend std::istream &operator>>(std::istream &is, Arbitrary_Modint &p) { int64_t y; is >> y; p = Arbitrary_Modint(y); return (is); } static int get_mod() { return MOD; } }; int Arbitrary_Modint::MOD = 998244353; using modint9 = Modint<998244353>; using modint1 = Modint<1000000007>; using modint = Arbitrary_Modint; using mint = modint9; template struct Combination { int N; std::vector fact, invfact; Combination(int N) : N(N) { fact.resize(N + 1); invfact.resize(N + 1); fact[0] = 1; for (int i = 1; i <= N; i++) { fact[i] = fact[i - 1] * i; } invfact[N] = T(1) / fact[N]; for (int i = N - 1; i >= 0; i--) { invfact[i] = invfact[i + 1] * (i + 1); } } void extend(int n) { int le = fact.size(); fact.resize(n + 1); invfact.resize(n + 1); for (int i = le; i <= n; i++) { fact[i] = fact[i - 1] * i; } invfact[n] = T(1) / fact[n]; for (int i = n - 1; i >= le; i--) { invfact[i] = invfact[i + 1] * (i + 1); } } T nCk(int n, int k) { if (k > n || k < 0) return T(0); if (n >= int(fact.size())) extend(n); return fact[n] * invfact[k] * invfact[n - k]; } T nPk(int n, int k) { if (k > n || k < 0) return T(0); if (n >= int(fact.size())) extend(n); return fact[n] * invfact[n - k]; } T nHk(int n, int k) { if (n == 0 && k == 0) return T(1); return nCk(n + k - 1, k); } T Catalan(int n) { return nCk(2 * n, n) - nCk(2 * n, n + 1); } // n 個の +1, m 個の -1, 累積和が常にk以下 T Catalan(int n, int m, int k) { if (n > m + k || k < 0) return T(0); else return nCk(n + m, n) - nCk(n + m, m + k + 1); } T Narayana(int n, int k) { return nCk(n, k) * nCk(n, k - 1) / n; } T inv(int n) { assert(n >= 1); if (n >= int(fact.size())) extend(n); return invfact[n] * fact[n - 1]; } }; template struct NumberTheoreticTransform { static std::vector roots, iroots, rate3, irate3; static int max_base; NumberTheoreticTransform() = default; static void init() { if (!roots.empty()) return; const unsigned mod = mint::get_mod(); auto tmp = mod - 1; max_base = 0; while (tmp % 2 == 0) { tmp >>= 1; max_base++; } mint root = 2; while (root.pow((mod - 1) >> 1) == 1) root++; roots.resize(max_base + 1); iroots.resize(max_base + 1); rate3.resize(max_base + 1); irate3.resize(max_base + 1); roots[max_base] = root.pow((mod - 1) >> max_base); iroots[max_base] = mint(1) / roots[max_base]; for (int i = max_base - 1; i >= 0; i--) { roots[i] = roots[i + 1] * roots[i + 1]; iroots[i] = iroots[i + 1] * iroots[i + 1]; } mint prod = 1, iprod = 1; for (int i = 0; i <= max_base - 3; i++) { rate3[i] = roots[i + 3] * prod; irate3[i] = iroots[i + 3] * iprod; prod *= iroots[i + 3]; iprod *= roots[i + 3]; } } static void ntt(std::vector &A) { init(); int n = int(A.size()); int h = __builtin_ctz(n); int le = 0; mint imag = roots[2]; if (h & 1) { int p = 1 << (h - 1); for (int i = 0; i < p; i++) { auto r = A[i + p]; A[i + p] = A[i] - r; A[i] += r; } le++; } for (; le + 1 < h; le += 2) { int p = 1 << (h - le - 2); for (int i = 0; i < p; i++) { auto a0 = A[i]; auto a1 = A[i + p]; auto a2 = A[i + 2 * p]; auto a3 = A[i + 3 * p]; auto a1na3imag = (a1 - a3) * imag; A[i] = a0 + a2 + a1 + a3; A[i + p] = a0 + a2 - (a1 + a3); A[i + 2 * p] = a0 - a2 + a1na3imag; A[i + 3 * p] = a0 - a2 - a1na3imag; } mint rot = rate3[0]; for (int s = 1; s < (1 << le); s++) { int offset = s << (h - le); mint rot2 = rot * rot; mint rot3 = rot2 * rot; for (int i = 0; i < p; i++) { auto a0 = A[i + offset]; auto a1 = A[i + offset + p] * rot; auto a2 = A[i + offset + 2 * p] * rot2; auto a3 = A[i + offset + 3 * p] * rot3; auto a1na3imag = (a1 - a3) * imag; A[i + offset] = a0 + a2 + a1 + a3; A[i + offset + p] = a0 + a2 - (a1 + a3); A[i + offset + 2 * p] = a0 - a2 + a1na3imag; A[i + offset + 3 * p] = a0 - a2 - a1na3imag; } rot *= rate3[__builtin_ctz(~s)]; } } } static void intt(std::vector &A, bool f = true) { init(); int n = int(A.size()); int h = __builtin_ctz(n); int le = h; mint iimag = iroots[2]; for (; le > 1; le -= 2) { int p = 1 << (h - le); for (int i = 0; i < p; i++) { auto a0 = A[i]; auto a1 = A[i + p]; auto a2 = A[i + 2 * p]; auto a3 = A[i + 3 * p]; auto a2na3iimag = (a2 - a3) * iimag; A[i] = a0 + a1 + a2 + a3; A[i + p] = a0 - a1 + a2na3iimag; A[i + 2 * p] = a0 + a1 - (a2 + a3); A[i + 3 * p] = a0 - a1 - a2na3iimag; } mint irot = irate3[0]; for (int s = 1; s < (1 << (le - 2)); s++) { int offset = s << (h - le + 2); mint irot2 = irot * irot; mint irot3 = irot2 * irot; for (int i = 0; i < p; i++) { auto a0 = A[i + offset]; auto a1 = A[i + offset + p]; auto a2 = A[i + offset + 2 * p]; auto a3 = A[i + offset + 3 * p]; auto a2na3iimag = (a2 - a3) * iimag; A[i + offset] = a0 + a1 + a2 + a3; A[i + offset + p] = (a0 - a1 + a2na3iimag) * irot; A[i + offset + 2 * p] = (a0 + a1 - (a2 + a3)) * irot2; A[i + offset + 3 * p] = (a0 - a1 - a2na3iimag) * irot3; } irot *= irate3[__builtin_ctz(~s)]; } } if (le >= 1) { int p = 1 << (h - 1); for (int i = 0; i < p; i++) { auto ajp = A[i] - A[i + p]; A[i] += A[i + p]; A[i + p] = ajp; } } if (f) { mint inv = mint(1) / n; for (int i = 0; i < n; i++) { A[i] *= inv; } } } static std::vector multiply(std::vector A, std::vector B) { int need = int(A.size() + B.size()) - 1; if (std::min(A.size(), B.size()) < 60u) { std::vector C(need, 0); for (size_t i = 0; i < A.size(); i++) for (size_t j = 0; j < B.size(); j++) { C[i + j] += A[i] * B[j]; } return C; } int sz = 1; while (sz < need) sz <<= 1; A.resize(sz, 0); B.resize(sz, 0); ntt(A); ntt(B); mint inv = mint(1) / sz; for (int i = 0; i < sz; i++) A[i] *= B[i] * inv; intt(A, false); A.resize(need); return A; } }; template std::vector NumberTheoreticTransform::roots = std::vector(); template std::vector NumberTheoreticTransform::iroots = std::vector(); template std::vector NumberTheoreticTransform::rate3 = std::vector(); template std::vector NumberTheoreticTransform::irate3 = std::vector(); template int NumberTheoreticTransform::max_base = 0; template T modpow(T a, S b, T MOD) { T ret = 1; while (b > 0) { if (b & 1) { ret *= a; ret %= MOD; } a *= a; a %= MOD; b >>= 1; } return ret; } template T cipolla(T a, T MOD) { if (MOD == 2) return a; else if (a == 0) return 0; else if (modpow(a, (MOD - 1) / 2, MOD) != 1) return -1; T b = 1; while (modpow((b * b + MOD - a) % MOD, (MOD - 1) / 2, MOD) == 1) { b++; } T base = (b * b + MOD - a) % MOD; auto multi = [&](T a0, T b0, T a1, T b1) -> std::pair { return {(a0 * a1 + (b0 * b1 % MOD) * base) % MOD, (a0 * b1 + b0 * a1) % MOD}; }; auto pow_ = [&](auto self, T a, T b, T n) -> std::pair { if (n == 0) return {1, 0}; auto tmp = multi(a, b, a, b); auto ret = self(self, tmp.first, tmp.second, n / 2); if (n & 1) { ret = multi(ret.first, ret.second, a, b); } return ret; }; return pow_(pow_, b, 1LL, (MOD + 1) / 2).first; } template struct FormalPowerSeries : std::vector { using std::vector::vector; using FPS = FormalPowerSeries; static std::vector inv_x; void shrink() { while (this->size() && this->back() == mint(0)) { this->pop_back(); } } FPS &operator+=(const FPS &A) { if (A.size() > this->size()) this->resize(A.size()); for (size_t i = 0; i < A.size(); i++) (*this)[i] += A[i]; return *this; } FPS &operator+=(const mint &x) { if (this->empty()) this->resize(1); (*this)[0] += x; return *this; } FPS &operator-=(const FPS &A) { if (A.size() > this->size()) this->resize(A.size()); for (size_t i = 0; i < A.size(); i++) (*this)[i] -= A[i]; return *this; } FPS &operator-=(const mint &x) { if (this->empty()) this->resize(1); (*this)[0] -= x; return *this; } FPS &operator*=(const FPS &A) { if (this->empty() || A.empty()) { this->clear(); return *this; } auto res = NumberTheoreticTransform::multiply(*this, A); return *this = {begin(res), end(res)}; } FPS &operator*=(const mint &x) { for (size_t i = 0; i < this->size(); i++) (*this)[i] *= x; return *this; } FPS operator+(const FPS &A) const { return FPS(*this) += A; } FPS operator+(const mint &x) const { return FPS(*this) += x; } FPS operator-(const FPS &A) const { return FPS(*this) -= A; } FPS operator-(const mint &x) const { return FPS(*this) -= x; } FPS operator*(const FPS &A) const { return FPS(*this) *= A; } FPS operator*(const mint &x) const { return FPS(*this) *= x; } FPS operator-() const { FPS ret(this->size); for (size_t i = 0; i < this->size(); i++) ret[i] = -(*this)[i]; return ret; } FPS inv(int deg = -1) { assert((*this)[0] != mint(0)); if (deg == -1) deg = this->size(); FPS g = {mint(1) / (*this)[0]}; int l = 1; while (l < deg) { FPS tmp = g * 2; l <<= 1; FPS tmp2; g *= g; if (int(this->size()) >= l) tmp2 = FPS({this->begin(), this->begin() + l}) * g; else tmp2 = (*this) * g; g = tmp - tmp2; g.resize(l); } g.resize(deg); return g; } void iinv(int deg = -1) { *this = inv(deg); } FPS differential() { FPS ret(this->size() - 1); for (size_t i = 0; i < this->size() - 1; i++) ret[i] = (*this)[i + 1] * (i + 1); return ret; } void idifferential() { *this = this->differential(); } void extend_inv(int n) { int bn = inv_x.size(); if (n >= bn) { inv_x.resize(n + 1, 0); if (bn == 0) { inv_x[0] = 0; inv_x[1] = 1; bn = 2; } long long mod = mint::get_mod(); for (int i = bn; i <= n; i++) { inv_x[i] = mod - inv_x[mod % i].x * (mod / i) % mod; } } } FPS integral() { extend_inv(this->size()); FPS ret(this->size() + 1); for (size_t i = 0; i < this->size(); i++) ret[i + 1] = (*this)[i] * inv_x[i + 1]; return ret; } void iintegral() { *this = this->integral(); } FPS log(int deg = -1) { assert((*this)[0] == mint(1)); if (deg == -1) deg = this->size(); FPS B = (this->differential()) * (this->inv()); B.resize(deg - 1); return B.integral(); } void ilog(int deg = -1) { *this = this->log(deg); } FPS exp(int deg = -1) { assert((*this)[0] == mint(0)); if (deg == -1) deg = this->size(); FPS g = {1}; int l = 1; while (l < deg * 2) { l *= 2; FPS tmp = {1}; tmp -= g.log(l); if (int(this->size()) >= l) tmp += FPS({this->begin(), this->begin() + l}); else tmp += (*this); g *= tmp; g.resize(l); } g.resize(deg); return g; } void iexp(int deg = -1) { *this = this->exp(deg); } FPS pow(long long k, int deg = -1) { if (deg == -1) deg = this->size(); if (k == 0) { FPS ret(deg, 0); ret[0] = 1; return ret; } int p = -1; for (int i = 0; i < deg; i++) { if ((*this)[i] != 0) { p = i; break; } } if (p == -1 || p > deg / k) { FPS ret(deg, 0); return ret; } mint inv = mint(1) / (*this)[p]; FPS A = FPS({(*this).begin() + p, (*this).end()}); A *= inv; A.ilog(deg); A *= k % mint::get_mod(); A.iexp(deg); FPS B(p * k, 0); B.insert(B.end(), A.begin(), A.begin() + (deg - p * k)); B *= (*this)[p].pow(k); return B; } void ipow(long long k, int deg = -1) { *this = this->pow(k, deg); } FPS sqrt(int deg = -1) { if (deg == -1) deg = this->size(); if (this->size() == 0u) { FPS ret(deg, 0); return ret; } if ((*this)[0] == mint(0)) { for (size_t i = 1; i < this->size(); i++) { if ((*this)[i] != 0) { if (i & 1) { FPS ret; return ret; } if (deg <= int(i / 2)) break; FPS ret = FPS({this->begin() + i, this->end()}).sqrt(deg - i / 2); if (ret.size() == 0u) return ret; FPS ret2(i / 2, 0); ret2.insert(ret2.end(), ret.begin(), ret.end()); std::swap(ret, ret2); if (int(ret.size()) < deg) ret.resize(deg); return ret; } } FPS ret(deg, 0); return ret; } long long sq = cipolla((*this)[0].x, mint::get_mod()); if (sq == -1) { FPS ret; return ret; } mint inv2 = mint(1) / 2; FPS g = {sq}; int l = 1; while (l < deg) { l *= 2; if (int(this->size()) >= l) g += FPS({this->begin(), this->begin() + l}) * g.inv(l); else g += (*this) * g.inv(l); g *= inv2; } g.resize(deg); return g; } void isqrt(int deg = -1) { *this = this->sqrt(deg); } FPS taylorshift(mint a) { auto A = (*this); int deg = A.size(); extend_inv(deg); mint fac = 1; for (int i = 0; i < deg; i++) { A[i] *= fac; fac *= (i + 1); } reverse(A.begin(), A.end()); FPS g(deg, 0); g[0] = 1; for (int i = 1; i < deg; i++) g[i] = g[i - 1] * a * inv_x[i]; A *= g; if (int(A.size()) > deg) A.resize(deg); reverse(A.begin(), A.end()); mint invfac = 1; for (int i = 0; i < deg; i++) { A[i] *= invfac; invfac *= inv_x[i + 1]; } return A; } void itaylorshift(mint a) { int deg = this->size(); extend_inv(deg); mint fac = 1; for (int i = 0; i < deg; i++) { (*this)[i] *= fac; fac *= (i + 1); } reverse(this->begin(), this->end()); FPS g(deg, 0); g[0] = 1; for (int i = 1; i < deg; i++) g[i] = g[i - 1] * a * inv_x[i]; (*this) *= g; if (int(this->size()) > deg) this->resize(deg); reverse(this->begin(), this->end()); mint invfac = 1; for (int i = 0; i < deg; i++) { (*this)[i] *= invfac; invfac *= inv_x[i + 1]; } } std::pair division_of_polynomial(FPS G) { FPS F = *this; if (F.size() < G.size()) { return {{}, F}; } reverse(F.begin(), F.end()); reverse(G.begin(), G.end()); int deg = F.size() - G.size() + 1u; auto Q = F * G.inv(deg); if (int(Q.size()) > deg) Q.resize(deg); reverse(Q.begin(), Q.end()); reverse(F.begin(), F.end()); reverse(G.begin(), G.end()); auto R = F - G * Q; R.shrink(); return {Q, R}; } std::vector multipoint_evaluation(std::vector &X) { int m = X.size(); int m2 = 1; while (m2 <= m - 1) m2 *= 2; std::vector G(m2 << 1, FPS(1, 1)); for (int i = 0; i < m; i++) G[m2 + i] = {-X[i], 1}; for (int i = m2 - 1; i >= 0; i--) G[i] = G[i << 1] * G[(i << 1) | 1]; G[1] = this->division_of_polynomial(G[1]).second; for (int i = 2; i < m2 + m; i++) G[i] = G[i >> 1].division_of_polynomial(G[i]).second; std::vector Y(m); for (int i = 0; i < m; i++) { if (G[m2 + i].empty()) Y[i] = 0; else Y[i] = G[m2 + i][0]; } return Y; } std::vector multipoint_evaluation(std::vector &X) { int m = X.size(); int m2 = 1; while (m2 <= m - 1) m2 *= 2; std::vector G(m2 << 1, FPS(1, 1)); for (int i = 0; i < m; i++) G[m2 + i] = {-X[i], 1}; for (int i = m2 - 1; i >= 0; i--) G[i] = G[i << 1] * G[(i << 1) | 1]; G[1] = this->division_of_polynomial(G[1]).second; for (int i = 2; i < m2 + m; i++) G[i] = G[i >> 1].division_of_polynomial(G[i]).second; std::vector Y(m); for (int i = 0; i < m; i++) { if (G[m2 + i].empty()) Y[i] = 0; else Y[i] = G[m2 + i][0].x; } return Y; } friend std::ostream &operator<<(std::ostream &os, const FPS &A) { for (size_t i = 0; i < A.size(); i++) { os << A[i]; if (i != A.size() - 1) os << ' '; } return os; } friend std::istream &operator>>(std::istream &is, FPS &A) { for (size_t i = 0; i < A.size(); i++) { is >> A[i]; } return (is); } }; template std::vector FormalPowerSeries::inv_x = std::vector(); using FPS = FormalPowerSeries; void solve() { INT(n); VEC(int, A, n); int m = A[0]; Combination Comb(m + 10); vec(ll, C, m + 1, 0); for (auto a : A) C[a]++; vec(mint, pow2, m + 1, 1); fori(i, 1, m + 1) pow2[i] = pow2[i - 1] * 2; FPS F(m + 1); FPS G(m + 1); fori(i, m + 1) G[i] = Comb.invfact[i]; auto dc = [&](auto &&self, int l, int r) -> void { if (l + 1 == r) { if (l == 0) { F[0] = 1; } else { F[l] *= pow2[C[l]] - 1; } return; } int m = (l + r) / 2; self(self, l, m); auto L = FPS(F.begin() + l, F.begin() + m); auto R = FPS(G.begin(), G.begin() + r - l); auto res = L * R; fori(i, m, r) { F[i] += res[i - l]; } self(self, m, r); }; dc(dc, 0, m + 1); mint ans = 0; fori(i, 1, m + 1) { ans += F[i] * Comb.fact[i]; } print(ans); } int main() { #ifndef INTERACTIVE cin.tie(0)->sync_with_stdio(0); #endif // cout << fixed << setprecision(12); int t; t = 1; // cin >> t; while (t--) solve(); return 0; }