Submission #67789700

---

Source Code Expand

#include <bits/stdc++.h>

#include <cstdint>
#include <istream>
#include <ostream>

using std::istream, std::ostream;

template <uint32_t base>
struct Montgomery {
	using i32 = int32_t;
	using u32 = uint32_t;
	using u64 = uint64_t;

	static constexpr u32 mod() {
		return base;
	}

	static constexpr u32 np = []() {
		u32 x = base;
		for (int i = 0; i < 4; ++i) {
			x *= 2u - base * x;
		}
		return -x;
	}();
	static constexpr u32 r2 = -(u64)(base) % base;

	static_assert(base < (1u << 30));
	static_assert(base * np + 1 == 0);

	static u32 reduce(u64 x) {
		return (x + (u64)((u32)x * np) * base) >> 32;
	}

	u32 x;
	Montgomery(): x(0) {}
	constexpr Montgomery(long long y): x(y ? reduce((u64)(y % base + base) * r2) : 0) {}

	Montgomery& operator +=(const Montgomery& ot) {
		if ((i32)(x += ot.x - 2 * base) < 0) {
			x += 2 * base;
		}
		return *this;
	}

	Montgomery& operator -=(const Montgomery& ot) {
		if ((i32)(x -= ot.x) < 0) {
			x += 2 * base;
		}
		return *this;
	}

	Montgomery& operator *=(const Montgomery& ot) {
		x = reduce((u64)x * ot.x);
		return *this;
	}

	Montgomery& operator /=(const Montgomery& ot) {
		return *this *= ot.inverse();
	}

	friend Montgomery operator +(Montgomery a, const Montgomery& b) {
		a += b;
		return a;
	}

	friend Montgomery operator -(Montgomery a, const Montgomery& b) {
		a -= b;
		return a;
	}

	friend Montgomery operator *(Montgomery a, const Montgomery& b) {
		a *= b;
		return a;
	}

	friend Montgomery operator /(Montgomery a, const Montgomery& b) {
		a /= b;
		return a;
	}

	Montgomery operator -() const {
		return Montgomery() - *this;
	}

	u32 get() const {
		u32 res = reduce(x);
		return res < base ? res : res - base;
	}

	u32 operator ()() const {
		return get();
	}

	Montgomery inverse() const {
		return pow(base - 2);
	}

	Montgomery inv() const {
		return inverse();
	}

	Montgomery pow(int64_t p) const {
		if (p < 0) {
			return pow(-p).inverse();
		}
		Montgomery res = 1;
		Montgomery a = *this;
		while (p) {
			if (p & 1) {
				res *= a;
			}
			p >>= 1;
			a *= a;
		}
		return res;
	}

	friend istream& operator >>(istream& istr, Montgomery& m) {
		long long x;
		istr >> x;
		m = Montgomery(x);
		return istr;
	}

	friend ostream& operator <<(ostream& ostr, const Montgomery& m) {
		return ostr << m.get();
	}

	bool operator ==(const Montgomery& ot) const {
		return (x >= base ? x - base : x) == (ot.x >= base ? ot.x - base : ot.x);
	}

	bool operator !=(const Montgomery& ot) const {
		return (x >= base ? x - base : x) != (ot.x >= base ? ot.x - base : ot.x);
	}

	explicit operator int64_t() const {
		return x;
	}

	explicit operator bool() const {
		return x;
	}
};


#include <vector>


using std::vector;

template <int mod>
struct InvfactStuff {
	using Mint = Montgomery<mod>;

	int n;
	vector<Mint> inv, fact, invfact;

	explicit InvfactStuff(int _n): n(_n + 1), inv(_n + 1, 1), fact(_n + 1, 1), invfact(_n + 1, 1) {
		for (int i = 2; i < n; ++i) {
			inv[i] = -inv[mod % i] * (mod / i);
			fact[i] = fact[i - 1] * i;
			invfact[i] = invfact[i - 1] * inv[i];
		}
	}

