#include using u32 = unsigned; using i64 = long long; using u64 = unsigned long long; using i128 = __int128; using u128 = unsigned __int128; template struct static_modint { private: static_assert((P >> (sizeof(U0) * 8 - 1)) == 0, "'Mod' must less than max(U0)/2"); static constexpr U0 Mod = P; U0 x; template static constexpr U0 safeMod(T x) { if constexpr (std::is_unsigned::value) { return static_cast(x % Mod); } else { ((x %= static_cast(Mod)) < 0) && (x += Mod); return static_cast(x); } } public: constexpr static_modint(): x(static_cast(0)) {} template constexpr static_modint(T _x): x(safeMod(_x)) {} static constexpr static_modint from_raw(U0 _x) noexcept { static_modint x; return x.x = _x, x; } static constexpr U0 getMod() { return Mod; } template explicit constexpr operator T() const { return static_cast(x); } constexpr static_modint &operator += (const static_modint &rhs) { x += rhs.x, (x - Mod) >> (sizeof(U0) * 8 - 1) || (x -= Mod); return *this; } constexpr static_modint &operator -= (const static_modint &rhs) { x -= rhs.x, x >> (sizeof(U0) * 8 - 1) && (x += Mod); return *this; } constexpr static_modint &operator *= (const static_modint &rhs) { x = static_cast(static_cast(x) * static_cast(rhs.x) % static_cast(Mod)); return *this; } constexpr static_modint &operator /= (const static_modint &rhs) { return (*this *= rhs.inv()); } friend constexpr static_modint fma(const static_modint &a, const static_modint &b, const static_modint &c) { return from_raw((static_cast(a.x) * b.x + c.x) % Mod); } friend constexpr static_modint fam(const static_modint &a, const static_modint &b, const static_modint &c) { return from_raw((a.x + static_cast(b.x) * c.x) % Mod); } friend constexpr static_modint fms(const static_modint &a, const static_modint &b, const static_modint &c) { return from_raw((static_cast(a.x) * b.x + Mod - c.x) % Mod); } friend constexpr static_modint fsm(const static_modint &a, const static_modint &b, const static_modint &c) { return from_raw((a.x + static_cast(Mod - b.x) * c.x) % Mod); } constexpr static_modint inv() const { U0 a = Mod, b = x; S0 y = 0, z = 1; while (b) { const U0 q = a / b; const U0 c = a - q * b; a = b, b = c; const S0 w = y - static_cast(q) * z; y = z, z = w; } return from_raw(y < 0 ? y + Mod : y); } friend constexpr static_modint inv(const static_modint &x) { return x.inv(); } friend constexpr static_modint operator + (const static_modint &x) { return x; } friend constexpr static_modint operator - (static_modint x) { x.x = x.x ? (Mod - x.x) : 0U; return x; } constexpr static_modint &operator ++ () { return *this += 1; } constexpr static_modint &operator -- () { return *this -= 1; } constexpr static_modint operator ++ (int) { static_modint v = *this; return *this += 1, v; } constexpr static_modint operator -- (int) { static_modint v = *this; return *this -= 1, v; } friend constexpr static_modint operator + (static_modint x, const static_modint &y) { return x += y; } friend constexpr static_modint operator - (static_modint x, const static_modint &y) { return x -= y; } friend constexpr static_modint operator * (static_modint x, const static_modint &y) { return x *= y; } friend constexpr static_modint operator / (static_modint x, const static_modint &y) { return x /= y; } template constexpr static_modint pow(T y) const { if (y < 0) return inv().pow(- y); static_modint x = *this, ans = from_raw(1U); for (; y; y >>= 1, x *= x) { if (y & 1) { ans *= x; } } return ans; } template friend constexpr static_modint pow(const static_modint &x, T y) { return x.pow(y); } std::pair sqrt() const { if (x == 0U) { return {true, from_raw(0)}; } if (Mod == 2U) { return {true, from_raw(1)}; } if (pow((Mod - 1) / 2) != from_raw(1)) { return {false, from_raw(0)}; } static std::mt19937_64 rnd(std::chrono::system_clock::now().time_since_epoch().count()); std::uniform_int_distribution uid(1U, Mod - 1); static_modint x, y; do { x = from_raw(uid(rnd)); y = x * x - *this; } while (y.pow((Mod - 1) / 2) == from_raw(1)); auto mul = [](std::pair &f, const std::pair &g, const static_modint &h) { f = {f.first * g.first + f.second * g.second * h, f.first * g.second + f.second * g.first}; }; std::pair f{x, 1}, g{1, 0}; auto exp = (Mod + 1) / 2; for (; exp; exp >>= 1, mul(f, f, y)) { if (exp & 1) { mul(g, f, y); } } return {true, from_raw(std::min(g.first.x, Mod - g.first.x))}; } friend std::pair sqrt(const static_modint &x) { return x.sqrt(); } friend constexpr std::istream& operator >> (std::istream& is, static_modint &x) { S0 y; is >> y, x = y; return is; } friend constexpr std::ostream& operator << (std::ostream& os, const static_modint &x) { return os << x.x; } friend constexpr bool operator == (const static_modint &x, const static_modint &y) { return x.x == y.x; } friend constexpr bool operator != (const static_modint &x, const static_modint &y) { return x.x != y.x; } friend constexpr bool operator <= (const static_modint &x, const static_modint &y) { return x.x <= y.x; } friend constexpr bool operator >= (const static_modint &x, const static_modint &y) { return x.x >= y.x; } friend constexpr bool operator < (const static_modint &x, const static_modint &y) { return x.x < y.x; } friend constexpr bool operator > (const static_modint &x, const static_modint &y) { return x.x > y.x; } }; template using sm32 = static_modint; template using sm64 = static_modint; using Z = sm32<998244353U>; struct Combination { int N; std::vector _fac, _ifac, _inv; Combination(int n): N(0), _fac{1}, _ifac{1}, _inv{0} { init(n); } Combination(): N(0), _fac{1}, _ifac{1}, _inv{0} {} void init(int n) { if (n <= N) return; _fac.resize(n + 1), _ifac.resize(n + 1), _inv.resize(n + 1); for (int i = N + 1; i <= n; i ++) _fac[i] = _fac[i - 1] * i; _ifac[n] = _fac[n].inv(); for (int i = n; i > N; i --) _ifac[i - 1] = _ifac[i] * i, _inv[i] = _ifac[i] * _fac[i - 1]; N = n; return; } Z fac(int n) { init(n << 1); return n < 0 ? Z::from_raw(0) : _fac[n]; } Z ifac(int n) { init(n << 1); return n < 0 ? Z::from_raw(0) : _ifac[n]; } Z inv(int n) { init(n << 1); return n < 0 ? Z::from_raw(0) : _inv[n]; } Z binom(int n, int m) { if (n < m || n < 0 || m < 0) return 0; return fac(n) * ifac(m) * ifac(n - m); } } comb; template struct polynomial: public std::vector { //private: static constexpr Z g = 3; static constexpr auto Mod = Z::getMod(); static constexpr int log_ord = []() { auto x = Mod - 1; int y = 0; while (~x & 1) { x >>= 1, ++ y; } return y; }(); static constexpr std::array invn = []() { std::array inv{}; for (int i = 0; i <= log_ord; i ++) { inv[i] = Z(1 << i).inv(); } return inv; }(); static std::pair get_root(const int &n) { static std::vector root{Z::from_raw(1)}; static std::vector inv_root{Z::from_raw(1)}; if (static_cast(root.size()) < n) { int i = root.size(); root.resize(n), inv_root.resize(n); for (; i != n; i <<= 1) { const Z w = g.pow(Mod / (i << 2)), iw = w.inv(); for (int j = 0; j != i; j ++) { root[i + j] = root[j] * w; inv_root[i + j] = inv_root[j] * iw; } } } return {root.data(), inv_root.data()}; } static constexpr int get_len(int n) { return n < 3 ? n : 2 << std::__lg(n - 1); } static void dif_n(Z *f, const int &n) { const Z* rt = get_root(n).first; for (int i = n; i >>= 1; ) { for (int j = 0, k = 0; j != n; j += i << 1, ++ k) { for (int p = j, q = j + i; p != j + i; ++ p, ++ q) { const Z u = f[p], v = f[q] * rt[k]; f[p] = u + v, f[q] = u - v; } } } } static void dit_n(Z *f, const int &n) { const Z* irt = get_root(n).second; for (int i = 1; i != n; i <<= 1) { for (int j = 0, k = 0; j != n; j += i << 1, ++ k) { for (int p = j, q = j + i; p != j + i; ++ p, ++ q) { const Z u = f[p], v = f[q]; f[p] = u + v, f[q] = (u - v) * irt[k]; } } } const Z inv = invn[std::__lg(n)]; for (int i = 0; i < n; i ++) { f[i] *= inv; } } static void dif_rhalf_n(Z *f, const int &n) { const Z* rt = get_root(n).first; for (int i = n, m = 1; i >>= 1; m <<= 1) { for (int j = 0, k = 0; j != n; j += i << 1, ++ k) { for (int p = j, q = j + i; p != j + i; ++ p, ++ q) { const Z u = f[p], v = f[q] * rt[m + k]; f[p] = u + v, f[q] = u - v; } } } } static void dit_rhalf_n(Z *f, const int &n) { const Z* irt = get_root(n).second; for (int i = 1, m = n; m >>= 1; i <<= 1) { for (int j = 0, k = 0; j != n; j += i << 1, ++ k) { for (int p = j, q = j + i; p != j + i; ++ p, ++ q) { const Z u = f[p], v = f[q]; f[p] = u + v, f[q] = (u - v) * irt[m + k]; } } } const Z inv = invn[std::__lg(n)]; for (int i = 0; i < n; i ++) { f[i] *= inv; } } static void neg_n(const Z *a, const int &n, Z *b) { for (int i = 0; i < n; i ++) { b[i] = -a[i]; } } static void add_n(const Z *a, const Z *b, const int &n, Z *c) { for (int i = 0; i < n; i ++) { c[i] = a[i] + b[i]; } } static void sub_n(const Z *a, const Z *b, const int &n, Z *c) { for (int i = 0; i < n; i ++) { c[i] = a[i] - b[i]; } } static void dot_n(const Z *a, const Z *b, const int &n, Z *c) { for (int i = 0; i < n; i ++) { c[i] = a[i] * b[i]; } } static void mul_c_n(const Z *a, const Z &c, const int &n, Z *b) { for (int i = 0; i < n; i ++) { b[i] = a[i] * c; } } static void deriv_n(const Z *a, const int &n, Z *b) { for (int i = 1; i != n; i ++) { b[i - 1] = a[i] * i; } b[n - 1] = Z::from_raw(0); } static void integ_n(const Z *a, const int &n, Z *b) { comb.init(n); for (int i = n - 1; i; i --) { b[i] = a[i - 1] * comb.inv(i); } b[0] = Z::from_raw(0); } static void zero_n(Z *a, const int &n) { for (int i = 0; i < n; i ++) { a[i] = Z::from_raw(0); } } static void copy_n(const Z *a, const int &n, Z *b) { memcpy(b, a, sizeof(Z) * n); } static polynomial convolution_fft(polynomial f, polynomial g) { const int n = f.size() + g.size() - 1, m = get_len(n); f.resize(m), dif_n(f.data(), m); g.resize(m), dif_n(g.data(), m); dot_n(f.data(), g.data(), m, f.data()); dit_n(f.data(), m), f.resize(n); return f; } static polynomial convolution_naive(const polynomial &f, const polynomial &g) { const int n = f.size(), m = g.size(); if (__builtin_expect(f.empty() || g.empty(), 0)) { return polynomial{}; } polynomial fg(n + m - 1); for (int i = 0; i != n; i ++) { for (int j = 0; j != m; j ++) { fg[i + j] += f[i] * g[j]; } } return fg; } polynomial pow_one(const int &n, const int &k) const { if (__builtin_expect(k == 0, 0)) { polynomial f(this->size()); f[0] = 1; return f; } return (k == 1) ? mod_xk(n) : (mod_xk(n).log() * k).exp(); } polynomial pow_ord_zero(const int &n, int k_mod_p, int k_mod_phi_p) const { if ((*this)[0] == Z::from_raw(1)) { return pow_one(n, k_mod_p); } return (*this * (*this)[0].inv()).pow_one(n, k_mod_p) * (*this)[0].pow(k_mod_phi_p); } polynomial sqrt_ord_zero(const Z &s) const { constexpr Z c = -Z(2).inv(); const int shrink_len = this->size(); const int n = get_len(shrink_len); polynomial res(shrink_len), inv_res(shrink_len); polynomial dft_res(n), dft_inv_res(n), f(n); int N = 1, N2 = 2; res[0] = dft_res[0] = s, inv_res[0] = s.inv(); while (N < shrink_len) { const int newN = (N2 == n) ? shrink_len : N2; dot_n(dft_res.data(), dft_res.data(), N, dft_res.data()), dit_n(dft_res.data(), N); sub_n(dft_res.data(), this->data(), N, dft_res.data() + N); sub_n(dft_res.data() + N, this->data() + N, newN - N, dft_res.data() + N); zero_n(dft_res.data(), N), dif_n(dft_res.data(), N2); copy_n(inv_res.data(), N, dft_inv_res.data()), dif_n(dft_inv_res.data(), N2); dot_n(dft_res.data(), dft_inv_res.data(), N2, dft_res.data()), dit_n(dft_res.data(), N2); mul_c_n(dft_res.data() + N, c, newN - N, res.data() + N); if (__builtin_expect(N2 < n, 1)) { copy_n(res.data(), N2, f.data()), dif_n(f.data(), N2); copy_n(f.data(), N2, dft_res.data()); dot_n(f.data(), dft_inv_res.data(), N2, f.data()), dit_n(f.data(), N2); zero_n(f.data(), N), dif_n(f.data(), N2); dot_n(f.data(), dft_inv_res.data(), N2, f.data()), dit_n(f.data(), N2); neg_n(f.data() + N, N, inv_res.data() + N); } N <<= 1, N2 <<= 1; } return res; } public: polynomial(): std::vector() {} explicit polynomial(const int &n): std::vector(n) {} explicit polynomial(const std::vector &a): std::vector(a) {} explicit polynomial(const std::initializer_list &a): std::vector(a) {} template> explicit polynomial(_InputIterator __first, _InputIterator __last): std::vector(__first, __last) {} template explicit polynomial(const int &n, const F &f): std::vector(n) { for (int i = 0; i != n; i ++) { (*this)[i] = f(i); } } int ord() const { int ord = 0; while (ord != static_cast(this->size()) && !(*this)[ord]) { ++ ord; } return ord; } int deg() const { int deg = static_cast(this->size()) - 1; while (deg >= 0 && !(*this)[deg]) { -- deg; } return deg; } void remove0() { while (this->size() && !this->back()) { this->pop_back(); } } polynomial mul_xk(const int &k) const { auto f = *this; f.insert(f.begin(), k, Z::from_raw(0)); return f; } polynomial div_xk(const int &k) const { return polynomial(this->begin() + std::min(static_cast(this->size()), k), this->end()); } polynomial mod_xk(const int &k) const { return polynomial(this->begin(), this->begin() + std::min(static_cast(this->size()), k)); } polynomial compose_cx(const Z &c) const { const int n = this->size(); auto f = *this; Z ci = Z::from_raw(1); for (int i = 1; i < n; i ++) { f[i] *= (ci *= c); } return f; } Z operator () (const Z &x) const { const int n = this->size(); Z y = Z::from_raw(0); for (int i = n - 1; i >= 0; i --) { (y *= x) += (*this)[i]; } return y; } polynomial operator + () const { return *this; } polynomial operator - () const { auto f = *this; neg_n(f.data(), f.size(), f.data()); return f; } polynomial& operator += (const Z &rhs) { if (__builtin_expect(this->empty(), 0)) { return polynomial(1, rhs); } return (*this)[0] += rhs, *this; } polynomial& operator -= (const Z &rhs) { if (__builtin_expect(this->empty(), 0)) { return polynomial(1, -rhs); } return (*this)[0] -= rhs, *this; } polynomial& operator *= (const Z &rhs) { mul_c_n(this->data(), rhs, this->size(), this->data()); return *this; } polynomial& operator /= (const Z &rhs) { mul_c_n(this->data(), rhs.inv(), this->size(), this->data()); return *this; } polynomial& operator += (const polynomial &rhs) { if (this->size() < rhs.size()) { this->resize(rhs.size()); } add_n(this->data(), rhs.data(), rhs.size(), this->data()); return *this; } polynomial& operator -= (const polynomial &rhs) { if (this->size() < rhs.size()) { this->resize(rhs.size()); } sub_n(this->data(), rhs.data(), rhs.size(), this->data()); return *this; } polynomial& operator *= (const polynomial &rhs) { return *this = *this * rhs; } friend polynomial operator + (polynomial lhs, const Z &rhs) { return lhs += rhs; } friend polynomial operator - (polynomial lhs, const Z &rhs) { return lhs -= rhs; } friend polynomial operator * (polynomial lhs, const Z &rhs) { return lhs *= rhs; } friend polynomial operator / (polynomial lhs, const Z &rhs) { return lhs /= rhs; } friend polynomial operator + (polynomial lhs, const polynomial &rhs) { return lhs += rhs; } friend polynomial operator - (polynomial lhs, const polynomial &rhs) { return lhs -= rhs; } friend polynomial operator * (const polynomial &f, const polynomial &g) { if (std::min(f.size(), g.size()) > 8 && std::max(f.size(), g.size()) > 128) { return convolution_fft(f, g); } else { return convolution_naive(f, g); } } polynomial deriv() const { if (__builtin_expect(this->emtpy())) { return *this; } auto f = *this; deriv_n(f.data(), f.size(), f.data()); return f.pop_back(), f; } polynomial integ() const { auto f = *this; f.resize(this->size() + 1); integ_n(f.data(), f.size(), f.data()); return f.pop_back(), f; } polynomial inv() const { if (__builtin_expect(this->empty(), 0)) { return polynomial{}; } assert((*this)[0] != Z::from_raw(0)); const int shrink_len = this->size(); const int n = get_len(shrink_len); polynomial res(shrink_len), f(n), g(n); res[0] = (*this)[0].inv(); int N = 1, N2 = 2; while (N < shrink_len) { const int newN = (N2 == n) ? shrink_len : N2; copy_n(this->data(), newN, f.data()), dif_n(f.data(), N2); copy_n(res.data(), N, g.data()), dif_n(g.data(), N2); dot_n(f.data(), g.data(), N2, f.data()), dit_n(f.data(), N2); zero_n(f.data(), N), dif_n(f.data(), N2); dot_n(f.data(), g.data(), N2, f.data()), dit_n(f.data(), N2); neg_n(f.data() + N, newN - N, res.data() + N); N = newN, N2 = N << 1; } return res; } polynomial div(const polynomial &g) const { if (g.size() > this->size()) { return polynomial{}; } assert(g.size()); const int n = this->size() - g.size() + 1; polynomial q(n), r(n); std::reverse_copy(this->begin() + g.size() - 1, this->end(), q.begin()); std::reverse_copy(g.begin() + std::max(0, static_cast(g.size()) - n), g.end(), r.begin()); q = (q * r.inv()).mod_xk(n), std::reverse(q.begin(), q.end()); return q; } std::pair div_mod(const polynomial &g) const { if (g.size() > this->size()) { return {polynomial{}, *this}; } assert(g.size()); const int n = g.size() - 1; auto q = div(g); return {q, mod_xk(n) - (g.mod_xk(n) * q.mod_xk(n)).mod_xk(n)}; } polynomial exp() const { if (__builtin_expect(this->empty(), 0)) { return polynomial{}; } assert((*this)[0] == Z::from_raw(0)); if (__builtin_expect(this->size() == 1, 0)) { return polynomial{Z::from_raw(1)}; } const int shrink_len = this->size(); const int n = get_len(shrink_len); int N = 1, N2 = 2; polynomial res(shrink_len), inv_res(shrink_len); polynomial dft_res(n), dft_inv_res(n), f(n); res[0] = inv_res[0] = dft_res[0] = dft_res[1] =1; while (N < shrink_len) { const int newN = (N2 == n) ? shrink_len : N2; deriv_n(this->data(), N, f.data()), dif_n(f.data(), N); dot_n(f.data(), dft_res.data(), N, f.data()), dit_n(f.data(), N); f[N - 1] = -f[N - 1], deriv_n(res.data(), N, f.data() + N); sub_n(f.data() + N, f.data(), N - 1, f.data() + N); zero_n(f.data(), N - 1), dif_n(f.data(), N2); copy_n(inv_res.data(), N, dft_inv_res.data()), dif_n(dft_inv_res.data(), N2); dot_n(f.data(), dft_inv_res.data(), N2, f.data()), dit_n(f.data(), N2); integ_n(f.data(), N2, f.data()), zero_n(f.data(), N); sub_n(f.data() + N, this->data() + N, newN - N, f.data() + N), dif_n(f.data(), N2); dot_n(f.data(), dft_res.data(), N2, f.data()), dit_n(f.data(), N2); neg_n(f.data() + N, newN - N, res.data() + N); if (__builtin_expect(N2 < n, 1)) { copy_n(res.data(), N2, dft_res.data()), dif_n(dft_res.data(), N2 + N2); copy_n(dft_res.data(), N2, f.data()); dot_n(f.data(), dft_inv_res.data(), N2, f.data()), dit_n(f.data(), N2); zero_n(f.data(), N), dif_n(f.data(), N2); dot_n(f.data(), dft_inv_res.data(), N2, f.data()), dit_n(f.data(), N2); neg_n(f.data() + N, N, inv_res.data() + N); } N <<= 1, N2 <<= 1; } return res; } polynomial log() const { if (__builtin_expect(this->empty(), 0)) { return polynomial{}; } assert((*this)[0] == Z::from_raw(1)); polynomial f(this->size()); deriv_n(this->data(), this->size(), f.data()); f *= inv(), f.resize(this->size()); integ_n(f.data(), this->size(), f.data()); return f; } polynomial pow(int k_chkmin_sz, int k_mod_p, int k_mod_phi_p) const { if (__builtin_expect(this->empty(), 0)) { return polynomial{}; } if (__builtin_expect(k_chkmin_sz == 0, 0)) { polynomial f(this->size()); return f[0] = Z::from_raw(1), f; } const int i = ord(); if (static_cast(i) * k_chkmin_sz >= static_cast(this->size())) { return polynomial(this->size()); } return div_xk(i).pow_ord_zero(this->size() - i * k_chkmin_sz, k_mod_p, k_mod_phi_p).mul_xk(i * k_chkmin_sz); } template inline polynomial pow(T x) const { if (x < 0) return inv().pow(-x); return pow(x < static_cast(this->size()) ? x : this->size(), x % Mod, x % (Mod - 1)); } std::pair sqrt() const { const int i = ord(); if (i == static_cast(this->size())) { return {true, polynomial(this->size())}; } if (i & 1) { return {false, polynomial{}}; } auto [o, s] = (*this)[i].sqrt(); if (o == false) { return {false, polynomial{}}; } auto x = div_xk(i); x.insert(x.end(), i >> 1, 0); return {true, x.sqrt_ord_zero(s).mul_xk(i >> 1)}; } polynomial composional_inv() const { assert(this->size() >= 2U); assert((*this)[0] == Z::from_raw(0)); assert((*this)[1] != Z::from_raw(0)); constexpr Z inv2 = Z(2).inv(); const int shrink_len = this->size(); const int n = get_len(shrink_len << 1); const Z *irt = get_root(2 * n).second; const Z invf1 = (*this)[1].inv(); polynomial P(n), Q(n); polynomial dftP(n << 1), dftQ(n << 1); P[0] = Z::from_raw(1); Z v = -1; for (int i = 1; i < shrink_len; i ++) { v *= invf1, Q[i] = (*this)[i] * v; } int N = shrink_len, M = 1, newN = 0; for (; N > 1; N = newN, M <<= 1) { newN = (N + 1) >> 1; const int len = get_len((N * M) << 2); zero_n(dftP.data(), len), zero_n(dftQ.data(), len); for (int j = 0; j < M; j ++) { copy_n(P.data() + j * N, N, dftP.data() + 2 * j * N); copy_n(Q.data() + j * N, N, dftQ.data() + 2 * j * N); } dftQ[2 * N * M] = 1; dif_n(dftP.data(), len); dif_n(dftQ.data(), len); P.resize(len >> 1), Q.resize(len >> 1); for (int i = 0; i < len; i += 2) { Q[i >> 1] = dftQ[i] * dftQ[i + 1]; } if (N & 1) { for (int i = 0; i < len; i += 2) { P[i >> 1] = (dftP[i] * dftQ[i + 1] + dftP[i + 1] * dftQ[i]) * inv2; } } else { for (int i = 0; i < len; i += 2) { P[i >> 1] = (dftP[i] * dftQ[i + 1] - dftP[i + 1] * dftQ[i]) * irt[i >> 1] * inv2; } } dit_n(P.data(), len >> 1); dit_n(Q.data(), len >> 1); if (N * M * 4 == len) { -- Q[0]; } for (int j = 1; j < M * 2; j ++) { copy_n(P.data() + j * N, newN, P.data() + j * newN); copy_n(Q.data() + j * N, newN, Q.data() + j * newN); } } P.resize(M * newN); std::reverse(P.begin(), P.end()); P.resize(shrink_len); for (int i = 1; i < shrink_len; i ++) { P[i] *= (shrink_len - 1) * comb.inv(i); } std::reverse(P.begin(), P.end()); P = P.pow_one(shrink_len - 1, (int)(-comb.inv(shrink_len - 1))) * invf1; P.insert(P.begin(), 0); return P; } template friend Z linear_recurrence(polynomial f, polynomial g, T k) { constexpr Z inv2 = Z(2).inv(); const int n = g.size(); const int N = get_len(2 * n); const Z* irt = get_root(N).second; polynomial P(N), Q(N), dft(N); copy_n(f.data(), f.size(), P.data()); copy_n(g.data(), n, Q.data()); dif_n(P.data(), N); dif_n(Q.data(), N); while (k) { if (k & 1) { for (int i = 0; i < N / 2; i ++) { dft[i] = (P[2 * i] * Q[2 * i + 1] - P[2 * i + 1] * Q[2 * i]) * irt[i] * inv2; } } else { for (int i = 0; i < N / 2; i ++) { dft[i] = (P[2 * i] * Q[2 * i + 1] + P[2 * i + 1] * Q[2 * i]) * inv2; } } copy_n(dft.data(), N / 2, P.data()); dit_n(dft.data(), N / 2); dif_rhalf_n(dft.data(), N / 2); copy_n(dft.data(), N / 2, P.data() + N / 2); for (int i = 0; i < N / 2; i ++) { dft[i] = Q[2 * i] * Q[2 * i + 1]; } copy_n(dft.data(), N / 2, Q.data()); dit_n(dft.data(), N / 2); dif_rhalf_n(dft.data(), N / 2); copy_n(dft.data(), N / 2, Q.data() + N / 2); k >>= 1; } dit_n(P.data(), N); dit_n(Q.data(), N); return P[0] / Q[0]; } }; using poly = polynomial; constexpr int K = 62; constexpr Z inv2 = Z(2).inv(); int main() { #ifdef LOCAL freopen("!in.in", "r", stdin); freopen("!out.out", "w", stdout); #endif std::ios::sync_with_stdio(false); std::cin.tie(nullptr); std::cout.tie(nullptr); i64 N; std::cin >> N; std::array op{}; op.fill(1); // for (int o = 0; o < 5; o ++) { // std::cin >> op[o]; // } std::array Q{}, Q2{}, Q3{}; Q.fill(poly{1}); auto get2 = [&](const poly &P, const poly& Q, int o) { static poly p, A[2], B[2]; A[0].clear(), A[1].clear(); for (int i = 0; i < (int)Q.size(); i ++) { A[i % 2].push_back(Q[i]); } B[0].clear(), B[1].clear(); for (int i = 0; i < (int)P.size(); i ++) { B[i % 2].push_back(P[i]); } if (o == 0) { const int n = (P.size() + Q.size()) / 2; const int N = poly::get_len(n); const Z* rt = poly::get_root(N).first; A[0].resize(N); A[1].resize(N); B[0].resize(N); B[1].resize(N); poly::dif_n(A[0].data(), N); poly::dif_n(A[1].data(), N); poly::dif_n(B[0].data(), N); poly::dif_n(B[1].data(), N); p.resize(N); for (int i = 0; i < N; i ++) { if (i % 2 == 0) { p[i] = A[0][i] * B[0][i] - A[1][i] * B[1][i] * rt[i / 2]; } else { p[i] = A[0][i] * B[0][i] + A[1][i] * B[1][i] * rt[i / 2]; } } poly::dit_n(p.data(), N); p.resize(n); } else { const int n = (P.size() + Q.size() - 1) / 2; const int N = poly::get_len(n); A[0].resize(N); A[1].resize(N); B[0].resize(N); B[1].resize(N); poly::dif_n(A[0].data(), N); poly::dif_n(A[1].data(), N); poly::dif_n(B[0].data(), N); poly::dif_n(B[1].data(), N); p.resize(N); for (int i = 0; i < N; i ++) { p[i] = A[0][i] * B[1][i] - A[1][i] * B[0][i]; } poly::dit_n(p.data(), N); p.resize(n); } return p; }; auto get3 = [&](const poly &P, const poly& Q, int o) { static poly p, A[3], B[3]; A[0].clear(), A[1].clear(), A[2].clear(); for (int i = 0; i < (int)Q.size(); i ++) { A[i % 3].push_back(Q[i]); } B[0].clear(), B[1].clear(), B[2].clear(); for (int i = 0; i < (int)P.size(); i ++) { B[i % 3].push_back(P[i]); } if (o == 0) { const int n = (P.size() + 2 * Q.size()) / 3; const int N = poly::get_len(n); const Z* rt = poly::get_root(N).first; A[0].resize(N); A[1].resize(N); A[2].resize(N); B[0].resize(N); B[1].resize(N); B[2].resize(N); poly::dif_n(A[0].data(), N); poly::dif_n(A[1].data(), N); poly::dif_n(A[2].data(), N); poly::dif_n(B[0].data(), N); poly::dif_n(B[1].data(), N); poly::dif_n(B[2].data(), N); p.resize(N); for (int i = 0; i < N; i ++) { p[i] = 0; if (i % 2 == 0) { p[i] += (A[0][i] * A[0][i] - A[1][i] * A[2][i] * rt[i / 2]) * B[0][i]; p[i] += (A[2][i] * A[2][i] * rt[i / 2] - A[0][i] * A[1][i]) * B[2][i] * rt[i / 2]; p[i] += (A[1][i] * A[1][i] - A[0][i] * A[2][i]) * B[1][i] * rt[i / 2]; } else { p[i] += (A[0][i] * A[0][i] + A[1][i] * A[2][i] * rt[i / 2]) * B[0][i]; p[i] += (A[2][i] * A[2][i] * rt[i / 2] + A[0][i] * A[1][i]) * B[2][i] * rt[i / 2]; p[i] += (A[0][i] * A[2][i] - A[1][i] * A[1][i]) * B[1][i] * rt[i / 2]; } } poly::dit_n(p.data(), N); p.resize(n); } else if (o == 1) { const int n = (P.size() + 2 * Q.size() - 1) / 3; const int N = poly::get_len(n); const Z* rt = poly::get_root(N).first; A[0].resize(N); A[1].resize(N); A[2].resize(N); B[0].resize(N); B[1].resize(N); B[2].resize(N); poly::dif_n(A[0].data(), N); poly::dif_n(A[1].data(), N); poly::dif_n(A[2].data(), N); poly::dif_n(B[0].data(), N); poly::dif_n(B[1].data(), N); poly::dif_n(B[2].data(), N); p.resize(N); for (int i = 0; i < N; i ++) { p[i] = 0; if (i % 2 == 0) { p[i] += (A[0][i] * A[0][i] - A[1][i] * A[2][i] * rt[i / 2]) * B[1][i]; p[i] += (A[2][i] * A[2][i] * rt[i / 2] - A[0][i] * A[1][i]) * B[0][i]; p[i] += (A[1][i] * A[1][i] - A[0][i] * A[2][i]) * B[2][i] * rt[i / 2]; } else { p[i] += (A[0][i] * A[0][i] + A[1][i] * A[2][i] * rt[i / 2]) * B[1][i]; p[i] -= (A[2][i] * A[2][i] * rt[i / 2] + A[0][i] * A[1][i]) * B[0][i]; p[i] += (A[0][i] * A[2][i] - A[1][i] * A[1][i]) * B[2][i] * rt[i / 2]; } } poly::dit_n(p.data(), N); p.resize(n); } else { const int n = (P.size() + 2 * Q.size() - 2) / 3; const int N = poly::get_len(n); const Z* rt = poly::get_root(N).first; A[0].resize(N); A[1].resize(N); A[2].resize(N); B[0].resize(N); B[1].resize(N); B[2].resize(N); poly::dif_n(A[0].data(), N); poly::dif_n(A[1].data(), N); poly::dif_n(A[2].data(), N); poly::dif_n(B[0].data(), N); poly::dif_n(B[1].data(), N); poly::dif_n(B[2].data(), N); p.resize(N); for (int i = 0; i < N; i ++) { p[i] = 0; if (i % 2 == 0) { p[i] += (A[0][i] * A[0][i] - A[1][i] * A[2][i] * rt[i / 2]) * B[2][i]; p[i] += (A[2][i] * A[2][i] * rt[i / 2] - A[0][i] * A[1][i]) * B[1][i]; p[i] += (A[1][i] * A[1][i] - A[0][i] * A[2][i]) * B[0][i]; } else { p[i] += (A[0][i] * A[0][i] + A[1][i] * A[2][i] * rt[i / 2]) * B[2][i]; p[i] -= (A[2][i] * A[2][i] * rt[i / 2] + A[0][i] * A[1][i]) * B[1][i]; p[i] += (A[1][i] * A[1][i] - A[0][i] * A[2][i]) * B[0][i]; } } poly::dit_n(p.data(), N); p.resize(n); } return p; }; auto getQ = [&](auto&& self, int i, int j, const poly& P) -> void { if (i + j == K) { return; } Q[i + j] *= P; if (i == 0) { Q2[i + j] = P; } if (j == 0) { Q3[i + j] = P; } if (!j) { self(self, i + 1, j, get2(P, P, 0)); } self(self, i, j + 1, get3(P, P, 0)); }; getQ(getQ, 0, 0, poly{1, -op[0], -op[1], -op[2]}); Q2[0] = poly{1}; Q3[0] = poly{1}; for (int i = 1; i < K; i ++) { Q[i] *= Q[i - 1]; Q2[i] *= Q2[i - 1]; Q3[i] *= Q3[i - 1]; } auto trans2 = [&](int k, i64 N, const poly& P) -> std::pair { return {N / 2, get2(P, Q[k], N % 2) * Q2[k + 1]}; }; auto trans2C = [&](int k, i64 N, const poly& P) -> std::pair { return {N / 2, get2(P, Q[k], N % 2) * Q2[k + 1] * Q[0]}; }; auto trans3 = [&](int k, i64 N, const poly& P) -> std::pair { return {N / 3, get3(P, Q[k], N % 3) * Q3[k + 1]}; }; auto solve = [&](auto&& self, int k, const std::vector>& Fs, const std::vector>& Gs) -> Z { if (Fs.empty() && Gs.empty()) { return 0; } assert(k < K); std::vector> fs, gs; Z ans = 0; for (auto [N, P] : Fs) { if (N == 0) { continue; } gs.push_back(trans2C(k, N, P.mul_xk(1))); if (op[3]) { fs.push_back(trans2(k, N, P)); } if (op[4]) { fs.push_back(trans3(k, N, P)); } } for (auto [N, P] : Gs) { if (N == 0) { ans += P[0] / Q[k][0]; continue; } gs.push_back(trans2C(k, N, P)); } auto reduce = [&](std::vector>& Ps) { std::sort(Ps.begin(), Ps.end(), [&](const std::pair& x, const std::pair& y) { return x.first < y.first; }); std::vector> ps; for (int i = 0; i < (int)Ps.size(); i ++) { if (ps.empty() || ps.back().first != Ps[i].first) { ps.push_back(Ps[i]); } else { ps.back().second += Ps[i].second; } } Ps = std::move(ps); }; reduce(fs); reduce(gs); ans += self(self, k + 1, fs, gs); return ans; }; std::cout << solve(solve, 0, {{N, poly{1}}}, {}) << '\n'; return 0; }