#1143431 (C++23) No.3399 One Two Three Two Three

• 問題: No.3399 One Two Three Two Three
• ユーザー: definieren
• 提出日時: 2026-01-11 11:03:25
• 言語: C++23 (gcc 15.2.0 + boost 1.89.0)
• 結果: AC
• 実行時間: 1,961 ms / 4,000 ms
• コード長: 32,288 bytes
• コンパイル時間: 5,402 ms
• コンパイル使用メモリ: 368,344 KB
• 実行使用メモリ: 31,744 KB
• 最終ジャッジ日時: 2026-01-11 11:03:51
• 合計ジャッジ時間: 24,802 ms
• ジャッジサーバーID （参考情報）: judge4 / judge2

ファイルパターン	結果
sample	AC * 5
other	AC * 40

ソースコード

#include <bits/stdc++.h>

using u32 = unsigned;
using i64 = long long;
using u64 = unsigned long long;
using i128 = __int128;
using u128 = unsigned __int128;

template<class U0, class U1, class S0, U0 P>
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<class T>
	static constexpr U0 safeMod(T x) {
		if constexpr (std::is_unsigned<T>::value) {
			return static_cast<U0>(x % Mod);
		} else {
			((x %= static_cast<S0>(Mod)) < 0) && (x += Mod);
			return static_cast<U0>(x);
		}
	}
public:
	constexpr static_modint(): x(static_cast<U0>(0)) {}
	template<class T>
	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<class T>
	explicit constexpr operator T() const {
		return static_cast<T>(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<U0>(static_cast<U1>(x) * static_cast<U1>(rhs.x) % static_cast<U1>(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<U1>(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<U1>(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<U1>(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<U1>(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<S0>(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<class T>
	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<class T>
	friend constexpr static_modint pow(const static_modint &x, T y) {
		return x.pow(y);
	}
	std::pair<bool, static_modint> 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<U0> 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<static_modint, static_modint> &f,
			const std::pair<static_modint, static_modint> &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<static_modint, static_modint> 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<bool, static_modint> 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<u32 P>
using sm32 = static_modint<u32, u64, int, P>;
template<u64 P>
using sm64 = static_modint<u64, u128, i64, P>;
using Z = sm32<998244353U>;

struct Combination {
	int N;
	std::vector<Z> _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<class Z> 
struct polynomial: public std::vector<Z> {
//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<Z, log_ord + 1> invn = []() {
		std::array<Z, log_ord + 1> inv{};
		for (int i = 0; i <= log_ord; i ++) {
			inv[i] = Z(1 << i).inv();
		}
		return inv;
	}();
	static std::pair<Z*, Z*> get_root(const int &n) {
		static std::vector<Z> root{Z::from_raw(1)};
		static std::vector<Z> inv_root{Z::from_raw(1)};
		if (static_cast<int>(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<Z>() {}
	explicit polynomial(const int &n): std::vector<Z>(n) {}
	explicit polynomial(const std::vector<Z> &a): std::vector<Z>(a) {}
	explicit polynomial(const std::initializer_list<Z> &a): std::vector<Z>(a) {}
	template<class _InputIterator, class = std::_RequireInputIter<_InputIterator>>
	explicit polynomial(_InputIterator __first, _InputIterator __last): std::vector<Z>(__first, __last) {}
	template<class F, class = Z(*)(int)>
	explicit polynomial(const int &n, const F &f): std::vector<Z>(n) {
		for (int i = 0; i != n; i ++) {
			(*this)[i] = f(i);
		}
	}
	
	int ord() const {
		int ord = 0;
		while (ord != static_cast<int>(this->size()) && !(*this)[ord]) {
			++ ord;
		}
		return ord;
	}
	int deg() const {
		int deg = static_cast<int>(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<int>(this->size()), k), this->end());
	}
	polynomial mod_xk(const int &k) const {
		return polynomial(this->begin(), this->begin() + std::min(static_cast<int>(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<int>(g.size()) - n), g.end(), r.begin());
		q = (q * r.inv()).mod_xk(n), std::reverse(q.begin(), q.end());
		return q;
	}
	std::pair<polynomial, polynomial> 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<i64>(i) * k_chkmin_sz >= static_cast<i64>(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<class T> inline polynomial pow(T x) const {
		if (x < 0) return inv().pow(-x);
		return pow(x < static_cast<T>(this->size()) ? x : this->size(), x % Mod, x % (Mod - 1));
	}
	std::pair<bool, polynomial> sqrt() const {
		const int i = ord();
		if (i == static_cast<int>(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<class T>
	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<Z>;

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<int, 5> op{};
	op.fill(1);
//	for (int o = 0; o < 5; o ++) {
//		std::cin >> op[o];
//	}
	
	std::array<poly, K> 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<i64, poly> {
		return {N / 2, get2(P, Q[k], N % 2) * Q2[k + 1]};
	};
	auto trans2C = [&](int k, i64 N, const poly& P) -> std::pair<i64, poly> {
		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<i64, poly> {
		return {N / 3, get3(P, Q[k], N % 3) * Q3[k + 1]};
	};
	auto solve = [&](auto&& self, int k, const std::vector<std::pair<i64, poly>>& Fs, const std::vector<std::pair<i64, poly>>& Gs) -> Z {
		if (Fs.empty() && Gs.empty()) {
			return 0;
		}
		
		assert(k < K);
		std::vector<std::pair<i64, poly>> 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<std::pair<i64, poly>>& Ps) {
			std::sort(Ps.begin(), Ps.end(), [&](const std::pair<i64, poly>& x, const std::pair<i64, poly>& y) {
				return x.first < y.first;
			});
			std::vector<std::pair<i64, poly>> 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;
}