	Mint C(int n, int k) const {
		if (k < 0 || k > n) {
			return 0;
		}
		assert(n < this->n);
		return fact[n] * invfact[k] * invfact[n - k];
	}

	Mint binom(int n, int k) const {
		return C(n, k);
	}

	Mint factorial(int n) const {
		assert(n < this->n);
		return fact[n];
	}

	Mint inverse_factorial(int n) const {
		assert(n < this->n);
		return invfact[n];
	}

	Mint inverse(int n) const {
		assert(n < this->n);
		return inv[n];
	}

	Mint falling(int n, int k) const {
		if (k > n) {
			return 0;
		}
		assert(n < this->n);
		return fact[n] * invfact[n - k];
	}
};

template <typename Func>
void rec_pythagorean(long long x, long long y, long long z, long long n, const Func& f) {
	if (z > n) {
		return;
	}
	f(x, y, z);
	rec_pythagorean(
		1 * x - 2 * y + 2 * z,
		2 * x - 1 * y + 2 * z,
		2 * x - 2 * y + 3 * z,
		n, f);
	rec_pythagorean(
		1 * x + 2 * y + 2 * z,
		2 * x + 1 * y + 2 * z,
		2 * x + 2 * y + 3 * z,
		n, f);
	rec_pythagorean(
		-1 * x + 2 * y + 2 * z,
		-2 * x + 1 * y + 2 * z,
		-2 * x + 2 * y + 3 * z,
		n, f);
}

template <typename Func>
void for_all_pythagorean_triples(long long n, const Func& f) {
	rec_pythagorean(3, 4, 5, n, f);
}


#include <cassert>
#include <algorithm>
#include <utility>
#include <type_traits>
#include <stdexcept>


#include <memory>	// to define make_unique

#define all(x) (x).begin(), (x).end()
#define make_unique(x) sort((x).begin(), (x).end()); (x).erase(unique((x).begin(), (x).end()), (x).end())
#define imie(x) #x << ": " << x



using li = long long;
using LI = __int128_t;
using ld = long double;
using uint = unsigned int;
using ull = unsigned long long;

template <typename, typename D = void> struct next_size;
template <typename T> using next_size_t = typename next_size<T>::type;
template <typename T> struct tag { using type = T; };

template <> struct next_size<int8_t> : tag<int16_t> {};
template <> struct next_size<int16_t> : tag<int32_t> {};
template <> struct next_size<int32_t> : tag<int64_t> {};
template <typename T> struct next_size<T, std::enable_if_t<4 < sizeof(T) && sizeof(T) <= 8, void>> : tag<__int128_t> {};


ld eps = 1e-8;
void set_eps(ld new_eps) {
	eps = new_eps;
}

template <typename T>
std::enable_if_t<std::is_integral_v<T>, int> sign(T x) {
	return x < 0 ? -1 : x > 0;
}

template <typename T>
std::enable_if_t<std::is_floating_point_v<T>, int> sign(T x) {
	return x < -eps ? -1 : x > eps;
}

#include <cmath>

using std::swap, std::vector;
using std::is_same_v, std::conditional_t;

#define _repeat_cat(a, b) a##b
#define _repeat_helper(ctr, n) for (int _repeat_cat(_mx_, ctr) = n, _repeat_cat(_i_, ctr) = 0; _repeat_cat(_i_, ctr) < _repeat_cat(_mx_, ctr); ++_repeat_cat(_i_, ctr))
#define repeat(n) _repeat_helper(__COUNTER__, n)

LI gcd(LI x, LI y) {
	while (y) {
		x %= y;
		swap(x, y);
	}
	return x;
}

LI lcm(LI x, LI y) {
	return x / gcd(x, y) * y;
}

// LI abs(LI x) {
// 	return x > 0 ? x : -x;
// }

template <typename T>
T pw(T a, li b) {
	T res = 1;
	while (b) {
		if (b & 1ll) {
			res *= a;
		}
		b >>= 1;
		a *= a;
	}
	return res;
}

li pw(li a, li b, int mod) {
	li res = 1;
	while (b) {
		if (b & 1ll) {
			res = res * a % mod;
		}
		b >>= 1;
		a = a * a % mod;
	}
	return res;
}

li mult_big(li x, li y, li mod) {
	li res = x * y - (li)((ld)x * (ld)y / (ld)mod) * mod;
	while (res < 0) {
		res += mod;
	}
	while (res >= mod) {
		res -= mod;
	}
	return res;
}

li pw_big(li a, li b, li mod) {
	li res = 1;
	while (b) {
		if (b & 1ll) {
			res = mult_big(res, a, mod);
		}
		b >>= 1;
		a = mult_big(a, a, mod);
	}
	return res;
}

template <typename T>
conditional_t<is_same_v<T, int>, long long, T> sqr(T x) {
	if constexpr (is_same_v<T, int>) {
		return 1ll * x * x;
	}
	return x * x;
}

template <typename T>
bool remin(T& x, T y) {
	if (x > y) {
		x = y;
		return true;
	} else {
		return false;
	}
}

template <typename T>
bool remax(T& x, T y) {
	if (x < y) {
		x = y;
		return true;
	} else {
		return false;
	}
}

int isqrt(long long n) {
	int x = sqrtl(n);
	while (sqr(x + 1) <= n) {
		++x;
	}
	return x;
}

int icbrt(long long n) {
	int x = pow(n, 1./3);
	while (sqr(x + 1) * (x + 1) <= n) {
		++x;
	}
	return x;
}

int mex(vector<int> a) {
	make_unique(a);
	int res = 0;
	while (res < (int)a.size() && a[res] == res) {
		++res;
	}
	return res;
}

template <typename T, typename Iter>
int find_position(Iter b, Iter e, const T& x) {
	auto it = lower_bound(b, e, x);
	if (it != e && *it == x) {
		return it - b;
	} else {
		return -1;
	}
}
using std::vector, std::pair;
using std::max, std::min, std::swap;
using std::is_convertible_v;
using std::runtime_error;

template <typename outer_type, typename inner_type, int N>
class IFFT {
public:
	using Poly = vector<outer_type>;

	IFFT(): initialized_(false) {}

	~IFFT() {}

	virtual Poly multiply(Poly a, Poly b) {
		if (a.empty() || b.empty()) {
			return {};
		}
		if ((int)a.size() + (int)b.size() > N) {
			Poly result(a.size() + b.size() - 1);
			const int low_len = (max(a.size(), b.size()) + 1) / 2;
			Poly a_low(a.begin(), min(a.begin() + low_len, a.end()));
			Poly a_high(min(a.begin() + low_len, a.end()), a.end());
			Poly b_low(b.begin(), min(b.begin() + low_len, b.end()));
			Poly b_high(min(b.begin() + low_len, b.end()), b.end());

			auto tmp = multiply(a_low, b_low);
			for (int i = 0; i < (int)tmp.size(); ++i) {
				result[i] += tmp[i];
				if (low_len + i < (int)result.size()) {
					result[low_len + i] -= tmp[i];
				}
			}
			tmp = multiply(a_high, b_high);
			for (int i = 0; i < (int)tmp.size(); ++i) {
				result[2 * low_len + i] += tmp[i];
				if (low_len + i < (int)result.size()) {
					result[low_len + i] -= tmp[i];
				}
			}
			for (int i = 0; i < (int)a_high.size(); ++i) {
				a_low[i] += a_high[i];
			}
			for (int i = 0; i < (int)b_high.size(); ++i) {
				b_low[i] += b_high[i];
			}
			tmp = multiply(a_low, b_low);
			for (int i = 0; i < (int)tmp.size(); ++i) {
				result[low_len + i] += tmp[i];
			}
			return result;
		}
		int n = 1;
		while (n < (int)a.size() + (int)b.size()) {
			n *= 2;
		}
		vector<inner_type> ar(n), br(n);
		if constexpr (is_convertible_v<outer_type, inner_type>) {
			copy(all(a), ar.begin());
			copy(all(b), br.begin());
		} else {
			throw runtime_error("please, implement your own child multiply function");
		}
		fft(ar);
		fft(br);
		for (int i = 0; i < (int)ar.size(); ++i) {
			ar[i] *= br[i];
		}
		ifft(ar);
		Poly res((int)a.size() + (int)b.size() - 1);
		assert(res.size() <= ar.size());
		if constexpr (is_convertible_v<inner_type, outer_type>) {
			copy_n(ar.begin(), res.size(), res.begin());
		} else {
			throw runtime_error("please, implement your own child multiply function");
		}
		return res;
	}

	virtual Poly square(const Poly& a) {
		int n = 1;
		while (n < (int)a.size()) {
			n *= 2;
		}
		vector<inner_type> ar(n + n);
		if constexpr (is_convertible_v<outer_type, inner_type>) {
			copy(all(a), ar.begin());
		} else {
			throw runtime_error("please, implement your own child square function");
		}
		fft(ar);
		for (int i = 0; i < (int)ar.size(); ++i) {
			ar[i] *= ar[i];
		}
		ifft(ar);
		Poly res((int)a.size() + (int)a.size() - 1);
		assert(res.size() <= ar.size());
		if constexpr (is_convertible_v<inner_type, outer_type>) {
			copy_n(ar.begin(), res.size(), res.begin());
		} else {
			throw runtime_error("please, implement your own child square function");
		}
		return res;
	}

	virtual Poly pow(const Poly& a, int k) {
		int n = 1;
		while (n < (int)a.size() * k - k + 1) {
			n *= 2;
		}
		vector<inner_type> ar(n);
		if constexpr (is_convertible_v<outer_type, inner_type>) {
			copy(all(a), ar.begin());
		} else {
			throw runtime_error("please, implement your own child pow function");
		}
		fft(ar);
		for (int i = 0; i < (int)ar.size(); ++i) {
			ar[i] = pw(ar[i], k);
		}
		ifft(ar);
		Poly res((int)a.size() * k - k + 1);
		assert(res.size() <= ar.size());
		if constexpr (is_convertible_v<inner_type, outer_type>) {
			copy_n(ar.begin(), res.size(), res.begin());
		} else {
			throw runtime_error("please, implement your own child pow function");
		}
		return res;
	}

	Poly inverse(const Poly& a, int prec) {
		assert(!a.empty());
		assert(a[0] != 0);
		Poly b = {1 / a[0]};
		for (int len = 1; len < prec; len *= 2) {
			auto tmp = multiply(b, b);
			if ((int)tmp.size() > prec) {
				tmp.resize(prec);
			}
			tmp = multiply(tmp, Poly{a.begin(), a.begin() + min(2 * len, (int)a.size())});
			tmp.resize(2 * len);
			for (int i = 0; i < len; ++i) {
				tmp[i] = b[i] + b[i] - tmp[i];
				tmp[len + i] = -tmp[len + i];
			}
			b.swap(tmp);
		}
		b.resize(prec);
		return b;
	}

	Poly derivative(Poly a) {
		if (a.empty()) {
			return a;
		}
		a.erase(a.begin());
		for (int i = 0; i < (int)a.size(); ++i) {
			a[i] *= i + 1;
		}
		return a;
	}

	Poly primitive(Poly a) {
		a.insert(a.begin(), 0);
		for (int i = 1; i < (int)a.size(); ++i) {
			a[i] /= i;
		}
		return a;
	}

	Poly log(const Poly& a, int prec) {
		assert(!a.empty());
		assert(a[0] == 1);
		auto res = primitive(multiply(derivative(a), inverse(a, prec)));
		res.resize(prec);
		return res;
	}

	Poly exp(const Poly& a, int prec) {
		assert(!a.empty());
		assert(a[0] == 0);
		Poly b = {1};
		for (int len = 1; len < prec; len *= 2) {
			auto tmp = Poly{a.begin(), a.begin() + min(2 * len, (int)a.size())};
			tmp.resize(2 * len);
			tmp[0] += 1;
			auto l = log(b, 2 * len);
			for (int i = 0; i < 2 * len; ++i) {
				tmp[i] -= l[i];
			}
			b = multiply(tmp, b);
			b.resize(2 * len);
		}
		b.resize(prec);
		return b;
	}

	pair<Poly, Poly> divmod(Poly a, Poly b) {
		assert(!b.empty());
		assert(b.back() != 0);
		if (a.size() < b.size()) {
			return {{0}, a};
		}
		reverse(all(a));
		reverse(all(b));
		auto q = inverse(b, a.size() - b.size() + 1);
		q = multiply(a, q);
		q.resize(a.size() - b.size() + 1);
		reverse(all(q));
		reverse(all(a));
		reverse(all(b));
		Poly r(b.size() - 1);
		auto bq = multiply(b, q);
		for (int i = 0; i < (int)r.size(); ++i) {
			r[i] = a[i] - bq[i];
		}
		return {q, r};
	}

	struct ProductTree {
		int n;
		vector<Poly> a;
		vector<outer_type> x;
		IFFT* owner;

		ProductTree(const vector<outer_type>& _x, IFFT* that): x(_x), owner(that) {
			n = 1;
			while (n < (int)x.size()) {
				n *= 2;
			}
			a.resize(n + n);
			outer_type one{1};
			for (int i = 0; i < (int)x.size(); ++i) {
				a[n + i] = {-x[i], one};
			}
			for (int i = n - 1; i > 0; --i) {
				if (a[2 * i].empty()) {
					continue;
				} else if (a[2 * i + 1].empty()) {
					a[i] = a[2 * i];
				} else {
					a[i] = owner->multiply(a[2 * i], a[2 * i + 1]);
				}
			}
		}

		const Poly& top() {
			return a[1];
		}

		void rec_multipoint(int v, int l, int r, Poly p, vector<outer_type>& ans) {
			if (l >= (int)x.size()) {
				return;
			}
			p = owner->divmod(p, a[v]).second;
			if (r <= l + 64) {
				for (int i = l; i < r && i < (int)x.size(); ++i) {
					outer_type& res = ans[i] = 0;
					for (int j = (int)p.size() - 1; j >= 0; --j) {
						res = res * x[i] + p[j];
					}
				}
				return;
			}
			int m = (l + r) / 2;
			rec_multipoint(v + v, l, m, p, ans);
			rec_multipoint(v + v + 1, m, r, p, ans);
		}

		vector<outer_type> multipoint(const Poly& p) {
			vector<outer_type> ans(x.size());
			rec_multipoint(1, 0, n, p, ans);
			return ans;
		}
	};

	Poly multipoint(Poly p, const vector<outer_type>& x) {
		ProductTree tree(x, this);
		return tree.multipoint(p);
	}

	Poly interpolate(const vector<outer_type>& x, const vector<outer_type>& y) {
		ProductTree tree(x, this);
		auto denom = tree.multipoint(derivative(tree.top()));
		auto inter = [&](const auto& self, int v, int l, int r) -> Poly {
			if (l >= (int)x.size()) {
				return {};
			}
			if (l + 1 == r) {
				return {y[l] / denom[l]};
			}
			int m = (l + r) / 2;
			auto left = self(self, v + v, l, m);
			auto right = self(self, v + v + 1, m, r);
			if (right.empty()) {
				return left;
			}
			left = multiply(left, tree.a[v + v + 1]);
			right = multiply(right, tree.a[v + v]);
			if (left.size() < right.size()) {
				left.resize(right.size());
			}
			for (int i = 0; i < (int)right.size(); ++i) {
				left[i] += right[i];
			}
			return left;
		};
		return inter(inter, 1, 0, tree.n);
	}

protected:
	static constexpr int L = 31 - __builtin_clz(N);
	static_assert(!(N & (N - 1)));
	vector<inner_type> angles;
	vector<int> bitrev;
	bool initialized_;

	void initialize() {
		fill_angles();
		fill_bitrev();
		initialized_ = true;
	}

	virtual void fill_angles() = 0;

	void fill_bitrev() {
		bitrev.assign(N, 0);
		for (int i = 0; i < N; ++i) {
			for (int j = 0; j < L; ++j) {
				bitrev[i] = (bitrev[i] << 1) | !!(i & (1 << j));
			}
		}
	}

	void butterfly(vector<inner_type>& a) {
		if (!initialized_) {
			initialize();
		}
		const int n = a.size();
		assert(!(n & (n - 1)));
		const int l = __builtin_ctz(n);

		for (int i = 0; i < n; ++i) {
			int j = revbit(i, l);
			if (j > i) {
				swap(a[i], a[j]);
			}
		}
	}

	int revbit(int num, int len) const {
		return bitrev[num] >> (L - len);
	}

	virtual void fft(vector<inner_type>& a) {
		if (!initialized_) {
			initialize();
		}

		const int n = a.size();
		assert(!(n & (n - 1)));
		butterfly(a);

		for (int len = 1; len < n; len *= 2) {
			for (int start = 0; start < n; start += 2 * len) {
				for (int i = 0; i < len; ++i) {
					const auto w = angles[N / 2 / len * i];
					const auto x = a[start + i], y = a[start + len + i] * w;
					a[start + i] = x + y;
					a[start + len + i] = x - y;
				}
			}
		}
	}

	void ifft(vector<inner_type>& a) {
		fft(a);
		inner_type to_div = inner_type(1) / a.size();
		for (auto& x : a) {
			x *= to_div;
		}
		reverse(1 + all(a));
	}
};

template <int mod, int N = (1 << __builtin_ctz(mod - 1))>
class NTT : public IFFT<Montgomery<mod>, Montgomery<mod>, N> {
	using Mint = Montgomery<mod>;
protected:
	void fill_angles() {
		vector<int> primes;
		for (int x = mod - 1, p = 2; x > 1; ++p) {
			if (p * p > x) {
				primes.push_back(x);
				break;
			}
			if (x % p == 0) {
				primes.push_back(p);
				while (x % p == 0) {
					x /= p;
				}
			}
		}
		auto isPrimitiveRoot = [&](Mint g) {
			for (int p : primes) {
				if (g.pow(mod / p) == 1) {
					return false;
				}
			}
			return true;
		};
		int g = 2;
		while (!isPrimitiveRoot(g)) {
			++g;
		}
		g = Mint(g).pow(mod / N).get();

		this->angles.assign(N, 1);
		for (int i = 1; i < N; ++i) {
			this->angles[i] = this->angles[i - 1] * g;
		}
		assert(this->angles.back() * g == 1);
	}

	void ntt(vector<Mint>& a) {
		fft(a);
	}
};

#define all(x) (x).begin(), (x).end()

using namespace std;

inline int nxt() {
	int x;
	cin >> x;
	return x;
}

constexpr int mod = 998244353;
using Mint = Montgomery<mod>;
NTT<mod, (1 << 20)> ntt;

struct Sumh {
	int m;
	int k;
	int n;
	bool count_mod;
	vector<Mint> values;

	explicit Sumh(int _m, int _k, int maxn): m(_m), k(_k), n(0) {
		count_mod = (m < 1000);
		if (count_mod) {
			values.assign(m, 0);
			values[0] = 1;
		} else {
			InvfactStuff<mod> stuff(maxn);
			vector<Mint> sh_hat(maxn);
			for (int i = 0; i < maxn; ++i) {
				int n = (2 * i - k + m) % m;
				for (; n < maxn; n += m) {
					if (n >= i) {
						sh_hat[n] += stuff.inverse_factorial(i) * stuff.inverse_factorial(n - i);
					}
				}
			}
			vector<Mint> e(maxn);
			for (int i = 0; i < maxn; ++i) {
				e[i] = stuff.inverse_factorial(i);
			}
			values = ntt.multiply(sh_hat, e);
			values.resize(maxn);
			for (int i = 0; i < maxn; ++i) {
				values[i] *= stuff.factorial(i);
			}
		}
	}

	void inc_n() {
		n += 1;
		if (count_mod) {
			Mint fst = values[0];
			Mint pr = values.back();
			for (int i = 0; i < m; ++i) {
				pr = exchange(values[i], pr + values[i]);
				values[i] += (i + 1 < m) ? values[i + 1] : fst;
			}
		} else {
			// values.push_back(0);
			// Mint old = values[0];
			// values[0] += 2 * values[1];
			// for (int i = 1; i < (int)values.size(); ++i) {
			// 	old = exchange(values[i], values[i] + old);
			// 	if (i + 1 < (int)values.size()) {
			// 		values[i] += values[i + 1];
			// 	}
			// }
		}
	}

	void advance_n(int next_n) {
		while (n < next_n) {
			inc_n();
		}
	}

	Mint get(int kk) const {
		assert(kk == k);
		if (count_mod) {
			kk = abs(kk) % m;
			return (kk < m) ? values[kk] : 0;
		} else {
			return values[n];
			// Mint ans = 0;
			// k += n;
			// k %= m;
			// if (k < 0) {
			// 	k += m;
			// }
			// while (k <= 2 * n) {
			// 	ans += values[abs(k - n)];
			// 	k += m;
			// }
			// return ans;
		}
	}
};

void solve() {
	int h = nxt(), w = nxt(), t = nxt(), a = nxt() - 1, b = nxt() - 1, c = nxt() - 1, d = nxt() - 1;
	InvfactStuff<mod> stuff(t);

	vector<Mint> curh = {1};
	int curn = 0;

	vector<Sumh> mp;

	auto inc_curn = [&]() {
		++curn;
		curh.push_back(0);
		Mint old = curh[0];
		curh[0] += 2 * curh[1];
		for (int i = 1; i < (int)curh.size(); ++i) {
			old = exchange(curh[i], curh[i] + old);
			if (i + 1 < (int)curh.size()) {
				curh[i] += curh[i + 1];
			}
		}
	};
	[[maybe_unused]] auto hlp = [&](int n, int k) -> Mint {
		// k = abs(k);
		// Mint ans = 0;
		// for (int i = 0; 2 * i + k <= n; ++i) {
		// 	ans += stuff.inverse_factorial(i) * stuff.inverse_factorial(i + k) * stuff.inverse_factorial(n - 2 * i - k);
		// }
		// return ans * stuff.factorial(n);

		assert(n >= curn);
		while (curn < n) {
			inc_curn();
		}
		k = abs(k);
		return k <= n ? curh[k] : 0;
	};
	auto sumh = [&](int n, int k, int m) {
		int idx = 0;
		while (idx < (int)mp.size() && (mp[idx].m != m || mp[idx].k != k)) {
			++idx;
		}
		if (idx == (int)mp.size()) {
			mp.emplace_back(m, k, t + 10);
		}
		mp[idx].advance_n(n);
		return mp[idx].get(k);

		// Mint ans = 0;
		// k += n;
		// k %= m;
		// if (k < 0) {
		// 	k += m;
		// }
		// while (k <= 2 * n) {
		// 	ans += hlp(n, k - n);
		// 	k += m;
		// }
		// return ans;
	};
	auto g = [&](int a, int b, int s, int n) {
		return sumh(n, b - a, 2 * (s + 1)) - sumh(n, -b - 2 - a, 2 * (s + 1));
	};
	auto f = [&](int n) {
		return g(a, c, h, n) * g(b, d, w, n);
	};

	Mint ans = 0;
	for (int k = t; k >= 0; --k) {
		ans += stuff.C(t, k) * (k % 2 ? -1 : 1) * f(t - k);
	}
	cout << ans << "\n";
}

int main() {
	ios_base::sync_with_stdio(false);
	cin.tie(nullptr);

	int t = 1; // nxt();
	while (t--) {
		solve();
	}

	return 0;
}

Submission Info

Judge Result