//https://qandori.github.io/cp-library/test/dijkstra_bigint_yuki.test.cpp.html
#line 1 "test/dijkstra_bigint_yuki.test.cpp"
#define PROBLEM "https://yukicoder.me/problems/no/3669"
#line 1 "graph/dijkstra.hpp"



#line 1 "components/RadixHeap.hpp"



#include <algorithm>
#include <bit>
#include <cassert>
#include <cstddef>
#include <cstdint>
#include <limits>
#include <type_traits>
#include <utility>
#include <vector>

template <class Key = std::uint64_t, class T = int>
class RadixHeap {
	static_assert(std::is_integral_v<Key> && !std::is_same_v<Key, bool> && sizeof(Key) <= 8);
	using U = std::make_unsigned_t<Key>;
	using State = std::pair<Key, T>;
	std::vector<State> bucket[std::numeric_limits<U>::digits + 1];
	Key lower;
	std::size_t count = 0;
	int index(Key key) const { return std::bit_width(U(U(key) ^ U(lower))); }
	void refill() {
		assert(count);
		if (!bucket[0].empty()) return;
		int i = 1;
		while (bucket[i].empty()) ++i;
		lower = bucket[i][0].first;
		for (const auto &s : bucket[i]) lower = std::min(lower, s.first);
		for (auto &s : bucket[i]) bucket[index(s.first)].push_back(std::move(s));
		bucket[i].clear();
	}

  public:
	explicit RadixHeap(Key minimum = 0) : lower(minimum) {}
	RadixHeap(const RadixHeap &) = default;
	RadixHeap &operator=(const RadixHeap &) = default;
	RadixHeap(RadixHeap &&other) noexcept : lower(other.lower), count(std::exchange(other.count, 0)) {
		for (int i = 0; i <= std::numeric_limits<U>::digits; ++i) bucket[i] = std::move(other.bucket[i]);
	}
	RadixHeap &operator=(RadixHeap &&other) noexcept {
		if (this != &other) {
			lower = other.lower;
			count = std::exchange(other.count, 0);
			for (int i = 0; i <= std::numeric_limits<U>::digits; ++i) bucket[i] = std::move(other.bucket[i]);
		}
		return *this;
	}

	bool empty() const { return count == 0; }
	std::size_t size() const { return count; }
	Key last() const { return lower; }
	void clear(Key minimum = 0) {
		for (auto &b : bucket) b.clear();
		lower = minimum;
		count = 0;
	}
	void push(State s) {
		assert(s.first >= lower);
		bucket[index(s.first)].push_back(std::move(s));
		++count;
	}
	const State &top() {
		refill();
		return bucket[0].back();
	}
	void pop() {
		refill();
		bucket[0].pop_back();
		--count;
	}
	State pop_value() {
		refill();
		State s = std::move(bucket[0].back());
		bucket[0].pop_back();
		--count;
		return s;
	}
};


#line 9 "graph/dijkstra.hpp"
#include <functional>
#line 11 "graph/dijkstra.hpp"
#include <queue>
#line 15 "graph/dijkstra.hpp"

inline constexpr std::int64_t DIJKSTRA_INF = 1'000'000'000'000'000'000;

template <class T>
std::vector<T> dijkstra(const std::vector<std::vector<std::pair<int, T>>> &g, int st, T inf = static_cast<T>(std::min<std::int64_t>(DIJKSTRA_INF, std::numeric_limits<T>::max()))) {
	assert(g.size() <= static_cast<std::size_t>(std::numeric_limits<int>::max()));
	assert(st >= 0 && static_cast<std::size_t>(st) < g.size());
	assert(inf > 0);
	std::vector<T> d(g.size(), inf);
	d[st] = 0;
	using State = std::pair<T, int>;
	constexpr bool radix = std::is_integral_v<T> && sizeof(T) <= 8;
	std::conditional_t<radix, RadixHeap<T>, std::priority_queue<State, std::vector<State>, std::greater<State>>> pq;
	pq.push({0, st});
	while (!pq.empty()) {
		const auto [cost, v] = [&] {
			if constexpr (radix) return pq.pop_value();
			else {
				State s = pq.top();
				pq.pop();
				return s;
			}
		}();
		if (d[v] != cost) continue;
		for (const auto &[to, c] : g[v]) {
			assert(to >= 0 && static_cast<std::size_t>(to) < g.size());
			assert(c >= 0);
			if (c >= inf - cost) continue;
			T next = cost + c;
			if (next < d[to]) d[to] = next, pq.push({next, to});
		}
	}
	for (T &val : d)
		if (val == inf) val = -1;
	return d;
}


#line 1 "math/bigint.hpp"



#line 1 "components/ntt.hpp"



#line 7 "components/ntt.hpp"
#include <initializer_list>
#include <memory>
#include <new>
#line 13 "components/ntt.hpp"
#if defined(__ARM_NEON)
#include <arm_neon.h>
#endif

#if defined(__x86_64__) || defined(__i386__)
#include <immintrin.h>
#endif

#line 1 "math/mint.hpp"


#line 4 "math/mint.hpp"
#include <istream>
#include <ostream>
#line 7 "math/mint.hpp"

#if defined(__GNUC__) && !defined(__clang__) && (defined(__x86_64__) || defined(__i386__))
#pragma GCC push_options
#pragma GCC optimize("O3")
#endif

namespace fps_detail {
template <class I>
constexpr bool mint_integral = std::is_integral_v<I>
#if defined(__SIZEOF_INT128__)
							   || std::is_same_v<I, __int128> || std::is_same_v<I, unsigned __int128>
#endif
	;
struct Mint {
	using mint = Mint;
	unsigned v = 0;
	static constexpr unsigned mod() { return 998244353; }
	static unsigned reduce(unsigned long long x) {
		return (x + (unsigned long long)(unsigned(x) * 998244351u) * mod()) >> 32;
	}
	Mint() = default;
	template <class I, std::enable_if_t<mint_integral<I>, int> = 0>
	Mint(I x) {
		long long r = x % (long long)mod();
		if (r < 0) r += mod();
		v = reduce((unsigned long long)r * 932051910);
	}
	static Mint raw(unsigned x) {
		Mint a;
		a.v = reduce((unsigned long long)x * 932051910);
		return a;
	}
	static Mint word(unsigned x) {
		Mint a;
		a.v = x;
		return a;
	}
	unsigned val() const {
		unsigned r = reduce(v);
		return r < mod() ? r : r - mod();
	}
	friend std::istream &operator>>(std::istream &is, Mint &x) {
		long long value;
		if (is >> value) x = Mint(value);
		return is;
	}
	friend std::ostream &operator<<(std::ostream &os, const Mint &x) {
		return os << x.val();
	}
	Mint &operator+=(Mint b) {
		v += b.v;
		if (v >= 2 * mod()) v -= 2 * mod();
		return *this;
	}
	Mint &operator-=(Mint b) {
		v -= b.v;
		if (v >> 31) v += 2 * mod();
		return *this;
	}
	Mint &operator*=(Mint b) {
		v = reduce((unsigned long long)v * b.v);
		return *this;
	}
	Mint &operator/=(Mint b) { return *this *= b.inv(); }
	Mint &operator++() { return *this += Mint(1); }
	Mint &operator--() { return *this -= Mint(1); }
	Mint operator++(int) {
		Mint a = *this;
		++*this;
		return a;
	}
	Mint operator--(int) {
		Mint a = *this;
		--*this;
		return a;
	}
	Mint operator+() const { return *this; }
	Mint operator-() const { return word(v ? 2 * mod() - v : 0); }
	Mint operator+(Mint b) const { return Mint(*this) += b; }
	Mint operator-(Mint b) const { return Mint(*this) -= b; }
	Mint operator*(Mint b) const { return Mint(*this) *= b; }
	Mint operator/(Mint b) const { return Mint(*this) /= b; }
	bool operator==(Mint b) const { return (v < mod() ? v : v - mod()) == (b.v < mod() ? b.v : b.v - mod()); }
	bool operator!=(Mint b) const { return !(*this == b); }
	template <class I, std::enable_if_t<mint_integral<I>, int> = 0>
	friend Mint operator+(I a, Mint b) { return Mint(a) + b; }
	template <class I, std::enable_if_t<mint_integral<I>, int> = 0>
	friend Mint operator-(I a, Mint b) { return Mint(a) - b; }
	template <class I, std::enable_if_t<mint_integral<I>, int> = 0>
	friend Mint operator*(I a, Mint b) { return Mint(a) * b; }
	template <class I, std::enable_if_t<mint_integral<I>, int> = 0>
	friend Mint operator/(I a, Mint b) { return Mint(a) / b; }
	template <class I, std::enable_if_t<mint_integral<I>, int> = 0>
	friend bool operator==(I a, Mint b) { return Mint(a) == b; }
	template <class I, std::enable_if_t<mint_integral<I>, int> = 0>
	friend bool operator!=(I a, Mint b) { return Mint(a) != b; }
	Mint pow(long long n) const {
		assert(n >= 0);
		if (!n) return Mint(1);
		Mint r(1), a = *this;
		for (; n > 1; n >>= 1, a *= a)
			if (n & 1) r *= a;
		return r * a;
	}
	[[gnu::noinline]] Mint inv() const {
		assert(*this != Mint(0));
		Mint a2 = *this * *this, a3 = a2 * *this;
		Mint a7 = a3 * a3 * *this;
		Mint a31 = a7;
		for (int i = 0; i < 2; ++i) a31 *= a31;
		a31 *= a3;
		Mint r = a7;
		for (int i = 0; i < 3; ++i) r *= r;
		r *= a3;
		r *= r;
		for (int j = 0; j < 4; ++j) {
			for (int i = 0; i < 5; ++i) r *= r;
			r *= a31;
		}
		for (int i = 0; i < 3; ++i) r *= r;
		return r * a7;
	}
};
} // namespace fps_detail

using fps_mint = fps_detail::Mint;
using mint = fps_mint;

#if defined(__GNUC__) && !defined(__clang__) && (defined(__x86_64__) || defined(__i386__))
#pragma GCC pop_options
#endif


#line 22 "components/ntt.hpp"

#if defined(__GNUC__) && !defined(__clang__) && (defined(__x86_64__) || defined(__i386__))
#pragma GCC push_options
#pragma GCC optimize("O3")
#endif

namespace fps_detail {
struct NTT {
	using T = Mint;
	typedef unsigned int U __attribute__((may_alias));
	using V = std::vector<unsigned int>;
	static constexpr U mod = 998244353, twice = 2 * mod, one = 301989884;
	V roots, inverse_roots, twists, inverse_twists;
	V root_quotients, inverse_root_quotients;
	mutable V cube_roots, inverse_cube_roots, cube_quotients, inverse_cube_quotients;
	void prepare_cubes(int n) const {
		int begin = cube_roots.size();
		if (begin >= n) return;
		cube_roots.resize(n);
		inverse_cube_roots.resize(n);
		cube_quotients.resize(n);
		inverse_cube_quotients.resize(n);
		for (int k = begin; k < n; ++k) {
			U x = (unsigned long long)roots[k] * roots[2 * k] % mod;
			U y = (unsigned long long)inverse_roots[k] * inverse_roots[2 * k] % mod;
			cube_roots[k] = x;
			inverse_cube_roots[k] = y;
			cube_quotients[k] = quotient(x);
			inverse_cube_quotients[k] = quotient(y);
		}
	}
	static U reduce(unsigned long long x) {
		return (x + (unsigned long long)(U(x) * 998244351u) * mod) >> 32;
	}
	static U mul(U a, U b) { return reduce((unsigned long long)a * b); }
	static U add(U a, U b) {
		U c = a + b;
		return c < twice ? c : c - twice;
	}
	static U sub(U a, U b) {
		U c = a - b;
		return c < twice ? c : c + twice;
	}
	static U in(U a) { return mul(a, 932051910); }
	static U out(U a) {
		a = reduce(a);
		return a < mod ? a : a - mod;
	}
	static U quotient(U r) { return ((unsigned long long)r << 32) / mod; }
	static U fixed_mul(U a, U r, U q) {
		return a * r - U(((unsigned long long)a * q) >> 32) * mod;
	}
	NTT(int size, bool composition = true)
		: roots(size, one),
		  inverse_roots(size, one),
		  twists(composition ? size : 0, one),
		  inverse_twists(composition ? size : 0, one),
		  root_quotients(size),
		  inverse_root_quotients(size) {
		for (int h = 1; h < size; h <<= 1) {
			U r = in(T(3).pow((mod - 1) / (4 * h)).val());
			U ir = in(T(3).pow(mod - 1 - (mod - 1) / (4 * h)).val());
			U t = in(T(3).pow((mod - 1) / (2 * h)).val());
			U it = in(T(3).pow(mod - 1 - (mod - 1) / (2 * h)).val());
			U a = in(mod - (mod - 1) / h), b = a;
			for (int i = 0; i < h; ++i) {
				roots[h + i] = mul(roots[i], r);
				inverse_roots[h + i] = mul(inverse_roots[i], ir);
				if (composition) twists[h + i] = a, inverse_twists[h + i] = b;
				a = mul(a, t);
				b = mul(b, it);
			}
		}
		for (int i = 0; i < size; ++i) {
			root_quotients[i] = (roots[i] < mod ? roots[i] : roots[i] - mod) * 998244351u;
			inverse_root_quotients[i] = (inverse_roots[i] < mod ? inverse_roots[i] : inverse_roots[i] - mod) * 998244351u;
			roots[i] = out(roots[i]);
			inverse_roots[i] = out(inverse_roots[i]);
			if (composition) twists[i] = out(twists[i]), inverse_twists[i] = out(inverse_twists[i]);
		}
	}

#if defined(__ARM_NEON)
	using I = uint32x4_t;
	static I mul4(I a, I b) {
		uint64x2_t x = vmull_u32(vget_low_u32(a), vget_low_u32(b));
		uint64x2_t y = vmull_u32(vget_high_u32(a), vget_high_u32(b));
		auto reduce2 = [](uint64x2_t z) {
			uint32x2_t q = vmul_u32(vmovn_u64(z), vdup_n_u32(998244351));
			return vshrn_n_u64(vaddq_u64(z, vmull_u32(q, vdup_n_u32(mod))), 32);
		};
		return vcombine_u32(reduce2(x), reduce2(y));
	}
	static I add4(I a, I b) {
		I c = vaddq_u32(a, b);
		return vminq_u32(c, vsubq_u32(c, vdupq_n_u32(twice)));
	}
	static I sub4(I a, I b) {
		I c = vsubq_u32(a, b);
		return vminq_u32(c, vaddq_u32(c, vdupq_n_u32(twice)));
	}
	static I fixed_mul4(I a, I r, I q) {
		uint64x2_t x = vmull_u32(vget_low_u32(a), vget_low_u32(q));
		uint64x2_t y = vmull_u32(vget_high_u32(a), vget_high_u32(q));
		I t = vcombine_u32(vshrn_n_u64(x, 32), vshrn_n_u64(y, 32));
		return vmlsq_u32(vmulq_u32(a, r), t, vdupq_n_u32(mod));
	}
	template <bool inverse>
	static void fixed_butterfly4(I &x, I &y, I r, I q, bool unit = false) {
		if constexpr (!inverse)
			if (!unit) y = fixed_mul4(y, r, q);
		I z = sub4(x, y);
		x = add4(x, y);
		y = inverse && !unit ? fixed_mul4(z, r, q) : z;
	}
#endif
#if defined(__x86_64__) || defined(__i386__)
	static bool vectorized() {
#if defined(__AVX2__)
		return true;
#else
		static const bool available = __builtin_cpu_supports("avx2");
		return available;
#endif
	}
	using I8 = __m256i;
	__attribute__((target("avx2"))) static I8 load8(const U *a) {
		return _mm256_loadu_si256(reinterpret_cast<const I8 *>(a));
	}
	__attribute__((target("avx2"))) static void store8(U *a, I8 x) {
		_mm256_storeu_si256(reinterpret_cast<I8 *>(a), x);
	}
	__attribute__((target("avx2"))) static I8 mul8(I8 a, I8 b) {
		I8 x = _mm256_mul_epu32(a, b);
		I8 y = _mm256_mul_epu32(_mm256_srli_epi64(a, 32), _mm256_srli_epi64(b, 32));
		I8 r = _mm256_set1_epi32(998244351), p = _mm256_set1_epi32(mod);
		x = _mm256_add_epi64(x, _mm256_mul_epu32(_mm256_mul_epu32(x, r), p));
		y = _mm256_add_epi64(y, _mm256_mul_epu32(_mm256_mul_epu32(y, r), p));
		return _mm256_blend_epi32(_mm256_srli_epi64(x, 32), y, 0xaa);
	}
	__attribute__((target("avx2"))) static I8 fixed_mul8(I8 a, I8 r, I8 q) {
		I8 x = _mm256_srli_epi64(_mm256_mul_epu32(a, q), 32);
		I8 y = _mm256_mul_epu32(_mm256_srli_epi64(a, 32), q);
		I8 t = _mm256_blend_epi32(x, y, 0xaa);
		return _mm256_sub_epi32(_mm256_mullo_epi32(a, r),
								_mm256_mullo_epi32(t, _mm256_set1_epi32(mod)));
	}
	__attribute__((target("avx2"))) static I8 fixed_even8(I8 a, I8 r, I8 q) {
		I8 t = _mm256_srli_epi64(_mm256_mul_epu32(a, q), 32);
		return _mm256_sub_epi64(_mm256_mul_epu32(a, r),
								_mm256_mul_epu32(t, _mm256_set1_epi32(mod)));
	}
	__attribute__((target("avx2"))) static I8 add8(I8 a, I8 b) {
		I8 x = _mm256_add_epi32(a, b);
		return _mm256_min_epu32(x, _mm256_sub_epi32(x, _mm256_set1_epi32(twice)));
	}
	__attribute__((target("avx2"))) static I8 sub8(I8 a, I8 b) {
		I8 x = _mm256_sub_epi32(a, b);
		return _mm256_min_epu32(x, _mm256_add_epi32(x, _mm256_set1_epi32(twice)));
	}

	template <int stride>
	__attribute__((target("avx2"))) static I8 rates8(const V &values, int k) {
		if constexpr (stride == 1)
			return _mm256_cvtepu32_epi64(
				_mm_loadu_si128(reinterpret_cast<const __m128i *>(values.data() + 4 * k)));
		if constexpr (stride == 2)
			return _mm256_permute4x64_epi64(
				_mm256_cvtepu32_epi64(
					_mm_loadl_epi64(reinterpret_cast<const __m128i *>(values.data() + 2 * k))),
				0x50);
		if constexpr (stride == 4) return _mm256_set1_epi32(values[k]);
	}
	template <int stride, bool inverse>
	__attribute__((target("avx2"))) I8 small8(I8 z, int k) const {
		if constexpr (stride == 4) {
			if (k == 0) {
				I8 x = _mm256_permute2x128_si256(z, z, 0x00),
				   y = _mm256_permute2x128_si256(z, z, 0x11);
				return _mm256_blend_epi32(add8(x, y), sub8(x, y), 0xf0);
			}
		}
		const V &w = inverse ? inverse_roots : roots;
		const V &quot = inverse ? inverse_root_quotients : root_quotients;
		I8 r = rates8<stride>(w, k), q = rates8<stride>(quot, k);
		if constexpr (inverse) {
			if constexpr (stride == 2) z = _mm256_shuffle_epi32(z, 0xd8);
			if constexpr (stride == 4)
				z = _mm256_permutevar8x32_epi32(z, _mm256_setr_epi32(0, 4, 1, 5, 2, 6, 3, 7));
			I8 y = _mm256_srli_epi64(z, 32), x = add8(z, y);
			I8 d = fixed_even8(sub8(z, y), r, q);
			z = _mm256_blend_epi32(x, _mm256_slli_epi64(d, 32), 0xaa);
		} else {
			if constexpr (stride == 1) z = _mm256_shuffle_epi32(z, 0xb1);
			if constexpr (stride == 2) z = _mm256_shuffle_epi32(z, 0x72);
			if constexpr (stride == 4)
				z = _mm256_permutevar8x32_epi32(z, _mm256_setr_epi32(4, 0, 5, 1, 6, 2, 7, 3));
			I8 y = fixed_even8(z, r, q);
			I8 x = _mm256_blend_epi32(_mm256_srli_epi64(z, 32), z, 0xaa);
			y = _mm256_blend_epi32(
				y, _mm256_sub_epi32(_mm256_set1_epi32(twice), _mm256_slli_epi64(y, 32)), 0xaa);
			z = add8(x, y);
		}
		if constexpr (stride == 2) z = _mm256_shuffle_epi32(z, 0xd8);
		if constexpr (stride == 4)
			z = _mm256_permutevar8x32_epi32(z, _mm256_setr_epi32(0, 2, 4, 6, 1, 3, 5, 7));
		return z;
	}
	template <bool inverse>
	__attribute__((target("avx2"))) void fft8(U *a, int length, int limit = 0) const {
		if (!limit) limit = length;
		if constexpr (!inverse) {
			int h = limit / 2;
			for (; h >= 16; h >>= 2) radix4_8<false>(a, h / 2, length);
			if (h == 8) {
				if (limit >= 32)
					for (int i = 0; i < length; i += 8) {
						I8 x = load8(a + i), p = _mm256_set1_epi32(twice);
						store8(a + i, _mm256_min_epu32(x, _mm256_sub_epi32(x, p)));
					}
				stage8<false>(a, h, length);
			}
		}
		for (int j = 0; j < length; j += 8) {
			I8 z = load8(a + j);
			if constexpr (!inverse)
				if ((limit >= 32) && (__builtin_ctz((unsigned int)limit) & 1)) {
					I8 p = _mm256_set1_epi32(twice);
					z = _mm256_min_epu32(z, _mm256_sub_epi32(z, p));
				}
			if constexpr (inverse)
				z = small8<4, true>(small8<2, true>(small8<1, true>(z, j / 8), j / 8), j / 8);
			else
				z = small8<1, false>(small8<2, false>(small8<4, false>(z, j / 8), j / 8),
									 j / 8);
			store8(a + j, z);
		}
		if constexpr (inverse) {
			int h = 8;
			for (; 4 * h <= limit; h <<= 2) radix4_8<true>(a, h, length);
			if (h < limit) stage8<true>(a, h, length);
		}
	}
	template <bool inverse>
	__attribute__((target("avx2"))) void radix4_8(U *a, int stride, int length) const {
		const V &w = inverse ? inverse_roots : roots;
		const V &q = inverse ? inverse_root_quotients : root_quotients;
		prepare_cubes(length / (4 * stride));
		const V &w3 = inverse ? inverse_cube_roots : cube_roots;
		const V &q3s = inverse ? inverse_cube_quotients : cube_quotients;
		I8 p = _mm256_set1_epi32(twice);
		I8 im = _mm256_set1_epi32(w[1]), iq = _mm256_set1_epi32(q[1]);
		for (int j = 0, k = 0; j < length; j += 4 * stride, ++k) {
			I8 r = _mm256_set1_epi32(w[k]), rq = _mm256_set1_epi32(q[k]);
			I8 r0 = _mm256_set1_epi32(w[2 * k]), q0 = _mm256_set1_epi32(q[2 * k]);
			I8 r3 = _mm256_set1_epi32(w3[k]), q3 = _mm256_set1_epi32(q3s[k]);
			for (int i = 0; i < stride; i += 8) {
				I8 x0 = load8(a + j + i), x1 = load8(a + j + stride + i);
				I8 x2 = load8(a + j + 2 * stride + i), x3 = load8(a + j + 3 * stride + i);
				if constexpr (inverse) {
					I8 u = add8(x0, x1), v = sub8(x0, x1);
					I8 s = add8(x2, x3);
					I8 t = fixed_mul8(_mm256_sub_epi32(_mm256_add_epi32(x2, p), x3), im, iq);
					x0 = add8(u, s);
					x1 = _mm256_add_epi32(v, t);
					x2 = _mm256_sub_epi32(_mm256_add_epi32(u, p), s);
					x3 = _mm256_sub_epi32(_mm256_add_epi32(v, p), t);
					if (k) {
						x1 = fixed_mul8(x1, r0, q0);
						x2 = fixed_mul8(x2, r, rq);
						x3 = fixed_mul8(x3, r3, q3);
					} else {
						x1 = _mm256_min_epu32(x1, _mm256_sub_epi32(x1, p));
						x2 = _mm256_min_epu32(x2, _mm256_sub_epi32(x2, p));
						x3 = _mm256_min_epu32(x3, _mm256_sub_epi32(x3, p));
					}
				} else {
					x0 = _mm256_min_epu32(x0, _mm256_sub_epi32(x0, p));
					if (!k) {
						x1 = _mm256_min_epu32(x1, _mm256_sub_epi32(x1, p));
						x2 = _mm256_min_epu32(x2, _mm256_sub_epi32(x2, p));
						x3 = _mm256_min_epu32(x3, _mm256_sub_epi32(x3, p));
					}
					if (k) {
						x1 = fixed_mul8(x1, r0, q0);
						x2 = fixed_mul8(x2, r, rq);
						x3 = fixed_mul8(x3, r3, q3);
					}
					I8 u = add8(x0, x2), v = sub8(x0, x2), s = add8(x1, x3);
					I8 t = fixed_mul8(_mm256_sub_epi32(_mm256_add_epi32(x1, p), x3), im, iq);
					x0 = _mm256_add_epi32(u, s);
					x1 = _mm256_sub_epi32(_mm256_add_epi32(u, p), s);
					x2 = _mm256_add_epi32(v, t);
					x3 = _mm256_sub_epi32(_mm256_add_epi32(v, p), t);
				}
				store8(a + j + i, x0);
				store8(a + j + stride + i, x1);
				store8(a + j + 2 * stride + i, x2);
				store8(a + j + 3 * stride + i, x3);
			}
		}
	}
	template <bool inverse>
	__attribute__((target("avx2"))) void stage8(U *a, int stride, int length, int first = 0) const {
		const V &w = inverse ? inverse_roots : roots;
		const V &quot = inverse ? inverse_root_quotients : root_quotients;
		if (stride < 8) {
			for (int j = 0; j < length; j += 8) {
				I8 z = load8(a + j);
				if (stride == 1)
					z = small8<1, inverse>(z, j / 8);
				else if (stride == 2)
					z = small8<2, inverse>(z, j / 8);
				else
					z = small8<4, inverse>(z, j / 8);
				store8(a + j, z);
			}
		} else {
			for (int j = 0, k = first; j < length; j += 2 * stride, ++k) {
				I8 r = _mm256_set1_epi32(w[k]), q = _mm256_set1_epi32(quot[k]);
				for (int i = 0; i < stride; i += 8) {
					I8 x = load8(a + j + i), y = load8(a + j + stride + i);
					if constexpr (!inverse)
						if (k) y = fixed_mul8(y, r, q);
					I8 z = sub8(x, y);
					x = add8(x, y);
					if constexpr (inverse)
						if (k) z = fixed_mul8(z, r, q);
					store8(a + j + i, x);
					store8(a + j + stride + i, z);
				}
			}
		}
	}

	template <bool inverse, bool lazy = false>
	__attribute__((target("avx2"), always_inline)) static void fixed_butterfly8(I8 &x, I8 &y, I8 r, I8 q, bool unit = false) {
		if constexpr (lazy && !inverse) {
			I8 p = _mm256_set1_epi32(twice);
			x = _mm256_min_epu32(x, _mm256_sub_epi32(x, p));
			y = unit ? _mm256_min_epu32(y, _mm256_sub_epi32(y, p)) : fixed_mul8(y, r, q);
			I8 z = _mm256_sub_epi32(_mm256_add_epi32(x, p), y);
			x = _mm256_add_epi32(x, y);
			y = z;
		} else {
			if constexpr (!inverse)
				if (!unit) y = fixed_mul8(y, r, q);
			I8 z = sub8(x, y);
			x = add8(x, y);
			y = inverse && !unit ? fixed_mul8(z, r, q) : z;
		}
	}
	template <bool inverse, bool lazy = false>
	__attribute__((target("avx2"))) void stage_pair8(U *a, int stride, int length, int first = 0) const {
		const V &w = inverse ? inverse_roots : roots;
		const V &q = inverse ? inverse_root_quotients : root_quotients;
		for (int j = 0, k = first; j < length; j += 4 * stride, ++k) {
			I8 r = _mm256_set1_epi32(w[k]), rq = _mm256_set1_epi32(q[k]);
			I8 r0 = _mm256_set1_epi32(w[2 * k]), q0 = _mm256_set1_epi32(q[2 * k]);
			I8 r1 = _mm256_set1_epi32(w[2 * k + 1]), q1 = _mm256_set1_epi32(q[2 * k + 1]);
			for (int i = 0; i < stride; i += 8) {
				I8 x0 = load8(a + j + i), x1 = load8(a + j + stride + i);
				I8 x2 = load8(a + j + 2 * stride + i), x3 = load8(a + j + 3 * stride + i);
				if constexpr (inverse) {
					fixed_butterfly8<inverse, lazy>(x0, x1, r0, q0, k == 0);
					fixed_butterfly8<inverse, lazy>(x2, x3, r1, q1, false);
					fixed_butterfly8<inverse, lazy>(x0, x2, r, rq, k == 0);
					fixed_butterfly8<inverse, lazy>(x1, x3, r, rq, k == 0);
				} else {
					fixed_butterfly8<inverse, lazy>(x0, x2, r, rq, k == 0);
					fixed_butterfly8<inverse, lazy>(x1, x3, r, rq, k == 0);
					fixed_butterfly8<inverse, lazy>(x0, x1, r0, q0, k == 0);
					fixed_butterfly8<inverse, lazy>(x2, x3, r1, q1, false);
				}
				store8(a + j + i, x0);
				store8(a + j + stride + i, x1);
				store8(a + j + 2 * stride + i, x2);
				store8(a + j + 3 * stride + i, x3);
			}
		}
	}
	__attribute__((target("avx2"))) static void fold8(U *a, int columns) {
		for (int i = 0; i < columns; i += 8) {
			I8 x = add8(load8(a + i), load8(a + columns + i));
			store8(a + i, x);
			store8(a + columns + i, x);
		}
	}
	__attribute__((target("avx2"))) void project8_8(U *a, int rows) const {
		for (int j = 0; j < rows; ++j) {
			I8 x0 = load8(a + 16 * j + 0);
			x0 = small8<4, true>(small8<2, true>(small8<1, true>(x0, 0), 0), 0);
			I8 x1 = load8(a + 16 * j + 8);
			x1 = small8<4, true>(small8<2, true>(small8<1, true>(x1, 1), 1), 1);
			x0 = add8(x0, x1);
			x1 = x0;
			x0 = small8<1, false>(small8<2, false>(small8<4, false>(x0, 0), 0), 0);
			store8(a + 16 * j + 0, x0);
			x1 = small8<1, false>(small8<2, false>(small8<4, false>(x1, 1), 1), 1);
			store8(a + 16 * j + 8, x1);
		}
	}
	__attribute__((target("avx2"))) void project16_8(U *a, int rows) const {
		for (int j = 0; j < rows; ++j) {
			I8 x0 = load8(a + 32 * j + 0);
			x0 = small8<4, true>(small8<2, true>(small8<1, true>(x0, 0), 0), 0);
			I8 x1 = load8(a + 32 * j + 8);
			x1 = small8<4, true>(small8<2, true>(small8<1, true>(x1, 1), 1), 1);
			I8 x2 = load8(a + 32 * j + 16);
			x2 = small8<4, true>(small8<2, true>(small8<1, true>(x2, 2), 2), 2);
			I8 x3 = load8(a + 32 * j + 24);
			x3 = small8<4, true>(small8<2, true>(small8<1, true>(x3, 3), 3), 3);
			fixed_butterfly8<true>(x0, x1, _mm256_set1_epi32(inverse_roots[0]), _mm256_set1_epi32(inverse_root_quotients[0]), true);
			fixed_butterfly8<true>(x2, x3, _mm256_set1_epi32(inverse_roots[1]), _mm256_set1_epi32(inverse_root_quotients[1]), false);
			x0 = add8(x0, x2);
			x2 = x0;
			x1 = add8(x1, x3);
			x3 = x1;
			fixed_butterfly8<false>(x0, x1, _mm256_set1_epi32(roots[0]), _mm256_set1_epi32(root_quotients[0]), true);
			fixed_butterfly8<false>(x2, x3, _mm256_set1_epi32(roots[1]), _mm256_set1_epi32(root_quotients[1]), false);
			x0 = small8<1, false>(small8<2, false>(small8<4, false>(x0, 0), 0), 0);
			store8(a + 32 * j + 0, x0);
			x1 = small8<1, false>(small8<2, false>(small8<4, false>(x1, 1), 1), 1);
			store8(a + 32 * j + 8, x1);
			x2 = small8<1, false>(small8<2, false>(small8<4, false>(x2, 2), 2), 2);
			store8(a + 32 * j + 16, x2);
			x3 = small8<1, false>(small8<2, false>(small8<4, false>(x3, 3), 3), 3);
			store8(a + 32 * j + 24, x3);
		}
	}
	__attribute__((target("avx2"))) static void scale_rows8(U *a, int rows, int columns, const U *factors) {
		if (columns == 2) {
			int j = 0;
			for (; j + 1 < rows; j += 2) {
				U w0 = factors[j], w1 = factors[j + 1], q0 = quotient(w0), q1 = quotient(w1);
				I8 w = _mm256_setr_epi32(w0, w0, w0, w0, w1, w1, w1, w1);
				I8 q = _mm256_setr_epi32(q0, q0, q0, q0, q1, q1, q1, q1);
				store8(a + 4 * j, fixed_mul8(load8(a + 4 * j), w, q));
			}
			if (j < rows)
				for (int i = 0; i < 4; ++i) a[4 * j + i] = fixed_mul(a[4 * j + i], factors[j], quotient(factors[j]));
		} else {
			for (int j = 0; j < rows; ++j) {
				U w = factors[j];
				I8 v = _mm256_set1_epi32(w), q = _mm256_set1_epi32(quotient(w));
				for (int i = 0; i < 2 * columns; i += 8) store8(a + 2 * columns * j + i, fixed_mul8(load8(a + 2 * columns * j + i), v, q));
			}
		}
	}
	__attribute__((target("avx2"))) void project4_8(U *a, int rows) const {
		for (int j = 0; j < rows; ++j) {
			I8 x = small8<2, true>(small8<1, true>(load8(a + 8 * j), 0), 0);
			x = add8(x, _mm256_permute2x128_si256(x, x, 1));
			x = small8<1, false>(small8<2, false>(x, 0), 0);
			store8(a + 8 * j, x);
		}
	}
	__attribute__((target("avx2"))) static void load_pairs(const U *a, I8 &x, I8 &y) {
		I8 index = _mm256_setr_epi32(0, 2, 4, 6, 1, 3, 5, 7);
		I8 u = _mm256_permutevar8x32_epi32(load8(a), index),
		   v = _mm256_permutevar8x32_epi32(load8(a + 8), index);
		x = _mm256_permute2x128_si256(u, v, 0x20);
		y = _mm256_permute2x128_si256(u, v, 0x31);
	}
	__attribute__((target("avx2"))) static void store_pairs(U *a, I8 x, I8 y) {
		I8 index = _mm256_setr_epi32(0, 4, 1, 5, 2, 6, 3, 7);
		store8(a, _mm256_permutevar8x32_epi32(_mm256_permute2x128_si256(x, y, 0x20), index));
		store8(a + 8,
			   _mm256_permutevar8x32_epi32(_mm256_permute2x128_si256(x, y, 0x31), index));
	}
	__attribute__((target("avx2"))) static void scale_copy8(U *dst, const U *src, int length, U r) {
		U quotient_r = quotient(r);
		I8 v = _mm256_set1_epi32(r), q = _mm256_set1_epi32(quotient_r);
		int i = 0;
		for (; i + 8 <= length; i += 8) store8(dst + i, fixed_mul8(load8(src + i), v, q));
		for (; i < length; ++i) dst[i] = fixed_mul(src[i], r, quotient_r);
	}
	__attribute__((target("avx2"))) static void scale_normalize8(U *a, int n, U r) {
		U qr = quotient(r);
		I8 v = _mm256_set1_epi32(r), q = _mm256_set1_epi32(qr), p = _mm256_set1_epi32(mod);
		int i = 0;
		for (; i + 8 <= n; i += 8) {
			I8 x = fixed_mul8(load8(a + i), v, q);
			store8(a + i, _mm256_min_epu32(x, _mm256_sub_epi32(x, p)));
		}
		for (; i < n; ++i) {
			U x = fixed_mul(a[i], r, qr);
			a[i] = x < mod ? x : x - mod;
		}
	}
	__attribute__((target("avx2"))) static void square_fold8(U *a, const U *q, int n, bool negative) {
		for (int i = 0; i < n; i += 8) {
			I8 x = load8(a + i), y = load8(a + n + i);
			I8 u = mul8(load8(q + i), load8(q + i));
			I8 v = mul8(load8(q + n + i), load8(q + n + i));
			x = add8(x, x);
			store8(a + i, negative ? sub8(x, u) : add8(x, u));
			store8(a + n + i, add8(add8(y, y), v));
		}
	}
	__attribute__((target("avx2"))) static void scale8(U *a, int length, U r) {
		I8 v = _mm256_set1_epi32(r), q = _mm256_set1_epi32(quotient(r));
		int i = 0;
		for (; i + 8 <= length; i += 8) store8(a + i, fixed_mul8(load8(a + i), v, q));
		for (; i < length; ++i) a[i] = fixed_mul(a[i], r, quotient(r));
	}
	__attribute__((target("avx2"))) static void pair_product8(U *a, const U *b, int length) {
		I8 ni = _mm256_set1_epi32(998244351), modulus = _mm256_set1_epi32(mod), p = _mm256_set1_epi32(twice);
		I8 index = _mm256_setr_epi32(0, 2, 4, 6, 1, 3, 5, 7);
		for (int i = 0; i < length; i += 8) {
			I8 x = load8(b + 2 * i), y = load8(b + 2 * i + 8);
			x = _mm256_min_epu32(x, _mm256_sub_epi32(x, p));
			y = _mm256_min_epu32(y, _mm256_sub_epi32(y, p));
			x = _mm256_mul_epu32(x, _mm256_srli_epi64(x, 32));
			y = _mm256_mul_epu32(y, _mm256_srli_epi64(y, 32));
			x = _mm256_add_epi64(x, _mm256_mul_epu32(_mm256_mul_epu32(x, ni), modulus));
			y = _mm256_add_epi64(y, _mm256_mul_epu32(_mm256_mul_epu32(y, ni), modulus));
			x = _mm256_blend_epi32(_mm256_srli_epi64(x, 32), y, 0xaa);
			store8(a + i, _mm256_permutevar8x32_epi32(x, index));
		}
	}
	__attribute__((target("avx2"))) static void pair_multiply8(U *a, const U *b, int length) {
		I8 index0 = _mm256_setr_epi32(0, 0, 1, 1, 2, 2, 3, 3), index1 = _mm256_setr_epi32(4, 4, 5, 5, 6, 6, 7, 7);
		I8 p = _mm256_set1_epi32(twice);
		for (int i = 0; i < length; i += 8) {
			I8 r = load8(b + i), x = load8(a + 2 * i), y = load8(a + 2 * i + 8);
			x = _mm256_min_epu32(x, _mm256_sub_epi32(x, p));
			y = _mm256_min_epu32(y, _mm256_sub_epi32(y, p));
			x = _mm256_shuffle_epi32(x, 0xb1);
			y = _mm256_shuffle_epi32(y, 0xb1);
			store8(a + 2 * i, mul8(x, _mm256_permutevar8x32_epi32(r, index0)));
			store8(a + 2 * i + 8, mul8(y, _mm256_permutevar8x32_epi32(r, index1)));
		}
	}
	__attribute__((target("avx2"))) static void difference8(U *a, const U *b, const U *c,
															int length) {
		for (int i = 0; i < length; i += 8) store8(a + i, sub8(load8(b + i), load8(c + i)));
	}
	__attribute__((target("avx2"))) static void project1_8(U *a, int rows) {
		for (int j = 0; j < rows; j += 8) {
			I8 x, y;
			load_pairs(a + 2 * j, x, y);
			x = add8(x, y);
			store_pairs(a + 2 * j, x, x);
		}
	}
	__attribute__((target("avx2"))) void project2_8(U *a, int rows) const {
		const I8 index = _mm256_setr_epi32(0, 4, 1, 5, 2, 6, 3, 7);
		for (int j = 0; j < rows; j += 8) {
			U *p = a + 4 * j;
			I8 u = load8(p), v = load8(p + 8), w = load8(p + 16), z = load8(p + 24);
			I8 uv0 = _mm256_unpacklo_epi32(u, v), uv1 = _mm256_unpackhi_epi32(u, v);
			I8 wz0 = _mm256_unpacklo_epi32(w, z), wz1 = _mm256_unpackhi_epi32(w, z);
			I8 a0 = _mm256_permutevar8x32_epi32(_mm256_unpacklo_epi64(uv0, wz0), index);
			I8 a1 = _mm256_permutevar8x32_epi32(_mm256_unpackhi_epi64(uv0, wz0), index);
			I8 a2 = _mm256_permutevar8x32_epi32(_mm256_unpacklo_epi64(uv1, wz1), index);
			I8 a3 = _mm256_permutevar8x32_epi32(_mm256_unpackhi_epi64(uv1, wz1), index);
			I8 x = add8(add8(a0, a1), add8(a2, a3));
			I8 y =
				add8(sub8(a0, a1), fixed_mul8(sub8(a2, a3), _mm256_set1_epi32(inverse_roots[1]),
											  _mm256_set1_epi32(inverse_root_quotients[1])));
			a0 = add8(x, y);
			a1 = sub8(x, y);
			y = fixed_mul8(y, _mm256_set1_epi32(roots[1]),
						   _mm256_set1_epi32(root_quotients[1]));
			a2 = add8(x, y);
			a3 = sub8(x, y);
			uv0 = _mm256_unpacklo_epi32(a0, a1);
			uv1 = _mm256_unpackhi_epi32(a0, a1);
			wz0 = _mm256_unpacklo_epi32(a2, a3);
			wz1 = _mm256_unpackhi_epi32(a2, a3);
			u = _mm256_unpacklo_epi64(uv0, wz0);
			v = _mm256_unpackhi_epi64(uv0, wz0);
			w = _mm256_unpacklo_epi64(uv1, wz1);
			z = _mm256_unpackhi_epi64(uv1, wz1);
			store8(p, _mm256_permute2x128_si256(u, v, 0x20));
			store8(p + 8, _mm256_permute2x128_si256(w, z, 0x20));
			store8(p + 16, _mm256_permute2x128_si256(u, v, 0x31));
			store8(p + 24, _mm256_permute2x128_si256(w, z, 0x31));
		}
	}
#endif

#if defined(__ARM_NEON)
	template <int stride, bool inverse>
	void small4(I &a, I &b, int k) const {
		const V &w = inverse ? inverse_roots : roots;
		const V &quot = inverse ? inverse_root_quotients : root_quotients;
		if constexpr (stride == 4) {
			if (k == 0) {
				I x = sub4(a, b);
				a = add4(a, b);
				b = x;
			} else
				fixed_butterfly4<inverse>(a, b, vdupq_n_u32(w[k]), vdupq_n_u32(quot[k]));
		}
		if constexpr (stride == 2) {
			I x = vcombine_u32(vget_low_u32(a), vget_low_u32(b));
			I y = vcombine_u32(vget_high_u32(a), vget_high_u32(b));
			fixed_butterfly4<inverse>(
				x, y, vcombine_u32(vdup_n_u32(w[2 * k]), vdup_n_u32(w[2 * k + 1])),
				vcombine_u32(vdup_n_u32(quot[2 * k]), vdup_n_u32(quot[2 * k + 1])));
			a = vcombine_u32(vget_low_u32(x), vget_low_u32(y));
			b = vcombine_u32(vget_high_u32(x), vget_high_u32(y));
		}
		if constexpr (stride == 1) {
			uint32x4x2_t z = vuzpq_u32(a, b);
			fixed_butterfly4<inverse>(z.val[0], z.val[1], vld1q_u32(w.data() + 4 * k),
									  vld1q_u32(quot.data() + 4 * k));
			z = vzipq_u32(z.val[0], z.val[1]);
			a = z.val[0];
			b = z.val[1];
		}
	}
	template <bool inverse>
	void fft4(U *a, int length, int limit = 0) const {
		if (!limit) limit = length;
		if constexpr (!inverse) {
			int h = limit / 2;
			for (; h >= 16; h >>= 2) stage_pair<false>(a, h / 2, length);
			if (h == 8) stage<false>(a, h, length);
		}
		for (int j = 0; j < length; j += 8) {
			I x = vld1q_u32(a + j), y = vld1q_u32(a + j + 4);
			if constexpr (inverse) {
				small4<1, true>(x, y, j / 8);
				small4<2, true>(x, y, j / 8);
				small4<4, true>(x, y, j / 8);
			} else {
				small4<4, false>(x, y, j / 8);
				small4<2, false>(x, y, j / 8);
				small4<1, false>(x, y, j / 8);
			}
			vst1q_u32(a + j, x);
			vst1q_u32(a + j + 4, y);
		}
		if constexpr (inverse) {
			int h = 8;
			for (; 4 * h <= limit; h <<= 2) stage_pair<true>(a, h, length);
			if (h < limit) stage<true>(a, h, length);
		}
	}
#endif
	template <bool inverse>
	void stage(U *a, int stride, int length, int first = 0) const {
		const V &w = inverse ? inverse_roots : roots;
		const V &quot = inverse ? inverse_root_quotients : root_quotients;
#if defined(__x86_64__) || defined(__i386__)
		if (length >= 8 && vectorized()) {
			stage8<inverse>(a, stride, length, first);
			return;
		}
#endif
		for (int j = 0, k = first; j < length; j += 2 * stride, ++k) {
			U r = w[k], q = quot[k];
			int i = 0;
#if defined(__ARM_NEON)
			for (; i + 4 <= stride; i += 4) {
				I x = vld1q_u32(a + j + i), y = vld1q_u32(a + j + stride + i);
				if constexpr (!inverse)
					if (k) y = fixed_mul4(y, vdupq_n_u32(r), vdupq_n_u32(q));
				I z = sub4(x, y);
				x = add4(x, y);
				if constexpr (inverse)
					if (k) z = fixed_mul4(z, vdupq_n_u32(r), vdupq_n_u32(q));
				y = z;
				vst1q_u32(a + j + i, x);
				vst1q_u32(a + j + stride + i, y);
			}
#endif
			for (; i < stride; ++i) {
				U x = a[j + i], y = a[j + stride + i];
				if constexpr (!inverse)
					if (k) y = fixed_mul(y, r, q);
				a[j + i] = add(x, y);
				a[j + stride + i] = inverse && k ? fixed_mul(sub(x, y), r, q) : sub(x, y);
			}
		}
	}

	template <bool inverse>
	void stage_pair(U *a, int stride, int length, int first = 0) const {
#if defined(__x86_64__) || defined(__i386__)
		if (stride >= 8 && vectorized()) {
			stage_pair8<inverse>(a, stride, length, first);
			return;
		}
#endif
#if defined(__ARM_NEON)
		if (stride >= 4) {
			const V &w = inverse ? inverse_roots : roots;
			const V &q = inverse ? inverse_root_quotients : root_quotients;
			for (int j = 0, k = first; j < length; j += 4 * stride, ++k) {
				I r = vdupq_n_u32(w[k]), rq = vdupq_n_u32(q[k]);
				I r0 = vdupq_n_u32(w[2 * k]), q0 = vdupq_n_u32(q[2 * k]);
				I r1 = vdupq_n_u32(w[2 * k + 1]), q1 = vdupq_n_u32(q[2 * k + 1]);
				for (int i = 0; i < stride; i += 4) {
					I x0 = vld1q_u32(a + j + i), x1 = vld1q_u32(a + j + stride + i);
					I x2 = vld1q_u32(a + j + 2 * stride + i),
					  x3 = vld1q_u32(a + j + 3 * stride + i);
					if constexpr (inverse) {
						fixed_butterfly4<true>(x0, x1, r0, q0, k == 0);
						fixed_butterfly4<true>(x2, x3, r1, q1);
						fixed_butterfly4<true>(x0, x2, r, rq, k == 0);
						fixed_butterfly4<true>(x1, x3, r, rq, k == 0);
					} else {
						fixed_butterfly4<false>(x0, x2, r, rq, k == 0);
						fixed_butterfly4<false>(x1, x3, r, rq, k == 0);
						fixed_butterfly4<false>(x0, x1, r0, q0, k == 0);
						fixed_butterfly4<false>(x2, x3, r1, q1);
					}
					vst1q_u32(a + j + i, x0);
					vst1q_u32(a + j + stride + i, x1);
					vst1q_u32(a + j + 2 * stride + i, x2);
					vst1q_u32(a + j + 3 * stride + i, x3);
				}
			}
			return;
		}
#endif
		if constexpr (inverse)
			stage<true>(a, stride, length, 2 * first), stage<true>(a, 2 * stride, length, first);
		else
			stage<false>(a, 2 * stride, length, first), stage<false>(a, stride, length, 2 * first);
	}
	template <bool inverse>
	void fft(U *a, int length, int limit = 0) const {
		if (!limit) limit = length;
#if defined(__ARM_NEON)
		if (limit >= 8) {
			fft4<inverse>(a, length, limit);
			return;
		}
#elif defined(__x86_64__) || defined(__i386__)
		if (limit >= 8 && vectorized()) {
			fft8<inverse>(a, length, limit);
			return;
		}
#endif
		if constexpr (inverse) {
			for (int h = 1; h < limit; h <<= 1) stage<true>(a, h, length);
		} else {
			for (int h = limit / 2; h; h >>= 1) stage<false>(a, h, length);
		}
	}
#if defined(__x86_64__) || defined(__i386__)
	__attribute__((target("avx2"))) static void product8(U *dst, const U *a, const U *b, int n) {
		int i = 0;
		for (; i + 8 <= n; i += 8) {
			I8 x = mul8(load8(a + i), load8(b + i));
			store8(dst + i, x);
		}
		for (; i < n; ++i) dst[i] = mul(a[i], b[i]);
	}
	__attribute__((target("avx2"))) static void weighted8(U *dst, const U *src, int n, int first) {
		alignas(32) U factors[8];
		for (int j = 0; j < 8; ++j) factors[j] = in(first + j);
		I8 r = load8(factors), step = _mm256_set1_epi32(in(8));
		int i = 0;
		for (; i + 8 <= n; i += 8) {
			store8(dst + i, mul8(load8(src + i), r));
			r = add8(r, step);
		}
		for (; i < n; ++i) dst[i] = mul(src[i], in(first + i));
	}
	__attribute__((target("avx2"), always_inline)) static I8 inverse8(I8 a) {
		I8 a2 = mul8(a, a), a3 = mul8(a2, a);
		I8 a7 = mul8(mul8(a3, a3), a), a31 = a7;
		for (int repeat_index = 0; repeat_index < 2; ++repeat_index) a31 = mul8(a31, a31);
		a31 = mul8(a31, a3);
		I8 r = a7;
		for (int repeat_index = 0; repeat_index < 3; ++repeat_index) r = mul8(r, r);
		r = mul8(r, a3);
		r = mul8(r, r);
		for (int j = 0; j < 4; ++j) {
			for (int repeat_index = 0; repeat_index < 5; ++repeat_index) r = mul8(r, r);
			r = mul8(r, a31);
		}
		for (int repeat_index = 0; repeat_index < 3; ++repeat_index) r = mul8(r, r);
		return mul8(r, a7);
	}
	__attribute__((target("avx2"))) static void batch_inverse8(U *dst, const U *src, int n) {
		if (n < 16) {
			for (int i = 0; i < n; ++i) dst[i] = T::word(src[i]).inv().v;
			return;
		}
		int count = n / 16 * 16;
		I8 p = load8(src), q = load8(src + 8);
		store8(dst, p);
		store8(dst + 8, q);
		for (int i = 16; i < count; i += 16) {
			p = mul8(p, load8(src + i));
			q = mul8(q, load8(src + i + 8));
			store8(dst + i, p);
			store8(dst + i + 8, q);
		}
		alignas(32) U lanes[16];
		store8(lanes, p);
		store8(lanes + 8, q);
		for (auto x : lanes) assert(out(x) != 0);
		I8 inverse = inverse8(p), other = inverse8(q);
		for (int i = count - 16; i >= 16; i -= 16) {
			store8(dst + i, mul8(inverse, load8(dst + i - 16)));
			store8(dst + i + 8, mul8(other, load8(dst + i - 8)));
			inverse = mul8(inverse, load8(src + i));
			other = mul8(other, load8(src + i + 8));
		}
		store8(dst, inverse);
		store8(dst + 8, other);
		for (int i = count; i < n; ++i) dst[i] = T::word(src[i]).inv().v;
	}
	__attribute__((target("avx2"))) static void evaluate_leaf8(U *dst, const U *cur, const U *tree, int degree, const U *points, int count) {
		I8 p = _mm256_set1_epi32(mod), pp = _mm256_set1_epi32(twice);
		I8 x = load8(cur);
		x = _mm256_min_epu32(x, _mm256_sub_epi32(x, pp));
		x = _mm256_min_epu32(x, _mm256_sub_epi32(x, p));
		alignas(32) U window[16] = {}, coefficients[8] = {}, polynomial[8], tail[8] = {};
		store8(window + 8, x);
		std::copy_n(tree, std::min(degree, 8), coefficients);
		I8 y = load8(coefficients);
		y = _mm256_min_epu32(y, _mm256_sub_epi32(y, pp));
		store8(coefficients, _mm256_min_epu32(y, _mm256_sub_epi32(y, p)));
		I8 e = _mm256_setzero_si256(), o = e;
		for (int i = 0; i < 8; ++i) {
			I8 z = load8(window + 8 - i), r = _mm256_set1_epi32(coefficients[i]);
			e = _mm256_add_epi64(e, _mm256_mul_epu32(z, r));
			o = _mm256_add_epi64(o, _mm256_mul_epu32(_mm256_srli_epi64(z, 32), r));
		}
		I8 ni = _mm256_set1_epi32(998244351);
		e = _mm256_add_epi64(e, _mm256_mul_epu32(_mm256_mul_epu32(e, ni), p));
		o = _mm256_add_epi64(o, _mm256_mul_epu32(_mm256_mul_epu32(o, ni), p));
		x = _mm256_blend_epi32(_mm256_srli_epi64(e, 32), o, 0xaa);
		x = _mm256_min_epu32(x, _mm256_sub_epi32(x, pp));
		store8(polynomial, _mm256_min_epu32(x, _mm256_sub_epi32(x, p)));
		if (count == 8) x = load8(points);
		else {
			std::copy_n(points, count, tail);
			x = load8(tail);
		}
		I8 value = _mm256_setzero_si256();
		for (int i = 0; i < 8; ++i) value = add8(mul8(value, x), _mm256_set1_epi32(polynomial[i]));
		if (count == 8) store8(dst, value);
		else {
			store8(tail, value);
			std::copy_n(tail, count, dst);
		}
	}
	__attribute__((target("avx2"), always_inline)) static U sum_product8(I8 a, I8 b) {
		I8 p = _mm256_set1_epi32(mod);
		a = _mm256_min_epu32(a, _mm256_sub_epi32(a, p));
		b = _mm256_min_epu32(b, _mm256_sub_epi32(b, p));
		I8 sum = _mm256_add_epi64(_mm256_mul_epu32(a, b), _mm256_mul_epu32(_mm256_srli_epi64(a, 32), _mm256_srli_epi64(b, 32)));
		sum = _mm256_add_epi64(sum, _mm256_permute2x128_si256(sum, sum, 1));
		sum = _mm256_add_epi64(sum, _mm256_shuffle_epi32(sum, 0x4e));
		U result = reduce((unsigned long long)_mm_cvtsi128_si64(_mm256_castsi256_si128(sum)));
		return result < twice ? result : result - twice;
	}
	__attribute__((target("avx2"))) static void interpolate_leaf8(U *dst, const U *q, const U *points, const U *weights) {
		I8 x = load8(points), w = load8(weights), c = _mm256_set1_epi32(one);
		I8 sum = add8(w, _mm256_permute2x128_si256(w, w, 1));
		sum = add8(sum, _mm256_shuffle_epi32(sum, 0x4e));
		sum = add8(sum, _mm256_shuffle_epi32(sum, 0xb1));
		dst[7] = _mm_cvtsi128_si32(_mm256_castsi256_si128(sum));
		for (int k = 7 - 1; k >= 0; --k) {
			c = add8(mul8(c, x), _mm256_set1_epi32(q[k + 1]));
			dst[k] = sum_product8(c, w);
		}
	}
	__attribute__((target("avx2"))) static void evaluate_middle8(U *left, U *right, const U *cur, const U *ql, int dl, const U *qr, int dr) {
		I8 p = _mm256_set1_epi32(mod), pp = _mm256_set1_epi32(twice), ni = _mm256_set1_epi32(998244351);
		alignas(32) U window[16], coefficients[2][9] = {};
		for (int i = 0; i < 16; i += 8) {
			I8 x = load8(cur + i);
			x = _mm256_min_epu32(x, _mm256_sub_epi32(x, pp));
			store8(window + i, _mm256_min_epu32(x, _mm256_sub_epi32(x, p)));
		}
		for (int i = 0; i < std::min(dl, 9); ++i) {
			U x = ql[i];
			if (x >= twice) x -= twice;
			if (x >= mod) x -= mod;
			coefficients[0][i] = x;
		}
		for (int i = 0; i < std::min(dr, 9); ++i) {
			U x = qr[i];
			if (x >= twice) x -= twice;
			if (x >= mod) x -= mod;
			coefficients[1][i] = x;
		}
		I8 e[2] = {}, o[2] = {};
		for (int i = 0; i < 9; ++i) {
			I8 z = load8(window + 8 - i), high = _mm256_srli_epi64(z, 32);
			for (int k = 0; k < 2; ++k) {
				I8 r = _mm256_set1_epi32(coefficients[k][i]);
				e[k] = _mm256_add_epi64(e[k], _mm256_mul_epu32(z, r));
				o[k] = _mm256_add_epi64(o[k], _mm256_mul_epu32(high, r));
			}
		}
		U *dst[2] = {left, right};
		for (int k = 0; k < 2; ++k) {
			e[k] = _mm256_add_epi64(e[k], _mm256_mul_epu32(_mm256_mul_epu32(e[k], ni), p));
			o[k] = _mm256_add_epi64(o[k], _mm256_mul_epu32(_mm256_mul_epu32(o[k], ni), p));
			I8 x = _mm256_blend_epi32(_mm256_srli_epi64(e[k], 32), o[k], 0xaa);
			store8(dst[k], _mm256_min_epu32(x, _mm256_sub_epi32(x, pp)));
		}
	}
	__attribute__((target("avx2"))) static void inverse_sequence8(U *dst, int first, int n) {

		alignas(32) U initial[8];
		for (int j = 0; j < 8; ++j) initial[j] = in(first + j);
		I8 factor = load8(initial), step = _mm256_set1_epi32(in(8));
		I8 p = factor;
		store8(dst, p);
		for (int i = 8; i < n; i += 8) {
			factor = add8(factor, step);
			p = mul8(p, factor);
			store8(dst + i, p);
		}
		I8 inverse = _mm256_set1_epi32(one);
		for (U e = mod - 2; e; e >>= 1, p = mul8(p, p))
			if (e & 1) inverse = mul8(inverse, p);
		for (int i = n - 8; i >= 8; i -= 8) {
			store8(dst + i, mul8(inverse, load8(dst + i - 8)));
			inverse = mul8(inverse, factor);
			factor = sub8(factor, step);
		}
		store8(dst, inverse);
	}
	__attribute__((target("avx2"))) static void normalize8(U *a, int n) {
		I8 p = _mm256_set1_epi32(mod);
		int i = 0;
		for (; i + 8 <= n; i += 8) {
			I8 x = load8(a + i);
			store8(a + i, _mm256_min_epu32(x, _mm256_sub_epi32(x, p)));
		}
		for (; i < n; ++i)
			if (a[i] >= mod) a[i] -= mod;
	}
#endif
	static void weighted(U *dst, const U *src, int n, int first) {
#if defined(__x86_64__) || defined(__i386__)
		if (n >= 8 && vectorized()) {
			weighted8(dst, src, n, first);
			return;
		}
#endif
		for (int i = 0; i < n; ++i) dst[i] = mul(src[i], in(first + i));
	}
	static void scale_normalize(U *a, int n, U r) {
#if defined(__x86_64__) || defined(__i386__)
		if (n >= 8 && vectorized()) {
			scale_normalize8(a, n, r);
			return;
		}
#endif
		U qr = quotient(r);
		int i = 0;
#if defined(__ARM_NEON)
		I v = vdupq_n_u32(r), q = vdupq_n_u32(qr), p = vdupq_n_u32(mod);
		for (; i + 4 <= n; i += 4) {
			I x = fixed_mul4(vld1q_u32(a + i), v, q);
			vst1q_u32(a + i, vminq_u32(x, vsubq_u32(x, p)));
		}
#endif
		for (; i < n; ++i) {
			U x = fixed_mul(a[i], r, qr);
			a[i] = x < mod ? x : x - mod;
		}
	}
	static void square_fold(U *a, const U *q, int n, bool negative) {
#if defined(__x86_64__) || defined(__i386__)
		if (n >= 8 && vectorized()) {
			square_fold8(a, q, n, negative);
			return;
		}
#endif
		int i = 0;
#if defined(__ARM_NEON)
		for (; i + 4 <= n; i += 4) {
			I x = vld1q_u32(a + i), y = vld1q_u32(a + n + i);
			I u = mul4(vld1q_u32(q + i), vld1q_u32(q + i));
			I v = mul4(vld1q_u32(q + n + i), vld1q_u32(q + n + i));
			x = add4(x, x);
			vst1q_u32(a + i, negative ? sub4(x, u) : add4(x, u));
			vst1q_u32(a + n + i, add4(add4(y, y), v));
		}
#endif
		for (; i < n; ++i) {
			U u = mul(q[i], q[i]), v = mul(q[n + i], q[n + i]);
			U x = add(a[i], a[i]);
			a[i] = negative ? sub(x, u) : add(x, u);
			a[n + i] = add(add(a[n + i], a[n + i]), v);
		}
	}
	static void normalize(U *a, int n) {
#if defined(__x86_64__) || defined(__i386__)
		if (n >= 8 && vectorized()) {
			normalize8(a, n);
			return;
		}
#endif
		int i = 0;
#if defined(__ARM_NEON)
		for (; i + 4 <= n; i += 4) {
			I x = vld1q_u32(a + i);
			vst1q_u32(a + i, vminq_u32(x, vsubq_u32(x, vdupq_n_u32(mod))));
		}
#endif
		for (; i < n; ++i)
			if (a[i] >= mod) a[i] -= mod;
	}
	static void product(U *dst, const U *a, const U *b, int n) {
#if defined(__x86_64__) || defined(__i386__)
		if (n >= 8 && vectorized()) {
			product8(dst, a, b, n);
			return;
		}
#endif
		int i = 0;
#if defined(__ARM_NEON)
		for (; i + 4 <= n; i += 4) {
			I x = mul4(vld1q_u32(a + i), vld1q_u32(b + i));
			vst1q_u32(dst + i, x);
		}
#endif
		for (; i < n; ++i) dst[i] = mul(a[i], b[i]);
	}
#if defined(__x86_64__) || defined(__i386__)
	template <bool inverse>
	__attribute__((target("avx2"), always_inline)) void convolution_stage8(U *a, int stride, int k) const {
		const V &w = inverse ? inverse_roots : roots;
		const V &q = inverse ? inverse_root_quotients : root_quotients;
		I8 r = _mm256_set1_epi32(w[k]), rq = _mm256_set1_epi32(q[k]);
		for (int i = 0; i < stride; i += 8) {
			I8 x = load8(a + i), y = load8(a + stride + i);
			fixed_butterfly8<inverse, true>(x, y, r, rq, k == 0);
			store8(a + i, x);
			store8(a + stride + i, y);
		}
	}
	__attribute__((target("avx2"), always_inline)) void leaf8(U *a, const U *b, int k) const {
		U w = roots[k / 2], q = root_quotients[k / 2];
		if (k & 1) w = mod - w, q = ~q;
		I8 p = _mm256_set1_epi32(mod), pp = _mm256_set1_epi32(twice);
		I8 x = load8(a), y = load8(b);
		x = _mm256_min_epu32(x, _mm256_sub_epi32(x, pp));
		y = _mm256_min_epu32(y, _mm256_sub_epi32(y, pp));
		x = _mm256_min_epu32(x, _mm256_sub_epi32(x, p));
		y = _mm256_min_epu32(y, _mm256_sub_epi32(y, p));
		alignas(32) U window[16], coefficients[8];
		I8 low = fixed_mul8(x, _mm256_set1_epi32(w), _mm256_set1_epi32(q));
		store8(window, _mm256_min_epu32(low, _mm256_sub_epi32(low, p)));
		store8(window + 8, x);
		store8(coefficients, y);
		I8 e = _mm256_setzero_si256(), o = e;
		for (int i = 0; i < 8; ++i) {
			I8 z = load8(window + 8 - i), r = _mm256_set1_epi32(coefficients[i]);
			e = _mm256_add_epi64(e, _mm256_mul_epu32(z, r));
			o = _mm256_add_epi64(o, _mm256_mul_epu32(_mm256_srli_epi64(z, 32), r));
		}
		I8 ni = _mm256_set1_epi32(998244351);
		e = _mm256_add_epi64(e, _mm256_mul_epu32(_mm256_mul_epu32(e, ni), p));
		o = _mm256_add_epi64(o, _mm256_mul_epu32(_mm256_mul_epu32(o, ni), p));
		x = _mm256_blend_epi32(_mm256_srli_epi64(e, 32), o, 0xaa);
		x = _mm256_min_epu32(x, _mm256_sub_epi32(x, pp));
		store8(a, _mm256_min_epu32(x, _mm256_sub_epi32(x, p)));
	}
	template <int count, bool second = false>
	__attribute__((target("avx2"), always_inline)) void inverse_leaf_products8(U *a, U *b, const U *ws, const U *qs) const {
		I8 p = _mm256_set1_epi32(mod), pp = _mm256_set1_epi32(twice);
		alignas(32) U window[count][16], coefficients[count][8];
		I8 e[count], o[count];
		for (int k = 0; k < count; ++k) {
			I8 x = load8(a + 8 * k), y = load8(b + 8 * k);
			if constexpr (!second) x = _mm256_min_epu32(x, _mm256_sub_epi32(x, pp));
			if constexpr (!second) y = _mm256_min_epu32(y, _mm256_sub_epi32(y, pp));
			x = _mm256_min_epu32(x, _mm256_sub_epi32(x, p));
			if constexpr (!second) y = _mm256_min_epu32(y, _mm256_sub_epi32(y, p));
			I8 low = fixed_mul8(x, _mm256_set1_epi32(ws[k]), _mm256_set1_epi32(qs[k]));
			store8(window[k], _mm256_min_epu32(low, _mm256_sub_epi32(low, p)));
			store8(window[k] + 8, x);
			store8(coefficients[k], y);
			if constexpr (!second) store8(b + 8 * k, y);
			e[k] = o[k] = _mm256_setzero_si256();
		}
		for (int i = 0; i < 8; ++i)
			for (int k = 0; k < count; ++k) {
				I8 z = load8(window[k] + 8 - i), r = _mm256_set1_epi32(coefficients[k][i]);
				e[k] = _mm256_add_epi64(e[k], _mm256_mul_epu32(z, r));
				o[k] = _mm256_add_epi64(o[k], _mm256_mul_epu32(_mm256_srli_epi64(z, 32), r));
			}
		I8 ni = _mm256_set1_epi32(998244351);
		for (int k = 0; k < count; ++k) {
			I8 x = e[k], y = o[k];
			x = _mm256_add_epi64(x, _mm256_mul_epu32(_mm256_mul_epu32(x, ni), p));
			y = _mm256_add_epi64(y, _mm256_mul_epu32(_mm256_mul_epu32(y, ni), p));
			x = _mm256_blend_epi32(_mm256_srli_epi64(x, 32), y, 0xaa);
			store8(a + 8 * k, _mm256_min_epu32(x, _mm256_sub_epi32(x, pp)));
		}
	}
	__attribute__((target("avx2"))) static U reverse_dot8(const U *a, const U *b, int n) {
		I8 p = _mm256_set1_epi32(mod), pp = _mm256_set1_epi32(twice), ni = _mm256_set1_epi32(998244351);
		I8 reversed = _mm256_setr_epi32(7, 6, 5, 4, 3, 2, 1, 0), total = _mm256_setzero_si256();
		int i = 0;
		for (; i + 8 <= n;) {
			I8 even = _mm256_setzero_si256(), odd = even;
			int end = std::min(n / 8 * 8, i + 64);
			for (; i < end; i += 8) {
				I8 x = load8(a + i), y = _mm256_permutevar8x32_epi32(load8(b - i - 7), reversed);
				x = _mm256_min_epu32(x, _mm256_sub_epi32(x, pp));
				y = _mm256_min_epu32(y, _mm256_sub_epi32(y, pp));
				x = _mm256_min_epu32(x, _mm256_sub_epi32(x, p));
				y = _mm256_min_epu32(y, _mm256_sub_epi32(y, p));
				even = _mm256_add_epi64(even, _mm256_mul_epu32(x, y));
				odd = _mm256_add_epi64(odd, _mm256_mul_epu32(_mm256_srli_epi64(x, 32), _mm256_srli_epi64(y, 32)));
			}
			even = _mm256_add_epi64(even, _mm256_mul_epu32(_mm256_mul_epu32(even, ni), p));
			odd = _mm256_add_epi64(odd, _mm256_mul_epu32(_mm256_mul_epu32(odd, ni), p));
			I8 value = _mm256_blend_epi32(_mm256_srli_epi64(even, 32), odd, 0xaa);
			value = _mm256_min_epu32(value, _mm256_sub_epi32(value, pp));
			total = add8(total, value);
		}
		total = add8(total, _mm256_permute2x128_si256(total, total, 1));
		total = add8(total, _mm256_shuffle_epi32(total, 0x4e));
		total = add8(total, _mm256_shuffle_epi32(total, 0xb1));
		U result = _mm_cvtsi128_si32(_mm256_castsi256_si128(total));
		for (; i < n; ++i) {
			U x = a[i], y = b[-i];
			if (x >= twice) x -= twice;
			if (y >= twice) y -= twice;
			if (x >= mod) x -= mod;
			if (y >= mod) y -= mod;
			result = add(result, mul(x, y));
		}
		return result;
	}

	__attribute__((target("avx2"))) static void inverse_fold8(U *a, const U *f, int n) {
		int i = 0;
		for (; i + 8 <= n; i += 8) store8(a + i, add8(load8(a + i), load8(f + i)));
		for (; i < n; ++i) {
			a[i] += f[i];
			if (a[i] >= twice) a[i] -= twice;
		}
	}
	__attribute__((target("avx2"))) static void inverse_recover8(U *g, const U *a, int m) {
		int n = 2 * m;
		U ratio = T(3).pow((mod - 1) / (2 * (n / 8))).inv().val();
		U first = mod - (mod - 1) / (2 * (n / 8));
		U phase = T(3).pow((mod - 1) / 4).inv().val();
		U powers[8];
		powers[0] = first;
		for (int j = 1; j < 8; ++j) powers[j] = (unsigned long long)powers[j - 1] * ratio % mod;
		U step = T(ratio).pow(8).val();
		I8 factors = load8(powers), s = _mm256_set1_epi32(step), sq = _mm256_set1_epi32(quotient(step));
		I8 vphase = _mm256_set1_epi32(phase), qphase = _mm256_set1_epi32(quotient(phase)), p = _mm256_set1_epi32(mod), z = _mm256_setzero_si256();
		for (int i = 0; i < m; i += 64) {
			store8(powers, factors);
			for (int j = 0; j < 8; ++j) {
				int at = i + 8 * j;
				U w = powers[j];
				I8 v = _mm256_set1_epi32(w), q = _mm256_set1_epi32(quotient(w));
				I8 r = _mm256_add_epi32(load8(a + at), fixed_mul8(load8(a + m + at), vphase, qphase));
				r = fixed_mul8(r, v, q);
				store8(g + m + at, sub8(z, add8(r, load8(g + m + at))));
			}
			factors = fixed_mul8(factors, s, sq);
			factors = _mm256_min_epu32(factors, _mm256_sub_epi32(factors, p));
		}
	}
	template <bool inverse>
	__attribute__((target("avx2"))) void inverse_partial8(U *a, int n, int limit = 0) const {
		if (!limit) limit = n;
		if constexpr (!inverse) {
			int stride = limit / 2;
			for (; stride >= 16; stride >>= 2) radix4_8<false>(a, stride / 2, n);
			if (stride == 8)
				for (int i = 0; i < n; i += 16) convolution_stage8<false>(a + i, 8, i / 16);
		} else {
			int stride = 8;
			for (; 4 * stride <= n; stride <<= 2) radix4_8<true>(a, stride, n);
			if (stride < n) stage8<true>(a, stride, n);
		}
	}
	__attribute__((target("avx2"))) static void inverse_twist8(U *a, U *b, int n, U ratio, U first = one) {
		U powers[8];
		powers[0] = out(first);
		U r = out(ratio);
		for (int j = 1; j < 8; ++j) powers[j] = (unsigned long long)powers[j - 1] * r % mod;
		U step = T::word(ratio).pow(8).val();
		I8 factors = load8(powers), s = _mm256_set1_epi32(step), sq = _mm256_set1_epi32(quotient(step)), p = _mm256_set1_epi32(mod);
		for (int i = 0; i < n; i += 64) {
			store8(powers, factors);
			for (int j = 0; j < 8 && i + 8 * j < n; ++j) {
				U w = powers[j];
				I8 v = _mm256_set1_epi32(w), q = _mm256_set1_epi32(quotient(w));
				store8(a + i + 8 * j, fixed_mul8(load8(a + i + 8 * j), v, q));
				if (b && i + 8 * j < n / 2) store8(b + i + 8 * j, fixed_mul8(load8(b + i + 8 * j), v, q));
			}
			factors = fixed_mul8(factors, s, sq);
			factors = _mm256_min_epu32(factors, _mm256_sub_epi32(factors, p));
		}
	}
	template <bool negative>
	__attribute__((target("avx2"))) void inverse_cubic8(U *a, U *b, int n) const {
		U shift = one;
		if constexpr (negative) {
			T psi = T(3).pow((mod - 1) / (2 * (n / 8)));
			shift = psi.v;
			if (shift >= mod) shift -= mod;
			inverse_twist8(a, b, n, shift);
		}
		inverse_partial8<false>(a, n);
		if constexpr (negative) {
			std::copy_n(b, n / 2, b + n / 2);
			inverse_partial8<false>(b, n, n / 2);
		} else inverse_partial8<false>(b, n);
		for (int i = 0; i < n; i += 8 * 8) {
			U ws[8], qs[8];
			for (int j = 0; j < 8; ++j) {
				int k = i / 8 + j;
				U w = roots[k / 2];
				if (k & 1) w = mod - w;
				if constexpr (negative) {
					w = mul(w, shift);
					if (w >= mod) w -= mod;
					qs[j] = quotient(w);
				} else qs[j] = (k & 1) ? ~root_quotients[k / 2] : root_quotients[k / 2];
				ws[j] = w;
			}
			inverse_leaf_products8<8>(a + i, b + i, ws, qs);
			inverse_leaf_products8<8, true>(a + i, b + i, ws, qs);
			for (int j = 0; j < 8; ++j) store8(a + i + 8 * j, sub8(load8(a + i + 8 * j), load8(b + i + 8 * j)));
		}
		inverse_partial8<true>(a, n);
		U normal = mod - (mod - 1) / (2 * (n / 8));
		if constexpr (!negative) scale8(a, n, normal);
	}
	template <int n>
	__attribute__((target("avx2"), always_inline)) void convolution_tile8(U *a, U *b, int k) const {
		if constexpr (n == 8) leaf8(a, b, k);
		else if constexpr (n == 16) {
			convolution_stage8<false>(a, n / 2, k);
			convolution_stage8<false>(b, n / 2, k);
			convolution_tile8<n / 2>(a, b, 2 * k);
			convolution_tile8<n / 2>(a + n / 2, b + n / 2, 2 * k + 1);
			convolution_stage8<true>(a, n / 2, k);
		} else {
			stage_pair8<false, true>(a, n / 4, n, k);
			stage_pair8<false, true>(b, n / 4, n, k);
			for (int i = 0; i < 4; ++i) convolution_tile8<n / 4>(a + i * (n / 4), b + i * (n / 4), 4 * k + i);
			stage_pair8<true>(a, n / 4, n, k);
		}
	}
	__attribute__((target("avx2"))) void convolution8(U *a, U *b, int n, int k = 0) const {
		if (n == 256) {
			convolution_tile8<256>(a, b, k);
			return;
		}
		if (n == 128) {
			convolution_tile8<128>(a, b, k);
			return;
		}
		if (n == 64) {
			convolution_tile8<64>(a, b, k);
			return;
		}
		stage_pair8<false, true>(a, n / 4, n, k);
		stage_pair8<false, true>(b, n / 4, n, k);
		for (int i = 0; i < 4; ++i) convolution8(a + i * (n / 4), b + i * (n / 4), n / 4, 4 * k + i);
		stage_pair8<true>(a, n / 4, n, k);
	}
#endif
#if defined(__x86_64__) || defined(__i386__)
	template <bool odd>
	__attribute__((target("avx2"))) static void reduce_pairs8(U *a, U *b, const U *twist, int n) {
		for (int i = 0; i < n; i += 8) {
			I8 x, y, u, v;
			load_pairs(a + 2 * i, x, y);
			load_pairs(b + 2 * i, u, v);
			x = mul8(x, v);
			y = mul8(y, u);
			if constexpr (odd) x = mul8(sub8(x, y), load8(twist + i));
			else x = add8(x, y);
			store8(a + i, x);
			store8(b + i, mul8(u, v));
		}
	}
	__attribute__((target("avx2"))) static void combine8(U *a, const U *b, const U *x, const U *y, int n) {
		for (int i = 0; i < n; i += 8) store8(a + i, add8(mul8(load8(a + i), load8(y + i)), mul8(load8(b + i), load8(x + i))));
	}
#endif

	template <bool odd>
	static void reduce_pairs(U *a, U *b, const U *twist, int n) {
#if defined(__x86_64__) || defined(__i386__)
		if (n >= 8 && vectorized()) {
			reduce_pairs8<odd>(a, b, twist, n);
			return;
		}
#endif
		int i = 0;
#if defined(__ARM_NEON)
		for (; i + 4 <= n; i += 4) {
			auto p = vld2q_u32(a + 2 * i), q = vld2q_u32(b + 2 * i);
			I x = mul4(p.val[0], q.val[1]), y = mul4(p.val[1], q.val[0]);
			if constexpr (odd) x = mul4(sub4(x, y), vld1q_u32(twist + i));
			else x = add4(x, y);
			vst1q_u32(a + i, x);
			vst1q_u32(b + i, mul4(q.val[0], q.val[1]));
		}
#endif
		for (; i < n; ++i) {
			U u = mul(a[2 * i], b[2 * i + 1]), v = mul(a[2 * i + 1], b[2 * i]);
			a[i] = odd ? mul(sub(u, v), twist[i]) : add(u, v);
			b[i] = mul(b[2 * i], b[2 * i + 1]);
		}
	}
	static void combine(U *a, const U *b, const U *x, const U *y, int n) {
#if defined(__x86_64__) || defined(__i386__)
		if (n >= 8 && vectorized()) {
			combine8(a, b, x, y, n);
			return;
		}
#endif
		int i = 0;
#if defined(__ARM_NEON)
		for (; i + 4 <= n; i += 4) vst1q_u32(a + i, add4(mul4(vld1q_u32(a + i), vld1q_u32(y + i)), mul4(vld1q_u32(b + i), vld1q_u32(x + i))));
#endif
		for (; i < n; ++i) a[i] = add(mul(a[i], y[i]), mul(b[i], x[i]));
	}
#if defined(__x86_64__) || defined(__i386__)
	__attribute__((target("avx2"))) static void geometric8(U *a, int n, U first, U step) {
		alignas(32) U factors[8];
		U power = one;
		for (int i = 0; i < 8; ++i) {
			factors[i] = mul(first, power);
			power = mul(power, step);
		}
		I8 r = load8(factors), ratio = _mm256_set1_epi32(power);
		int i = 0;
		for (; i + 8 <= n; i += 8) {
			store8(a + i, mul8(load8(a + i), r));
			r = mul8(r, ratio);
		}
		store8(factors, r);
		for (int j = 0; i < n; ++i, ++j) a[i] = mul(a[i], factors[j]);
	}
#endif
	static void geometric(U *a, int n, U first, U step) {
#if defined(__x86_64__) || defined(__i386__)
		if (n >= 8 && vectorized()) {
			geometric8(a, n, first, step);
			return;
		}
#endif
		int i = 0;
#if defined(__ARM_NEON)
		U factors[4], power = one;
		for (int j = 0; j < 4; ++j) {
			factors[j] = mul(first, power);
			power = mul(power, step);
		}
		I r = vld1q_u32(factors), ratio = vdupq_n_u32(power);
		for (; i + 4 <= n; i += 4) {
			vst1q_u32(a + i, mul4(vld1q_u32(a + i), r));
			r = mul4(r, ratio);
		}
		vst1q_u32(factors, r);
		first = factors[0];
#endif
		for (; i < n; ++i) {
			a[i] = mul(a[i], first);
			first = mul(first, step);
		}
	}
#if defined(__x86_64__) || defined(__i386__)
	__attribute__((target("avx2"))) static void add_scaled8(U *a, const U *b, int n, U r) {
		U q = quotient(r);
		I8 v = _mm256_set1_epi32(r), vq = _mm256_set1_epi32(q);
		int i = 0;
		for (; i + 8 <= n; i += 8) store8(a + i, add8(load8(a + i), fixed_mul8(load8(b + i), v, vq)));
		for (; i < n; ++i) a[i] = add(a[i], fixed_mul(b[i], r, q));
	}
#endif
	static void add_scaled(U *a, const U *b, int n, U r) {
#if defined(__x86_64__) || defined(__i386__)
		if (n >= 8 && vectorized()) {
			add_scaled8(a, b, n, r);
			return;
		}
#endif
		int i = 0;
		U q = quotient(r);
#if defined(__ARM_NEON)
		for (; i + 4 <= n; i += 4) vst1q_u32(a + i, add4(vld1q_u32(a + i), fixed_mul4(vld1q_u32(b + i), vdupq_n_u32(r), vdupq_n_u32(q))));
#endif
		for (; i < n; ++i) a[i] = add(a[i], fixed_mul(b[i], r, q));
	}
	void leaf(U *a, const U *b, int k) const {
		U w = roots[k / 2], q = root_quotients[k / 2];
		if (k & 1) w = mod - w, q = ~q;
		U window[16], coefficients[8];
		for (int i = 0; i < 8; ++i) {
			U x = a[i] < mod ? a[i] : a[i] - mod;
			U z = fixed_mul(x, w, q);
			window[i] = z < mod ? z : z - mod;
			window[8 + i] = x;
			coefficients[i] = b[i] < mod ? b[i] : b[i] - mod;
		}
#if defined(__ARM_NEON)
		for (int j = 0; j < 8; j += 4) {
			uint64x2_t lo = vdupq_n_u64(0), hi = lo;
			for (int i = 0; i < 8; ++i) {
				I z = vld1q_u32(window + 8 - i + j);
				lo = vmlal_n_u32(lo, vget_low_u32(z), coefficients[i]);
				hi = vmlal_n_u32(hi, vget_high_u32(z), coefficients[i]);
			}
			uint32x2_t x = vmul_n_u32(vmovn_u64(lo), 998244351u), y = vmul_n_u32(vmovn_u64(hi), 998244351u);
			I z = vcombine_u32(vshrn_n_u64(vmlal_n_u32(lo, x, mod), 32), vshrn_n_u64(vmlal_n_u32(hi, y, mod), 32));
			z = vminq_u32(z, vsubq_u32(z, vdupq_n_u32(twice)));
			vst1q_u32(a + j, vminq_u32(z, vsubq_u32(z, vdupq_n_u32(mod))));
		}
#else
		for (int j = 0; j < 8; ++j) {
			unsigned long long sum = 0;
			for (int i = 0; i < 8; ++i) sum += (unsigned long long)window[8 - i + j] * coefficients[i];
			U x = reduce(sum);
			if (x >= twice) x -= twice;
			a[j] = x < mod ? x : x - mod;
		}
#endif
	}
	template <int n>
	void convolution_tile(U *a, U *b, int k) const {
		if constexpr (n == 8) leaf(a, b, k);
		else if constexpr (n == 16) {
			stage<false>(a, n / 2, n, k);
			stage<false>(b, n / 2, n, k);
			convolution_tile<n / 2>(a, b, 2 * k);
			convolution_tile<n / 2>(a + n / 2, b + n / 2, 2 * k + 1);
			stage<true>(a, n / 2, n, k);
		} else {
			stage_pair<false>(a, n / 4, n, k);
			stage_pair<false>(b, n / 4, n, k);
			for (int i = 0; i < 4; ++i) convolution_tile<n / 4>(a + i * (n / 4), b + i * (n / 4), 4 * k + i);
			stage_pair<true>(a, n / 4, n, k);
		}
	}
	void convolution(U *a, U *b, int n, int k = 0) const {
#if defined(__x86_64__) || defined(__i386__)
		if (vectorized()) {
			convolution8(a, b, n, k);
			return;
		}
#endif
		if (n == 256) {
			convolution_tile<256>(a, b, k);
			return;
		}
		if (n == 128) {
			convolution_tile<128>(a, b, k);
			return;
		}
		if (n == 64) {
			convolution_tile<64>(a, b, k);
			return;
		}
		stage_pair<false>(a, n / 4, n, k);
		stage_pair<false>(b, n / 4, n, k);
		for (int i = 0; i < 4; ++i) convolution(a + i * (n / 4), b + i * (n / 4), n / 4, 4 * k + i);
		stage_pair<true>(a, n / 4, n, k);
	}
#if defined(__x86_64__) || defined(__i386__)
	__attribute__((target("avx2"))) static void accumulate8(U *dst, const U *a, const U *b, const U *c, int n, bool minus) {
		for (int i = 0; i < n; i += 8) {
			I8 x = load8(b + i), y = load8(c + i);
			x = minus ? _mm256_sub_epi32(_mm256_add_epi32(x, _mm256_set1_epi32(twice)), y)
					  : _mm256_add_epi32(x, y);
			store8(dst + i, add8(load8(dst + i), mul8(load8(a + i), x)));
		}
	}
#endif
#if defined(__x86_64__) || defined(__i386__)
	template <int count>
	__attribute__((target("avx2"))) static void accumulate_batch8(U *dst, const U *const *a, const U *const *b,
																  const U *const *c, int n, bool minus) {
		static_assert(count >= 1 && count <= 3);
		I8 p = _mm256_set1_epi32(mod), pp = _mm256_set1_epi32(twice), ni = _mm256_set1_epi32(998244351);
		int offset = minus ? n : 0;
		for (int i = 0; i < n; i += 8) {
			I8 even = _mm256_setzero_si256(), odd = even;
#pragma GCC unroll 3
			for (int k = 0; k < count; ++k) {
				I8 x = load8(b[k] + offset + i), y = load8(c[k] + offset + i);
				x = minus ? _mm256_sub_epi32(_mm256_add_epi32(x, pp), y) : _mm256_add_epi32(x, y);
				I8 z = load8(a[k] + offset + i);
				even = _mm256_add_epi64(even, _mm256_mul_epu32(z, x));
				odd = _mm256_add_epi64(odd, _mm256_mul_epu32(_mm256_srli_epi64(z, 32), _mm256_srli_epi64(x, 32)));
			}
			even = _mm256_add_epi64(even, _mm256_mul_epu32(_mm256_mul_epu32(even, ni), p));
			odd = _mm256_add_epi64(odd, _mm256_mul_epu32(_mm256_mul_epu32(odd, ni), p));
			I8 value = _mm256_blend_epi32(_mm256_srli_epi64(even, 32), odd, 0xaa);
			value = _mm256_min_epu32(value, _mm256_sub_epi32(value, pp));
			store8(dst + i, add8(load8(dst + i), value));
		}
	}
#endif
#if defined(__x86_64__) || defined(__i386__)
	template <int count, bool first, bool canonical>
	__attribute__((target("avx2"))) static void accumulate_wide8(U *dst, const U *const *a, const U *const *b, const U *const *c, int n, bool minus) {
		static_assert(count >= 1 && count <= 7);
		I8 p = _mm256_set1_epi32(mod), pp = _mm256_set1_epi32(twice), ni = _mm256_set1_epi32(998244351);
		int offset = minus ? n : 0;
		for (int i = 0; i < n; i += 8) {
			I8 even = _mm256_setzero_si256(), odd = even;
#pragma GCC unroll 7
			for (int k = 0; k < count; ++k) {
				I8 x = load8(b[k] + offset + i), y = load8(c[k] + offset + i);
				x = minus ? _mm256_sub_epi32(_mm256_add_epi32(x, canonical ? p : pp), y) : _mm256_add_epi32(x, y);
				if constexpr (!canonical) x = _mm256_min_epu32(x, _mm256_sub_epi32(x, pp));
				I8 z = load8(a[k] + offset + i);
				even = _mm256_add_epi64(even, _mm256_mul_epu32(z, x));
				odd = _mm256_add_epi64(odd, _mm256_mul_epu32(_mm256_srli_epi64(z, 32), _mm256_srli_epi64(x, 32)));
			}
			even = _mm256_add_epi64(even, _mm256_mul_epu32(_mm256_mul_epu32(even, ni), p));
			odd = _mm256_add_epi64(odd, _mm256_mul_epu32(_mm256_mul_epu32(odd, ni), p));
			I8 value = _mm256_blend_epi32(_mm256_srli_epi64(even, 32), odd, 0xaa);
			value = _mm256_min_epu32(value, _mm256_sub_epi32(value, pp));
			if constexpr (count == 7) value = _mm256_min_epu32(value, _mm256_sub_epi32(value, pp));
			if constexpr (!first) value = add8(load8(dst + i), value);
			store8(dst + i, value);
		}
	}
#endif
	template <int count, bool first, bool canonical>
	static void accumulate_wide(U *dst, const U *const *a, const U *const *b, const U *const *c, int n, bool minus) {
#if defined(__x86_64__) || defined(__i386__)
		if (n >= 8 && vectorized()) {
			accumulate_wide8<count, first, canonical>(dst, a, b, c, n, minus);
			return;
		}
#endif
		if constexpr (first) std::fill_n(dst, n, 0u);
		int offset = minus ? n : 0;
		for (int k = 0; k < count; ++k) accumulate(dst, a[k] + offset, b[k] + offset, c[k] + offset, n, minus);
	}
	template <int count>
	static void accumulate_batch(U *dst, const U *const *a, const U *const *b, const U *const *c, int n, bool minus) {
#if defined(__x86_64__) || defined(__i386__)
		if (n >= 8 && vectorized()) {
			accumulate_batch8<count>(dst, a, b, c, n, minus);
			return;
		}
#endif
		int offset = minus ? n : 0;
		for (int k = 0; k < count; ++k) accumulate(dst, a[k] + offset, b[k] + offset, c[k] + offset, n, minus);
	}
	static void accumulate(U *dst, const U *a, const U *b, const U *c, int n, bool minus) {
#if defined(__x86_64__) || defined(__i386__)
		if (n >= 8 && vectorized()) {
			accumulate8(dst, a, b, c, n, minus);
			return;
		}
#endif
		int i = 0;
#if defined(__ARM_NEON)
		for (; i + 4 <= n; i += 4) {
			I x = vld1q_u32(b + i), y = vld1q_u32(c + i);
			x = minus ? sub4(x, y) : add4(x, y);
			vst1q_u32(dst + i, add4(vld1q_u32(dst + i), mul4(vld1q_u32(a + i), x)));
		}
#endif
		for (; i < n; ++i) dst[i] = add(dst[i], mul(a[i], minus ? sub(b[i], c[i]) : add(b[i], c[i])));
	}
	static void scale_copy(U *dst, const U *src, int length, U r) {
#if defined(__x86_64__) || defined(__i386__)
		if (length >= 8 && vectorized()) {
			scale_copy8(dst, src, length, r);
			return;
		}
#endif
		U q = quotient(r);
		int i = 0;
#if defined(__ARM_NEON)
		for (; i + 4 <= length; i += 4) vst1q_u32(dst + i, fixed_mul4(vld1q_u32(src + i), vdupq_n_u32(r), vdupq_n_u32(q)));
#endif
		for (; i < length; ++i) dst[i] = fixed_mul(src[i], r, q);
	}
	void scale(U *a, int length, U r) const {
#if defined(__x86_64__) || defined(__i386__)
		if (length >= 8 && vectorized()) {
			scale8(a, length, r);
			return;
		}
#endif
		int i = 0;
		U q = quotient(r);
#if defined(__ARM_NEON)
		for (; i + 4 <= length; i += 4)
			vst1q_u32(a + i, fixed_mul4(vld1q_u32(a + i), vdupq_n_u32(r), vdupq_n_u32(q)));
#endif
		for (; i < length; ++i) a[i] = fixed_mul(a[i], r, q);
	}

	static void pair_product(U *a, const U *b, int length) {
#if defined(__x86_64__) || defined(__i386__)
		if (length >= 8 && vectorized()) {
			pair_product8(a, b, length);
			return;
		}
#endif
		int i = 0;
#if defined(__ARM_NEON)
		for (; i + 4 <= length; i += 4) {
			uint32x4x2_t v = vld2q_u32(b + 2 * i);
			vst1q_u32(a + i, mul4(v.val[0], v.val[1]));
		}
#endif
		for (; i < length; ++i) a[i] = mul(b[2 * i], b[2 * i + 1]);
	}
	static void pair_multiply(U *a, const U *b, int length) {
#if defined(__x86_64__) || defined(__i386__)
		if (length >= 8 && vectorized()) {
			pair_multiply8(a, b, length);
			return;
		}
#endif
		int i = 0;
#if defined(__ARM_NEON)
		for (; i + 4 <= length; i += 4) {
			uint32x4x2_t v = vld2q_u32(a + 2 * i);
			I x = mul4(v.val[1], vld1q_u32(b + i));
			v.val[1] = mul4(v.val[0], vld1q_u32(b + i));
			v.val[0] = x;
			vst2q_u32(a + 2 * i, v);
		}
#endif
		for (; i < length; ++i) {
			U x = a[2 * i];
			a[2 * i] = mul(b[i], a[2 * i + 1]);
			a[2 * i + 1] = mul(b[i], x);
		}
	}
	static void difference(U *a, const U *b, const U *c, int length) {
#if defined(__x86_64__) || defined(__i386__)
		if (length >= 8 && vectorized()) {
			difference8(a, b, c, length);
			return;
		}
#endif
		int i = 0;
#if defined(__ARM_NEON)
		for (; i + 4 <= length; i += 4)
			vst1q_u32(a + i, sub4(vld1q_u32(b + i), vld1q_u32(c + i)));
#endif
		for (; i < length; ++i) a[i] = sub(b[i], c[i]);
	}
	static void fold(U *a, int columns) {
#if defined(__x86_64__) || defined(__i386__)
		if (columns >= 8 && vectorized()) {
			fold8(a, columns);
			return;
		}
#endif
		int i = 0;
#if defined(__ARM_NEON)
		for (; i + 4 <= columns; i += 4) {
			I x = add4(vld1q_u32(a + i), vld1q_u32(a + columns + i));
			vst1q_u32(a + i, x);
			vst1q_u32(a + columns + i, x);
		}
#endif
		for (; i < columns; ++i) a[i] = a[columns + i] = add(a[i], a[columns + i]);
	}
	void project(U *a, int rows, int columns) const {
		if (columns == 1) {
#if defined(__x86_64__) || defined(__i386__)
			if (rows >= 8 && vectorized()) {
				project1_8(a, rows);
				return;
			}
#endif
			int j = 0;
#if defined(__ARM_NEON)
			for (; j + 4 <= rows; j += 4) {
				uint32x4x2_t v = vld2q_u32(a + 2 * j);
				v.val[0] = v.val[1] = add4(v.val[0], v.val[1]);
				vst2q_u32(a + 2 * j, v);
			}
#endif
			for (; j < rows; ++j) a[2 * j] = a[2 * j + 1] = add(a[2 * j], a[2 * j + 1]);
			return;
		}
		if (columns == 2) {
#if defined(__x86_64__) || defined(__i386__)
			if (rows >= 8 && vectorized()) {
				project2_8(a, rows);
				return;
			}
#endif
			int j = 0;
#if defined(__ARM_NEON)
			for (; j + 4 <= rows; j += 4) {
				uint32x4x4_t v = vld4q_u32(a + 4 * j);
				I x = add4(add4(v.val[0], v.val[1]), add4(v.val[2], v.val[3]));
				I y = add4(sub4(v.val[0], v.val[1]),
						   fixed_mul4(sub4(v.val[2], v.val[3]), vdupq_n_u32(inverse_roots[1]),
									  vdupq_n_u32(inverse_root_quotients[1])));
				v.val[0] = add4(x, y);
				v.val[1] = sub4(x, y);
				y = fixed_mul4(y, vdupq_n_u32(roots[1]), vdupq_n_u32(root_quotients[1]));
				v.val[2] = add4(x, y);
				v.val[3] = sub4(x, y);
				vst4q_u32(a + 4 * j, v);
			}
#endif
			for (; j < rows; ++j) {
				U *r = a + 4 * j;
				U x = add(add(r[0], r[1]), add(r[2], r[3]));
				U y = add(sub(r[0], r[1]), fixed_mul(sub(r[2], r[3]), inverse_roots[1],
													 inverse_root_quotients[1]));
				r[0] = add(x, y);
				r[1] = sub(x, y);
				y = fixed_mul(y, roots[1], root_quotients[1]);
				r[2] = add(x, y);
				r[3] = sub(x, y);
			}
			return;
		}
		if (columns == 4) {
#if defined(__x86_64__) || defined(__i386__)
			if (vectorized()) {
				project4_8(a, rows);
				return;
			}
#endif
#if defined(__ARM_NEON)
			for (int j = 0; j < rows; ++j) {
				I x = vld1q_u32(a + 8 * j), y = vld1q_u32(a + 8 * j + 4);
				small4<1, true>(x, y, 0);
				small4<2, true>(x, y, 0);
				x = y = add4(x, y);
				small4<2, false>(x, y, 0);
				small4<1, false>(x, y, 0);
				vst1q_u32(a + 8 * j, x);
				vst1q_u32(a + 8 * j + 4, y);
			}
			return;
#endif
		}
#if defined(__x86_64__) || defined(__i386__)
		if (columns == 8 && vectorized()) {
			project8_8(a, rows);
			return;
		}
#endif
#if defined(__x86_64__) || defined(__i386__)
		if (columns == 16 && vectorized()) {
			project16_8(a, rows);
			return;
		}
#endif
		for (int j = 0; j < rows; ++j) {
			U *row = a + 2 * columns * j;
			fft<true>(row, 2 * columns, columns);
			fold(row, columns);
			fft<false>(row, 2 * columns, columns);
		}
	}
	void composition_difference(U *a, const U *b, const U *c, int n) const {
#if defined(__x86_64__) || defined(__i386__)
		if (vectorized()) {
			composition_difference8(a, b, c, n);
			return;
		}
#endif
		difference(a, b, c, n);
	}
#if defined(__x86_64__) || defined(__i386__)
	__attribute__((target("avx2"))) static void composition_difference8(U *a, const U *b, const U *c, int n) {
		I8 p = _mm256_set1_epi32(twice);
		for (int i = 0; i < n; i += 8) {
			I8 x = load8(c + i);
			x = _mm256_min_epu32(x, _mm256_sub_epi32(x, p));
			store8(a + i, sub8(load8(b + i), x));
		}
	}
#endif
	void composition_finish(U *p, const U *v, int columns, U factor) const {
#if defined(__x86_64__) || defined(__i386__)
		if (columns >= 8 && vectorized()) {
			composition_finish8(p, v, columns, out(factor));
			return;
		}
#endif
		for (int i = 0; i < columns; ++i) p[i] = mul(sub(add(v[i], v[columns + i]), add(v[2 * columns + i], v[3 * columns + i])), factor);
	}
#if defined(__x86_64__) || defined(__i386__)
	__attribute__((target("avx2"))) static void composition_finish8(U *dst, const U *v, int n, U factor) {
		I8 r = _mm256_set1_epi32(factor), q = _mm256_set1_epi32(quotient(factor));
		for (int i = 0; i < n; i += 8) {
			I8 x = add8(load8(v + i), load8(v + n + i)), y = add8(load8(v + 2 * n + i), load8(v + 3 * n + i));
			store8(dst + i, fixed_mul8(sub8(x, y), r, q));
		}
	}
#endif
	struct CompositionDelete {
		void operator()(unsigned int *p) const { ::operator delete[](p, std::align_val_t(64)); }
	};
	using CompositionWords = std::unique_ptr<unsigned int[], CompositionDelete>;
	void composition_forward_y(U *a, int size, int columns) const {
#if defined(__x86_64__) || defined(__i386__)
		if (vectorized()) {
			int h = size;
			bool lazy = false;
			for (; h / 2 > columns && h / 2 >= 8; h >>= 2) {
				radix4_8<false>(a, h / 2, 2 * size);
				lazy = true;
			}
			if (lazy && h > columns) composition_shrink8(a, 2 * size);
			for (; h / 2 > columns; h >>= 2) stage_pair<false>(a, h / 2, 2 * size);
			if (h > columns) stage<false>(a, h, 2 * size);
			return;
		}
#endif
		int h = size;
		for (; h / 2 > columns; h >>= 2) stage_pair<false>(a, h / 2, 2 * size);
		if (h > columns) stage<false>(a, h, 2 * size);
	}
#if defined(__x86_64__) || defined(__i386__)
	__attribute__((target("avx2"))) static void composition_shrink8(U *a, int n) {
		I8 p = _mm256_set1_epi32(twice);
		for (int i = 0; i < n; i += 8) {
			I8 x = load8(a + i);
			store8(a + i, _mm256_min_epu32(x, _mm256_sub_epi32(x, p)));
		}
	}
#endif
	void double_y(U *a, int rows, int columns, bool inverse) const {
		const int size = rows * columns;
		int h = 2 * columns;
		for (; h <= size / 2; h <<= 2) {
#if defined(__x86_64__) || defined(__i386__)
			if (h >= 8 && vectorized()) radix4_8<true>(a, h, 2 * size);
			else
#endif
				stage_pair<true>(a, h, 2 * size);
		}
		if (h <= size) stage<true>(a, h, 2 * size);
#if defined(__x86_64__) || defined(__i386__)
		if (columns >= 2 && vectorized()) {
			scale_rows8(a, rows, columns, (inverse ? inverse_twists : twists).data() + rows);
			return;
		}
#endif
		for (int j = 0; j < rows; ++j)
			scale(a + 2 * columns * j, 2 * columns,
				  (inverse ? inverse_twists : twists)[rows + j]);
	}
#if defined(__x86_64__) || defined(__i386__)
	__attribute__((target("avx2"))) static void projection_dot8(U *p, const U *q, int n) {
		I8 ni = _mm256_set1_epi32(998244351), modulus = _mm256_set1_epi32(mod), pp = _mm256_set1_epi32(twice);
		I8 index = _mm256_setr_epi32(0, 2, 4, 6, 1, 3, 5, 7);
		for (int i = 0; i < n; i += 8) {
			I8 x = load8(p + 2 * i), y = load8(p + 2 * i + 8), a = load8(q + 2 * i), b = load8(q + 2 * i + 8);
			x = _mm256_min_epu32(x, _mm256_sub_epi32(x, pp));
			y = _mm256_min_epu32(y, _mm256_sub_epi32(y, pp));
			a = _mm256_min_epu32(a, _mm256_sub_epi32(a, pp));
			b = _mm256_min_epu32(b, _mm256_sub_epi32(b, pp));
			I8 lo = _mm256_add_epi64(_mm256_mul_epu32(x, _mm256_srli_epi64(a, 32)), _mm256_mul_epu32(_mm256_srli_epi64(x, 32), a));
			I8 hi = _mm256_add_epi64(_mm256_mul_epu32(y, _mm256_srli_epi64(b, 32)), _mm256_mul_epu32(_mm256_srli_epi64(y, 32), b));
			lo = _mm256_add_epi64(lo, _mm256_mul_epu32(_mm256_mul_epu32(lo, ni), modulus));
			hi = _mm256_add_epi64(hi, _mm256_mul_epu32(_mm256_mul_epu32(hi, ni), modulus));
			I8 r = _mm256_permutevar8x32_epi32(_mm256_blend_epi32(_mm256_srli_epi64(lo, 32), hi, 0xaa), index);
			store8(p + i, _mm256_min_epu32(r, _mm256_sub_epi32(r, pp)));
		}
	}
#endif

	static void projection_dot(U *p, const U *q, int n) {
#if defined(__x86_64__) || defined(__i386__)
		if (n >= 8 && vectorized()) {
			projection_dot8(p, q, n);
			return;
		}
#endif
		for (int i = 0; i < n; ++i) {
			U x = p[2 * i], y = p[2 * i + 1], a = q[2 * i], b = q[2 * i + 1];
			if (x >= twice) x -= twice;
			if (y >= twice) y -= twice;
			if (a >= twice) a -= twice;
			if (b >= twice) b -= twice;
			U r = reduce((unsigned long long)x * b + (unsigned long long)y * a);
			p[i] = r < twice ? r : r - twice;
		}
	}
	void transpose_roots() {
		roots.swap(inverse_roots);
		root_quotients.swap(inverse_root_quotients);
		cube_roots.swap(inverse_cube_roots);
		cube_quotients.swap(inverse_cube_quotients);
	}
	U reversion_projection(U *p, U *q, int size) {
		int levels = __builtin_ctz((unsigned int)size);
		V saved(levels, one);
		U q_one = one, p_one = one;
		for (int i = 1, columns = size / 2; columns > 1; ++i, columns >>= 1) {
			q_one = mul(mul(q_one, q_one), in(2 * columns));
			saved[i] = q_one;
		}
		for (int i = levels - 1, columns = 2; i >= 0; --i, columns <<= 1) {
			p_one = mul(p_one, saved[i]);
			if (i) p_one = mul(p_one, in(4 * columns));
		}
		U factor = T::word(mul(p_one, in(4 * size))).inv().v;
		fft<false>(q, 2 * size);
		for (int i = 0; i < 2 * size; ++i) {
			U x = q[i];
			q[i] = add(x, one);
			q[2 * size + i] = sub(x, one);
		}

		for (int i = 0; i < size; ++i) {
			U x = fixed_mul(sub(q[2 * i + 1], q[2 * i]), roots[i], root_quotients[i]);
			p[i] = i < size / 2 ? x : sub(0, x);
			p[size + i] = sub(0, p[i]);
		}
		q_one = one;
		for (int columns = size / 2, rows = 2; columns > 1; columns >>= 1, rows <<= 1) {
			pair_product(q, q, 2 * size);
			q_one = mul(q_one, q_one);
			project(q, rows, columns);
			std::copy_n(q, 2 * size, q + 2 * size);
			double_y(q + 2 * size, rows, columns, false);
			q_one = mul(q_one, in(2 * columns));
			U correction = add(q_one, q_one);
			for (int i = 0; i < 2 * columns; ++i) q[2 * size + i] = sub(q[2 * size + i], correction);
			composition_forward_y(q + 2 * size, size, columns);
			transpose_roots();
			project(p, rows, columns);
			for (int i = 0; i < 2 * size; ++i) p[2 * size + i] = sub(0, p[i]);
			double_y(p + 2 * size, rows, columns, true);
			composition_forward_y(p + 2 * size, size, columns);
			transpose_roots();
			projection_dot(p, q, 2 * size);
		}
		for (int i = 0; i < size; ++i) p[i] = add(p[2 * i], p[2 * i + 1]);
		transpose_roots();
		fft<true>(p, size);
		transpose_roots();
		return factor;
	}

	void compose(U *p, const U *q, int rows, int columns, U &p_one, U &q_one, U *workspace = nullptr, int level = 0) const {
		const int size = rows * columns;
		if (columns == 1) {
			fft<false>(p, rows);
			for (int i = rows - 1; i >= 0; --i) p[2 * i] = p[2 * i + 1] = p[i];
			return;
		}
		CompositionWords storage;
		if (!workspace) {
			int levels = __builtin_ctz((unsigned int)columns);
			storage.reset(static_cast<unsigned int *>(::operator new[]((4 * size + 16) * levels * sizeof(unsigned int), std::align_val_t(64))));
			workspace = storage.get();
		}
		U *v = workspace + (4 * size + 16) * level;
		if (rows == 1) {
			std::fill_n(v + columns, columns, U(0));
			for (int i = 0; i < columns; ++i) v[i] = sub(0, q[i]);
			fft<false>(v, 2 * columns);
			for (int i = 0; i < 2 * columns; ++i) {
				U x = v[i];
				v[i] = add(x, one);
				v[i + 2 * columns] = sub(x, one);
			}
		} else {
			pair_product(v, q, 2 * size);
			q_one = mul(q_one, q_one);
			project(v, rows, columns);
			std::copy_n(v, 2 * size, v + 2 * size);
			U *upper = v + 2 * size;
			double_y(upper, rows, columns, false);
			q_one = mul(q_one, in(2 * columns));
			U correction = add(q_one, q_one);
			for (int i = 0; i < 2 * columns; ++i) upper[i] = sub(upper[i], correction);
			composition_forward_y(upper, size, columns);
		}
		U saved = q_one;
		compose(p, v, 2 * rows, columns / 2, p_one, q_one, workspace, level + 1);
		pair_multiply(v, p, 2 * size);
		p_one = mul(p_one, saved);
		if (rows == 1) {
			fft<true>(v, 2 * columns, columns);
			fft<true>(v + 2 * columns, 2 * columns, columns);
			U factor = in(T(out(mul(p_one, in(4 * size)))).inv().val());
			composition_finish(p, v, columns, factor);
		} else {
			U *upper = v + 2 * size;
			double_y(upper, rows, columns, true);
			composition_forward_y(upper, size, columns);
			p_one = add(p_one, p_one);
			composition_difference(p, v, upper, 2 * size);
			project(p, rows, columns);
			p_one = mul(p_one, in(2 * columns));
		}
	}
};
inline NTT &ntt_plan(int size) {
	static NTT t(1, false);
	if (static_cast<int>(t.roots.size()) < size) t = NTT(size, false);
	return t;
}

}

inline void ntt(std::vector<fps_mint> &a, bool inverse = false) {
	assert(!a.empty() && !(a.size() & (a.size() - 1)) && a.size() <= (1U << 23));
	int n = static_cast<int>(a.size());
	auto &t = fps_detail::ntt_plan(std::max(1, n / 2));
	auto reorder = [&] {
		for (int i = 1, j = 0; i < n; ++i) {
			int bit = n >> 1;
			for (; j & bit; bit >>= 1) j ^= bit;
			j ^= bit;
			if (i < j) std::swap(a[i], a[j]);
		}
	};
	static_assert(sizeof(fps_mint) == sizeof(fps_detail::NTT::U) && std::is_trivially_copyable_v<fps_mint>);
	auto *words = reinterpret_cast<fps_detail::NTT::U *>(a.data());
	if (inverse) {
		reorder();
		t.fft<true>(words, n);
		t.scale(words, n, fps_mint(n).inv().val());
	} else {
		t.fft<false>(words, n);
		reorder();
	}
}

#if defined(__GNUC__) && !defined(__clang__) && (defined(__x86_64__) || defined(__i386__))
#pragma GCC pop_options
#endif


#line 1 "components/small_vector.hpp"



#line 7 "components/small_vector.hpp"
#include <cstring>
#line 9 "components/small_vector.hpp"
#include <iterator>
#line 14 "components/small_vector.hpp"

template <class T, std::size_t N = 2>
class SmallVector {
	static_assert(N && std::is_trivially_copyable_v<T> && std::is_trivially_destructible_v<T>);
	static constexpr std::size_t flag = std::size_t(1) << (std::numeric_limits<std::size_t>::digits - 1);
	struct Heap { T *data; std::size_t capacity; };
	union Storage {
		T local[N];
		Heap heap;
		Storage() : local{} {}
		~Storage() {}
	} storage;
	std::size_t length = 0;
	bool large() const { return length & flag; }
	void set_size(std::size_t n) { length = (length & flag) | n; }

public:
	SmallVector() = default;
	explicit SmallVector(std::size_t n) { resize(n); }
	SmallVector(std::initializer_list<T> values) { assign(values.begin(), values.end()); }
	template <std::input_iterator It>
	SmallVector(It first, It last) { assign(first, last); }
	SmallVector(const SmallVector &b) { assign(b.begin(), b.end()); }
	SmallVector(SmallVector &&b) noexcept : length(b.length) {
		if (large()) storage.heap = b.storage.heap;
		else std::memcpy(storage.local, b.storage.local, size() * sizeof(T));
		b.length = 0;
		b.storage.local[0] = T{};
	}
	SmallVector &operator=(const SmallVector &b) {
		if (this != &b) assign(b.begin(), b.end());
		return *this;
	}
	SmallVector &operator=(SmallVector &&b) noexcept {
		if (this == &b) return *this;
		if (large()) delete[] storage.heap.data;
		length = b.length;
		if (large()) storage.heap = b.storage.heap;
		else std::memcpy(storage.local, b.storage.local, size() * sizeof(T));
		b.length = 0;
		b.storage.local[0] = T{};
		return *this;
	}
	~SmallVector() { if (large()) delete[] storage.heap.data; }
	std::size_t size() const { return length & ~flag; }
	std::size_t capacity() const { return large() ? storage.heap.capacity : N; }
	bool empty() const { return size() == 0; }
	T *data() { return large() ? storage.heap.data : storage.local; }
	const T *data() const { return large() ? storage.heap.data : storage.local; }
	T *begin() { return data(); }
	const T *begin() const { return data(); }
	T *end() { return data() + size(); }
	const T *end() const { return data() + size(); }
	T &operator[](std::size_t i) { return data()[i]; }
	const T &operator[](std::size_t i) const { return data()[i]; }
	T &back() { return data()[size() - 1]; }
	const T &back() const { return data()[size() - 1]; }
	void reserve(std::size_t n) {
		if (n <= capacity()) return;
		assert(n < flag / sizeof(T));
		std::size_t capacity = std::max(n, this->capacity() * 2);
		auto memory = std::make_unique_for_overwrite<T[]>(capacity);
		std::memcpy(memory.get(), data(), size() * sizeof(T));
		if (large()) delete[] storage.heap.data;
		storage.heap = {memory.release(), capacity};
		length |= flag;
	}
	void resize(std::size_t n) {
		std::size_t old = size();
		reserve(n);
		if (n > old) std::fill(data() + old, data() + n, T{});
		set_size(n);
	}
	void resize_for_overwrite(std::size_t n) { reserve(n); set_size(n); }
	void clear() { set_size(0); }
	void push_back(T value) {
		std::size_t n = size();
		if (n == capacity()) reserve(n + 1);
		data()[n] = value;
		set_size(n + 1);
	}
	void pop_back() { assert(!empty()); set_size(size() - 1); }
	void assign(std::size_t n, T value) {
		reserve(n);
		std::fill_n(n > N ? storage.heap.data : data(), n, value);
		set_size(n);
	}
	template <std::input_iterator It>
	void assign(It first, It last) {
		if constexpr (std::forward_iterator<It>) {
			std::size_t n = std::distance(first, last);
			reserve(n);
			std::copy(first, last, data());
			set_size(n);
		} else {
			clear();
			for (; first != last; ++first) push_back(*first);
		}
	}
	T *insert(const T *position, std::size_t count, T value) {
		std::size_t offset = position - data(), old = size();
		reserve(old + count);
		T *p = data();
		std::memmove(p + offset + count, p + offset, (old - offset) * sizeof(T));
		std::fill(p + offset, p + offset + count, value);
		set_size(old + count);
		return p + offset;
	}
	template <std::forward_iterator It>
	T *insert(const T *position, It first, It last) {
		std::size_t offset = position - data(), old = size(), count = std::distance(first, last);
		reserve(old + count);
		T *p = data();
		std::memmove(p + offset + count, p + offset, (old - offset) * sizeof(T));
		std::copy(first, last, p + offset);
		set_size(old + count);
		return p + offset;
	}
	T *erase(const T *first, const T *last) {
		std::size_t offset = first - data(), count = last - first;
		std::memmove(data() + offset, last, (end() - last) * sizeof(T));
		set_size(size() - count);
		return data() + offset;
	}
	friend bool operator==(const SmallVector &a, const SmallVector &b) {
		return a.size() == b.size() && std::equal(a.begin(), a.end(), b.begin());
	}
};


#line 1 "components/limb_arithmetic.hpp"



#line 6 "components/limb_arithmetic.hpp"
#if defined(__x86_64__) || defined(__i386__)
#include <immintrin.h>
#endif

#if defined(__GNUC__) && !defined(__clang__)
#pragma GCC push_options
#pragma GCC optimize("O3")
#endif

namespace limb_arithmetic {
using U = std::uint64_t;

#if defined(__x86_64__) || defined(__i386__)
template <U base, bool subtract>
__attribute__((target("avx2"))) inline U blocks(U *a, const U *b, std::size_t &offset, std::size_t n, U carry) {
	const auto limit = _mm256_set1_epi64x(base - 1), modulus = _mm256_set1_epi64x(base);
	const auto one = _mm256_set1_epi64x(1), zero = _mm256_setzero_si256();
	const auto indices = _mm256_setr_epi64x(0, 1, 2, 3);
	for (; offset + 16 <= n; offset += 16) {
		__m256i values[4];
		unsigned generate = 0, propagate = 0;
		for (int j = 0; j < 4; ++j) {
			auto x = _mm256_loadu_si256(reinterpret_cast<const __m256i *>(a + offset + 4 * j));
			auto y = _mm256_loadu_si256(reinterpret_cast<const __m256i *>(b + offset + 4 * j));
			__m256i g, p;
			if constexpr (subtract) {
				values[j] = _mm256_sub_epi64(x, y);
				g = _mm256_cmpgt_epi64(y, x);
				p = _mm256_cmpeq_epi64(x, y);
			} else {
				values[j] = _mm256_add_epi64(x, y);
				g = _mm256_cmpgt_epi64(values[j], limit);
				p = _mm256_cmpeq_epi64(values[j], limit);
			}
			generate |= unsigned(_mm256_movemask_pd(_mm256_castsi256_pd(g))) << (4 * j);
			propagate |= unsigned(_mm256_movemask_pd(_mm256_castsi256_pd(p))) << (4 * j);
		}
		unsigned incoming = (propagate + (generate << 1) + unsigned(carry)) ^ propagate;
		for (int j = 0; j < 4; ++j) {
			auto in = _mm256_and_si256(_mm256_srlv_epi64(_mm256_set1_epi64x(incoming >> (4 * j)), indices), one);
			auto out = _mm256_and_si256(_mm256_srlv_epi64(_mm256_set1_epi64x(incoming >> (4 * j + 1)), indices), one);
			auto correction = _mm256_and_si256(_mm256_sub_epi64(zero, out), modulus);
			auto result = subtract ? _mm256_add_epi64(_mm256_sub_epi64(values[j], in), correction)
								   : _mm256_sub_epi64(_mm256_add_epi64(values[j], in), correction);
			_mm256_storeu_si256(reinterpret_cast<__m256i *>(a + offset + 4 * j), result);
		}
		carry = incoming >> 16;
	}
	return carry;
}
#endif

template <U base>
inline U add(U *a, const U *b, std::size_t n, U carry = 0) {
	static_assert(base > 1 && base <= (U(1) << 62));
	std::size_t i = 0;
#if defined(__x86_64__) || defined(__i386__)
	if (n >= 32 && __builtin_cpu_supports("avx2")) carry = blocks<base, false>(a, b, i, n, carry);
#endif
	for (; i < n; ++i) {
		U value = a[i] + b[i] + carry;
		carry = value >= base;
		a[i] = value - carry * base;
	}
	return carry;
}

template <U base>
inline U subtract(U *a, const U *b, std::size_t n, U borrow = 0) {
	static_assert(base > 1 && base <= (U(1) << 62));
	std::size_t i = 0;
#if defined(__x86_64__) || defined(__i386__)
	if (n >= 32 && __builtin_cpu_supports("avx2")) borrow = blocks<base, true>(a, b, i, n, borrow);
#endif
	for (; i < n; ++i) {
		U value = b[i] + borrow;
		borrow = a[i] < value;
		a[i] = a[i] + borrow * base - value;
	}
	return borrow;
}
}

#if defined(__GNUC__) && !defined(__clang__)
#pragma GCC pop_options
#endif


#line 8 "math/bigint.hpp"
#include <array>
#line 11 "math/bigint.hpp"
#include <charconv>
#include <compare>
#line 18 "math/bigint.hpp"
#include <streambuf>
#include <string>
#include <string_view>
#line 24 "math/bigint.hpp"

#if defined(__GNUC__) && !defined(__clang__)
#pragma GCC push_options
#pragma GCC optimize("O3")
#endif

class BigInt {
	using U = std::uint64_t;
	using W = __uint128_t;
	using D = SmallVector<U>;
	static constexpr U base = 1000000000000000000ULL;
	D a;
	bool negative = false;
	template <class I>
	static constexpr bool integral = std::is_integral_v<I> || std::is_same_v<I, __int128_t> || std::is_same_v<I, __uint128_t>;
	static unsigned read8(const char *p) {
		U x;
		std::memcpy(&x, p, 8);
		if constexpr (std::endian::native == std::endian::big) x = std::byteswap(x);
		assert(!(((x + 0x4646464646464646ULL) | (x - 0x3030303030303030ULL)) & 0x8080808080808080ULL));
		x &= 0x0f0f0f0f0f0f0f0fULL;
		x = (x * 10 + (x >> 8)) & 0x00ff00ff00ff00ffULL;
		x = (x * 100 + (x >> 16)) & 0x0000ffff0000ffffULL;
		return unsigned(x * 10000 + (x >> 32));
	}
	static constexpr auto pairs = [] {
		std::array<char, 200> result{};
		for (unsigned i = 0; i < 100; ++i) result[2 * i] = '0' + i / 10, result[2 * i + 1] = '0' + i % 10;
		return result;
	}();
	static void write4(char *p, unsigned value) {
		unsigned q = value / 100;
		std::memcpy(p, pairs.data() + 2 * q, 2);
		std::memcpy(p + 2, pairs.data() + 2 * (value - q * 100), 2);
	}
	static void write8(char *p, unsigned value) {
		unsigned q = value / 10000;
		write4(p, q);
		write4(p + 4, value - q * 10000);
	}
	template <bool subtract>
	bool combine_small(const BigInt &b) {
		if (a.size() > 1 || b.a.size() > 1) return false;
		std::int64_t x = a.empty() ? 0 : a[0], y = b.a.empty() ? 0 : b.a[0];
		if (negative) x = -x;
		if (b.negative != subtract) y = -y;
		std::int64_t sum = x + y;
		negative = sum < 0;
		U value = negative ? -sum : sum;
		a.resize_for_overwrite(value >= base ? 2 : value != 0);
		if (value >= base) a[0] = value - base, a[1] = 1;
		else if (value) a[0] = value;
		return true;
	}
	static void trim(D &v) {
		while (!v.empty() && !v.back()) v.pop_back();
	}
	void normalize() {
		trim(a);
		if (a.empty()) negative = false;
	}
	static int compare(const D &a, const D &b) {
		if (a.size() != b.size()) return a.size() < b.size() ? -1 : 1;
		for (std::size_t i = a.size(); i--;) if (a[i] != b[i]) return a[i] < b[i] ? -1 : 1;
		return 0;
	}
	static void add(D &a, const D &b) {
		std::size_t m = b.size();
		a.resize(std::max(a.size(), m));
		U carry = limb_arithmetic::add<base>(a.data(), b.data(), m);
		for (std::size_t i = m; carry; ++i) {
			if (i == a.size()) a.push_back(0);
			U value = a[i] + carry;
			carry = value >= base;
			a[i] = value - carry * base;
		}
	}
	static void subtract(D &a, const D &b) {
		U borrow = limb_arithmetic::subtract<base>(a.data(), b.data(), b.size());
		for (std::size_t i = b.size(); borrow; ++i) {
			borrow = !a[i];
			a[i] = a[i] + borrow * base - 1;
		}
		trim(a);
	}
	static void multiply_small(D &a, U b) {
		if (!b) return a.clear();
		W carry = 0;
		for (auto &x : a) {
			W value = W(x) * b + carry;
			carry = base_quotient(value);
			x = value - carry * base;
		}
		if (carry) a.push_back(carry);
	}
	static W high_product(W a, W b) {
		W low = W(U(a)) * U(b);
		W middle = W(U(a)) * U(b >> 64) + (low >> 64);
		W upper = W(U(a >> 64)) * U(b) + U(middle);
		return (a >> 64) * (b >> 64) + (middle >> 64) + (upper >> 64);
	}
	static W base_quotient(W value) {
		constexpr W inverse = std::numeric_limits<W>::max() / base;
		W q = high_product(value, inverse);
		return q + (value - q * base >= base);
	}
	static U divide_small(D &a, U b) {
		assert(b);
		if (b == 1) { trim(a); return 0; }
		W inverse = std::numeric_limits<W>::max() / b, remainder = 0;
		for (std::size_t i = a.size(); i--;) {
			W value = remainder * base + a[i];
			U q = high_product(value, inverse);
			remainder = value - W(q) * b;
			if (remainder >= b) ++q, remainder -= b;
			a[i] = q;
		}
		trim(a);
		return remainder;
	}
	static D schoolbook(const D &a, const D &b) {
		D c(a.size() + b.size());
		W carry = 0;
		bool square = &a == &b;
		for (std::size_t k = 0; k + 1 < c.size(); ++k) {
			std::size_t begin = k >= b.size() ? k + 1 - b.size() : 0;
			std::size_t end = std::min(k + 1, a.size());
			W value = carry;
			if (square) {
				for (std::size_t i = begin; i < std::min(end, (k + 1) / 2); ++i) value += W(a[i]) * a[k - i] * 2;
				if (!(k & 1) && k / 2 < a.size()) value += W(a[k / 2]) * a[k / 2];
			} else for (std::size_t i = begin; i < end; ++i) value += W(a[i]) * b[k - i];
			carry = base_quotient(value);
			c[k] = value - carry * base;
		}
		c.back() = carry;
		trim(c);
		return c;
	}
	static void karatsuba(const U *a, const U *b, W *c, int n, W *work, U *scratch, int na, int nb) {
		std::fill(c, c + 2 * n, W(0));
		if (!na || !nb) return;
		if (n <= 16) {
			if (a == b) {
				for (int i = 0; i < na; ++i) {
					c[2 * i] += W(a[i]) * a[i];
					for (int j = 0; j < i; ++j) c[i + j] += W(a[i]) * a[j] * 2;
				}
			} else for (int i = 0; i < na; ++i)
				for (int j = 0; j < nb; ++j) c[i + j] += W(a[i]) * b[j];
			return;
		}
		int m = n / 2;
		karatsuba(a, b, c, m, work, scratch, std::min(na, m), std::min(nb, m));
		karatsuba(a + m, b + m, c + n, m, work, scratch, std::max(0, na - m), std::max(0, nb - m));
		auto *x = scratch;
		auto *y = a == b ? x : x + m;
		for (int i = 0; i < m; ++i) x[i] = a[i] + a[i + m];
		if (x != y) for (int i = 0; i < m; ++i) y[i] = b[i] + b[i + m];
		karatsuba(x, y, work, m, work + n, scratch + n, std::min(na, m), std::min(nb, m));
		for (int i = 0; i < n; ++i) work[i] -= c[i] + c[i + n];
		for (int i = 0; i < n; ++i) c[i + m] += work[i];
	}
	static D karatsuba(const D &a, const D &b) {
		int n = std::bit_ceil(std::max(a.size(), b.size()));
		auto words = std::make_unique_for_overwrite<U[]>(4 * n);
		auto buffer = std::make_unique_for_overwrite<W[]>(4 * n);
		U *x = words.get(), *y = &a == &b ? x : x + n;
		std::copy(a.begin(), a.end(), x);
		std::fill(x + a.size(), x + n, U(0));
		if (x != y) {
			std::copy(b.begin(), b.end(), y);
			std::fill(y + b.size(), y + n, U(0));
		}
		W *result = buffer.get();
		karatsuba(x, y, result, n, result + 2 * n, x + 2 * n, a.size(), b.size());
		D c(a.size() + b.size());
		W carry = 0;
		for (std::size_t i = 0; i < c.size(); ++i) {
			carry += result[i];
			W q = base_quotient(carry);
			c[i] = carry - q * base;
			carry = q;
		}
		trim(c);
		return c;
	}
	static void add_shifted(D &a, const D &b, std::size_t shift) {
		a.resize(std::max(a.size(), b.size() + shift));
		U carry = 0;
		for (std::size_t i = 0; i < b.size() || carry; ++i) {
			if (i + shift == a.size()) a.push_back(0);
			U value = a[i + shift] + (i < b.size() ? b[i] : 0) + carry;
			carry = value >= base;
			a[i + shift] = value - carry * base;
		}
	}
	template <unsigned radix, int width>
	static bool centered_fits(const D &a, const D &b) {
		U sa = 0, sb = 0, ma = 0, mb = 0, ca = 0, cb = 0;
		auto scan = [](U value, U &carry, U &sum, U &maximum) {
			value += carry;
			for (int j = 0; j < width; ++j) {
				int digit = value % radix;
				value /= radix;
				int adjust = digit > int(radix / 2);
				digit -= adjust * int(radix);
				value += adjust;
				U magnitude = digit < 0 ? -digit : digit;
				sum += magnitude;
				maximum = std::max(maximum, magnitude);
			}
			carry = value;
		};
		for (std::size_t i = 0; i < std::max(a.size(), b.size()); ++i) {
			if (i < a.size()) scan(a[i], ca, sa, ma);
			if (i < b.size()) scan(b[i], cb, sb, mb);
			if (std::min(sa, sb) > (i + 1) * width * (radix / 16)) return false;
			if (std::min(W(sa) * mb, W(sb) * ma) * 2 >= fps_detail::NTT::mod) return false;
		}
		sa += ca; sb += cb;
		ma = std::max(ma, ca); mb = std::max(mb, cb);
		return std::min(W(sa) * mb, W(sb) * ma) * 2 < fps_detail::NTT::mod;
	}
	template <unsigned radix, int width, bool centered = false>
	static D transform_product(const D &a, const D &b) {
		using N = fps_detail::NTT;
		auto chunks = [](const D &v) {
			std::size_t n = (v.size() - 1) * width;
			for (U x = v.back(); x; x /= radix) ++n;
			return n;
		};
		std::size_t length = centered ? (a.size() + b.size()) * width + 1 : chunks(a) + chunks(b) - 1;
		if (length > (1U << 23)) {
			const D &u = a.size() >= b.size() ? a : b, &v = a.size() >= b.size() ? b : a;
			std::size_t half = u.size() / 2;
			D low(u.begin(), u.begin() + half), high(u.begin() + half, u.end());
			trim(low);
			D c = multiply(low, v), upper = multiply(high, v);
			add_shifted(c, upper, half);
			return c;
		}
		int n = std::bit_ceil(length);
		bool square = &a == &b;
		std::vector<std::uint32_t> p(n), q(square ? 0 : n);
		auto fill = [](auto &p, const D &v) {
			std::size_t offset = 0, length = 0;
			U carry = 0;
			for (U limb : v) {
				U value = limb + carry;
				for (int j = 0; j < width; ++j) {
					int digit = value % radix;
					value /= radix;
					if constexpr (centered) {
						int adjust = digit > int(radix / 2);
						digit -= adjust * int(radix);
						value += adjust;
					}
					p[offset++] = digit < 0 ? digit + int(N::mod) : digit;
					if (digit) length = offset;
				}
				carry = value;
			}
			if (carry) p[offset++] = carry, length = offset;
			return length;
		};
		std::size_t na = fill(p, a), nb = square ? na : fill(q, b);
		length = na + nb - 1;
		n = std::bit_ceil(length);
		p.resize(n);
		if (!square) q.resize(n);
		auto &t = fps_detail::ntt_plan(std::max(1, n / 2));
		t.fft<false>(p.data(), n);
		if (!square) t.fft<false>(q.data(), n);
		N::product(p.data(), p.data(), square ? p.data() : q.data(), n);
		t.fft<true>(p.data(), n);
		N::scale_normalize(p.data(), n, std::uint64_t(fps_mint(n).inv().val()) * N::one % N::mod);
		D c;
		c.reserve(a.size() + b.size());
		using C = std::conditional_t<centered, std::int64_t, U>;
		C carry = 0;
		U value = 0, scale = 1;
		for (std::size_t i = 0; i < length || carry; ++i) {
			C digit = i < length ? p[i] : 0;
			if constexpr (centered) digit -= C(digit > N::mod / 2) * N::mod;
			carry += digit;
			C remainder = carry % radix;
			carry /= radix;
			if constexpr (centered) {
				C borrow = remainder < 0;
				remainder += borrow * radix;
				carry -= borrow;
			}
			value += remainder * scale;
			scale *= radix;
			if (scale == base) c.push_back(value), value = 0, scale = 1;
		}
		if (value) c.push_back(value);
		trim(c);
		return c;
	}
	static D multiply(const D &a, const D &b) {
		if (a.empty() || b.empty()) return {};
		if (!a[0] || !b[0]) {
			std::size_t na = 0, nb = 0;
			while (!a[na]) ++na;
			while (!b[nb]) ++nb;
			D x(a.begin() + na, a.end()), y(b.begin() + nb, b.end());
			D c = multiply(x, y);
			c.insert(c.begin(), na + nb, 0);
			return c;
		}
		if (std::min(a.size(), b.size()) == 1) {
			D c = a.size() >= b.size() ? a : b;
			multiply_small(c, a.size() < b.size() ? a[0] : b[0]);
			return c;
		}
		if (std::min(a.size(), b.size()) <= 8 || std::max(a.size(), b.size()) <= 32) return schoolbook(a, b);
		if (std::max(a.size(), b.size()) <= 256) return karatsuba(a, b);
		std::size_t shorter = std::min(a.size(), b.size());
		if (centered_fits<1000, 6>(a, b)) return transform_product<1000, 6, true>(a, b);
		if (shorter * 9 <= 998244352 / 9801) return transform_product<100, 9>(a, b);
		if (centered_fits<100, 9>(a, b)) return transform_product<100, 9, true>(a, b);
		return transform_product<10, 18>(a, b);
	}
	static std::pair<D, D> divide_knuth(D a, D b) {
		std::size_t n = a.size(), m = b.size();
		U factor = base / (b.back() + 1);
		multiply_small(a, factor);
		multiply_small(b, factor);
		a.resize(n + 1);
		D q(n - m + 1);
		W inverse = std::numeric_limits<W>::max() / b.back();
		for (std::size_t j = q.size(); j--;) {
			W value = W(a[j + m]) * base + a[j + m - 1];
			U digit = high_product(value, inverse);
			if (value - W(digit) * b.back() >= b.back()) ++digit;
			digit = std::min(digit, base - 1);
			W remainder = value - W(digit) * b.back();
			while (remainder < base && W(digit) * b[m - 2] > remainder * base + a[j + m - 2]) --digit, remainder += b.back();
			W borrow = 0;
			for (std::size_t i = 0; i < m; ++i) {
				W product = W(digit) * b[i] + borrow;
				W next = base_quotient(product);
				U low = product - next * base;
				borrow = next + (a[j + i] < low);
				a[j + i] = a[j + i] + (a[j + i] < low ? base : 0) - low;
			}
			bool over = a[j + m] < borrow;
			a[j + m] = a[j + m] + (over ? base : 0) - borrow;
			if (over) {
				--digit;
				U carry = 0;
				for (std::size_t i = 0; i < m; ++i) {
					U sum = a[j + i] + b[i] + carry;
					carry = sum >= base;
					a[j + i] = sum - carry * base;
				}
				a[j + m] += carry;
			}
			q[j] = digit;
		}
		a.resize(m);
		divide_small(a, factor);
		trim(q);
		return {std::move(q), std::move(a)};
	}
	static D reciprocal(const D &a) {
		std::size_t n = a.size();
		if (n <= 8) {
			D power(2 * n + 1);
			power.back() = 1;
			return divide_knuth(std::move(power), a).first;
		}
		std::size_t half = (n + 1) / 2;
		D high(a.end() - half, a.end());
		D x = reciprocal(high);
		x.insert(x.begin(), n - half, 0);
		D product = multiply(a, x), power(2 * n + 1);
		power.back() = 1;
		int order = compare(power, product);
		if (order >= 0) subtract(power, product);
		else subtract(product, power), power = std::move(product);
		D correction = multiply(x, power);
		bool rounded = false;
		if (order < 0) {
			for (std::size_t i = 0; i < std::min(2 * n, correction.size()); ++i) rounded |= correction[i] != 0;
		}
		if (correction.size() <= 2 * n) correction.clear();
		else correction.erase(correction.begin(), correction.begin() + 2 * n);
		if (rounded) add(correction, D{1});
		if (order >= 0) add(x, correction);
		else subtract(x, correction);
		product = multiply(a, x);
		power.assign(2 * n + 1, 0);
		power.back() = 1;
		assert(compare(power, product) >= 0);
		subtract(power, product);
		while (compare(power, a) >= 0) subtract(power, a), add(x, D{1});
		return x;
	}
	static std::pair<D, D> divide(const D &a, const D &b) {
		assert(!b.empty());
		if (compare(a, b) < 0) return {{}, a};
		if (b.size() == 1) {
			D q = a;
			U remainder = divide_small(q, b[0]);
			return {std::move(q), remainder ? D{remainder} : D{}};
		}
		if (b.size() <= 16 || a.size() - b.size() <= 16) return divide_knuth(a, b);
		U factor = base / (b.back() + 1);
		D u = a, v = b;
		multiply_small(u, factor);
		multiply_small(v, factor);
		D inverse = reciprocal(v), quotient(u.size()), remainder;
		std::size_t n = v.size(), blocks = (u.size() + n - 1) / n;
		for (std::size_t block = blocks; block--;) {
			std::size_t begin = block * n, end = std::min(begin + n, u.size());
			remainder.insert(remainder.begin(), u.begin() + begin, u.begin() + end);
			trim(remainder);
			if (compare(remainder, v) < 0) continue;
			D q = multiply(remainder, inverse);
			if (q.size() <= 2 * n) q.clear();
			else q.erase(q.begin(), q.begin() + 2 * n);
			D product = multiply(q, v);
			subtract(remainder, product);
			while (compare(remainder, v) >= 0) subtract(remainder, v), add(q, D{1});
			std::copy(q.begin(), q.end(), quotient.begin() + begin);
		}
		trim(quotient);
		divide_small(remainder, factor);
		return {std::move(quotient), std::move(remainder)};
	}

public:
	BigInt() = default;
	BigInt(const BigInt &) = default;
	BigInt &operator=(const BigInt &) = default;
	BigInt(BigInt &&b) noexcept : a(std::move(b.a)), negative(std::exchange(b.negative, false)) {}
	BigInt &operator=(BigInt &&b) noexcept {
		if (this != &b) a = std::move(b.a), negative = std::exchange(b.negative, false);
		return *this;
	}
	template <class I> requires integral<I>
	BigInt(I value) {
		if constexpr (std::numeric_limits<I>::is_signed) negative = value < 0;
		__uint128_t magnitude = value;
		if (negative) magnitude = -magnitude;
		while (magnitude) a.push_back(magnitude % base), magnitude /= base;
	}
	explicit BigInt(std::string_view value) {
		assert(!value.empty());
		std::size_t begin = value[0] == '-' || value[0] == '+';
		negative = value[0] == '-';
		assert(begin < value.size());
		a.resize_for_overwrite((value.size() - begin + 17) / 18);
		std::size_t word = 0, end = value.size();
		for (; end - begin >= 18; end -= 18) {
			const char *p = value.data() + end - 18;
			assert(p[16] >= '0' && p[16] <= '9' && p[17] >= '0' && p[17] <= '9');
			a[word++] = (U(read8(p)) * 100000000 + read8(p + 8)) * 100 + (p[16] - '0') * 10 + p[17] - '0';
		}
		if (end != begin) {
			U digit = 0;
			std::size_t i = begin;
			for (; i + 8 <= end; i += 8) digit = digit * 100000000 + read8(value.data() + i);
			for (; i < end; ++i) {
				assert(value[i] >= '0' && value[i] <= '9');
				digit = digit * 10 + value[i] - '0';
			}
			a[word] = digit;
		}
		normalize();
	}
	std::size_t decimal_capacity() const { return a.empty() ? 1 : a.size() * 18 + negative; }
	char *to_chars(char *out) const {
		if (a.empty()) { *out++ = '0'; return out; }
		if (negative) *out++ = '-';
		out = std::to_chars(out, out + 18, a.back()).ptr;
		for (std::size_t i = a.size() - 1; i--;) {
			U value = a[i], q = value / 100000000;
			write8(out + 10, unsigned(value - q * 100000000));
			unsigned high = unsigned(q / 100000000);
			write8(out + 2, unsigned(q - U(high) * 100000000));
			std::memcpy(out, pairs.data() + 2 * high, 2);
			out += 18;
		}
		return out;
	}
	std::string str() const {
		std::string result;
		result.resize_and_overwrite(decimal_capacity(), [this](char *out, std::size_t) {
			return std::size_t(to_chars(out) - out);
		});
		return result;
	}
	int sign() const { return a.empty() ? 0 : negative ? -1 : 1; }
	explicit operator bool() const { return !a.empty(); }
	template <class I> requires integral<I>
	explicit operator I() const {
		__uint128_t value = 0;
		for (std::size_t i = a.size(); i--;) value = value * base + a[i];
		return I(negative ? -value : value);
	}
	friend bool operator==(const BigInt &x, const BigInt &y) { return x.negative == y.negative && x.a == y.a; }
	friend std::strong_ordering operator<=>(const BigInt &x, const BigInt &y) {
		int order = x.negative != y.negative ? (x.negative ? -1 : 1) : compare(x.a, y.a) * (x.negative ? -1 : 1);
		return order <=> 0;
	}
	BigInt operator-() const {
		BigInt result = *this;
		if (!a.empty()) result.negative = !negative;
		return result;
	}
	BigInt operator+() const { return *this; }
	BigInt &operator+=(const BigInt &b) {
		if (combine_small<false>(b)) return *this;
		if (negative == b.negative) add(a, b.a);
		else if (compare(a, b.a) >= 0) subtract(a, b.a);
		else {
			D result = b.a;
			subtract(result, a);
			a = std::move(result);
			negative = b.negative;
		}
		if (a.empty()) negative = false;
		return *this;
	}
	BigInt &operator-=(const BigInt &b) {
		if (combine_small<true>(b)) return *this;
		if (negative != b.negative) add(a, b.a);
		else if (compare(a, b.a) >= 0) subtract(a, b.a);
		else {
			D result = b.a;
			subtract(result, a);
			a = std::move(result);
			negative = !negative;
		}
		if (a.empty()) negative = false;
		return *this;
	}
	BigInt &operator*=(const BigInt &b) {
		bool sign = negative != b.negative;
		a = multiply(a, b.a);
		negative = !a.empty() && sign;
		return *this;
	}
	template <class I> requires integral<I>
	BigInt &operator*=(I value) {
		bool sign = false;
		if constexpr (std::numeric_limits<I>::is_signed) sign = value < 0;
		__uint128_t magnitude = value;
		if (sign) magnitude = -magnitude;
		if (magnitude >= base) return *this *= BigInt(value);
		multiply_small(a, U(magnitude));
		negative = !a.empty() && (negative != sign);
		return *this;
	}
	friend std::pair<BigInt, BigInt> divmod(const BigInt &x, const BigInt &y) {
		auto [q, r] = divide(x.a, y.a);
		BigInt quotient, remainder;
		quotient.a = std::move(q);
		remainder.a = std::move(r);
		quotient.negative = !quotient.a.empty() && x.negative != y.negative;
		remainder.negative = !remainder.a.empty() && x.negative;
		return {std::move(quotient), std::move(remainder)};
	}
	BigInt &operator/=(const BigInt &b) {
		*this = divmod(*this, b).first;
		return *this;
	}
	BigInt &operator%=(const BigInt &b) {
		*this = divmod(*this, b).second;
		return *this;
	}
	friend BigInt operator+(BigInt a, const BigInt &b) { a += b; return a; }
	friend BigInt operator-(BigInt a, const BigInt &b) { a -= b; return a; }
	friend BigInt operator*(const BigInt &a, const BigInt &b) {
		BigInt result;
		result.a = multiply(a.a, b.a);
		result.negative = !result.a.empty() && a.negative != b.negative;
		return result;
	}
	template <class I> requires integral<I>
	friend BigInt operator*(BigInt a, I b) { a *= b; return a; }
	template <class I> requires integral<I>
	friend BigInt operator*(I a, BigInt b) { b *= a; return b; }
	friend BigInt operator/(const BigInt &a, const BigInt &b) { return divmod(a, b).first; }
	friend BigInt operator%(const BigInt &a, const BigInt &b) { return divmod(a, b).second; }
	BigInt &operator++() { return *this += 1; }
	BigInt operator++(int) { BigInt old = *this; ++*this; return old; }
	BigInt &operator--() { return *this -= 1; }
	BigInt operator--(int) { BigInt old = *this; --*this; return old; }
	friend BigInt abs(BigInt value) { value.negative = false; return value; }
	friend BigInt pow(BigInt value, std::uint64_t exponent) {
		BigInt result = 1;
		for (; exponent; exponent >>= 1) {
			if (exponent & 1) result *= value;
			if (exponent > 1) value *= value;
		}
		return result;
	}
	friend std::istream &operator>>(std::istream &in, BigInt &value) {
		std::istream::sentry sentry(in);
		if (!sentry) return in;
		auto state = std::ios_base::goodbit;
		try {
			std::string text;
			auto *buffer = in.rdbuf();
			auto c = buffer->sgetc();
			if (c == '+' || c == '-') {
				text.push_back(static_cast<char>(c));
				c = buffer->snextc();
			}
			auto start = text.size();
			while (c >= '0' && c <= '9') {
				text.push_back(static_cast<char>(c));
				c = buffer->snextc();
			}
			if (std::char_traits<char>::eq_int_type(c, std::char_traits<char>::eof())) state |= std::ios_base::eofbit;
			if (text.size() == start) {
				value = 0;
				state |= std::ios_base::failbit;
			} else value = BigInt(text);
		} catch (...) {
			bool rethrow = (in.exceptions() & std::ios_base::badbit) != 0;
			try { in.setstate(std::ios_base::badbit); } catch (...) {}
			if (rethrow) throw;
			return in;
		}
		if (state != std::ios_base::goodbit) in.setstate(state);
		return in;
	}
	friend std::ostream &operator<<(std::ostream &out, const BigInt &value) { return out << value.str(); }
};

#if defined(__GNUC__) && !defined(__clang__)
#pragma GCC pop_options
#endif


#line 1 "components/fast_io.hpp"



#line 9 "components/fast_io.hpp"
#include <cstdio>
#line 12 "components/fast_io.hpp"
#include <stdexcept>
#line 16 "components/fast_io.hpp"
#if defined(__unix__) || defined(__APPLE__)
#include <sys/mman.h>
#include <sys/stat.h>
#include <unistd.h>
#endif

class FastInput {
	static constexpr std::size_t capacity = 1 << 16;
	FILE *file = nullptr;
	std::unique_ptr<char[]> buffer;
	const char *current = "", *end = current;
	std::string scratch;
	void *mapping = nullptr;
	std::size_t mapped_size = 0;
	bool refill() {
		if (!file) return false;
		std::size_t n = std::fread(buffer.get(), 1, capacity, file);
		current = buffer.get();
		end = current + n;
		return n;
	}
	static const char *scan(const char *p, const char *end) {
		while (end - p >= 8) {
			std::uint64_t x;
			std::memcpy(&x, p, 8);
			if ((x - 0x2121212121212121ULL) & ~x & 0x8080808080808080ULL) break;
			p += 8;
		}
		while (p != end && static_cast<unsigned char>(*p) > ' ') ++p;
		return p;
	}

public:
	explicit FastInput(FILE *input = stdin) {
#if defined(__unix__) || defined(__APPLE__)
		struct stat status;
		int fd = fileno(input);
		auto offset = ftello(input);
		if (fd >= 0 && offset >= 0 && !fstat(fd, &status) && S_ISREG(status.st_mode) && status.st_size > offset) {
			mapped_size = status.st_size;
			mapping = mmap(nullptr, mapped_size, PROT_READ, MAP_PRIVATE, fd, 0);
			if (mapping != MAP_FAILED) {
				current = static_cast<const char *>(mapping) + offset;
				end = static_cast<const char *>(mapping) + mapped_size;
				return;
			}
			mapping = nullptr;
		}
#endif
		file = input;
		buffer = std::make_unique_for_overwrite<char[]>(capacity);
	}
	explicit FastInput(std::string_view input) : current(input.empty() ? "" : input.data()), end(current + input.size()) {}
	FastInput(const FastInput &) = delete;
	FastInput &operator=(const FastInput &) = delete;
	~FastInput() {
#if defined(__unix__) || defined(__APPLE__)
		if (mapping) munmap(mapping, mapped_size);
#endif
	}
	std::string_view token() {
		while (true) {
			while (current != end && static_cast<unsigned char>(*current) <= ' ') ++current;
			if (current != end) break;
			if (!refill()) return {};
		}
		const char *begin = current;
		current = scan(current, end);
		if (current != end || !file) return {begin, std::size_t(current - begin)};
		scratch.assign(begin, end);
		while (refill()) {
			current = scan(current, end);
			scratch.append(buffer.get(), std::size_t(current - buffer.get()));
			if (current != end) break;
		}
		return scratch;
	}
	template <class T>
	bool read(T &value) {
		if constexpr (std::is_integral_v<T> && !std::is_same_v<T, bool>) {
			while (true) {
				while (current != end && static_cast<unsigned char>(*current) <= ' ') ++current;
				if (current != end) break;
				if (!refill()) return false;
			}
			T parsed{};
			auto result = std::from_chars(current, end, parsed);
			if (result.ec == std::errc{} && (result.ptr != end ? static_cast<unsigned char>(*result.ptr) <= ' ' : !file)) {
				current = result.ptr;
				value = parsed;
				return true;
			}
		}
		auto s = token();
		if (s.empty()) return false;
		if constexpr (std::is_same_v<T, bool>) {
			assert(s == "0" || s == "1");
			value = s[0] == '1';
		} else if constexpr (std::is_integral_v<T>) {
			auto result = std::from_chars(s.data(), s.data() + s.size(), value);
			assert(result.ec == std::errc{} && result.ptr == s.data() + s.size());
			(void)result;
		} else value = T(s);
		return true;
	}
};

class FastOutput {
	FILE *file;
	std::unique_ptr<char[]> buffer = std::make_unique_for_overwrite<char[]>(1 << 16);
	std::size_t capacity = 1 << 16, length = 0;
	char *reserve(std::size_t n) {
		if (n > capacity - length) {
			if (!flush()) throw std::runtime_error("output failed");
			if (n > capacity) {
				std::size_t next = std::max(n, 2 * capacity);
				buffer = std::make_unique_for_overwrite<char[]>(next);
				capacity = next;
			}
		}
		return buffer.get() + length;
	}

public:
	explicit FastOutput(FILE *output = stdout) : file(output) {}
	FastOutput(const FastOutput &) = delete;
	FastOutput &operator=(const FastOutput &) = delete;
	~FastOutput() { flush(); }
	bool flush() {
		std::size_t offset = 0;
		while (offset < length) {
			auto n = std::fwrite(buffer.get() + offset, 1, length - offset, file);
			if (!n) {
				std::memmove(buffer.get(), buffer.get() + offset, length - offset);
				length -= offset;
				return false;
			}
			offset += n;
		}
		length = 0;
		return std::fflush(file) == 0;
	}
	void write(char c) { *reserve(1) = c; ++length; }
	void write(__uint128_t value) {
		char digits[40];
		char *end = digits + sizeof(digits), *begin = end;
		do { *--begin = char('0' + value % 10); value /= 10; } while (value);
		write(std::string_view(begin, end - begin));
	}
	void write(__int128_t value) {
		__uint128_t magnitude = __uint128_t(value);
		if (value < 0) write('-'), magnitude = 0 - magnitude;
		write(magnitude);
	}
	void write(std::string_view s) {
		char *out = reserve(s.size());
		if (!s.empty()) std::memcpy(out, s.data(), s.size());
		length += s.size();
	}
	template <class T> requires std::is_integral_v<T>
	void write(T value) {
		char *out = reserve(32);
		if constexpr (std::is_same_v<T, bool>) *out = '0' + value, ++length;
		else length += std::to_chars(out, out + 32, value).ptr - out;
	}
	template <class T> requires requires(const T &value, char *out) { value.decimal_capacity(); value.to_chars(out); }
	void write(const T &value) {
		char *out = reserve(value.decimal_capacity());
		length += value.to_chars(out) - out;
	}
	template <class T>
	void println(const T &value) { write(value); write('\n'); }
};


#line 5 "test/dijkstra_bigint_yuki.test.cpp"
#include <numeric>

int main() {
	FastInput input;
	FastOutput output;
	int n = 0, m = 0;
	input.read(n), input.read(m);
	struct Edge { int u, v, a, b; };
	std::vector<Edge> edges(m);
	int max_b = 1;
	for (auto& e : edges) input.read(e.u), input.read(e.v), input.read(e.a), input.read(e.b), max_b = std::max(max_b, e.b);
	BigInt denominator = 1;
	for (int b = 2; b <= max_b; ++b) denominator *= b / std::gcd(b, int(denominator % BigInt(b)));
	std::vector<std::vector<std::pair<int, BigInt>>> graph(n);
	std::vector<BigInt> unit(max_b + 1);
	for (int b = 1; b <= max_b; ++b) unit[b] = denominator / BigInt(b);
	for (auto e : edges) {
		auto weight = unit[e.b] * BigInt(e.a);
		graph[e.u - 1].emplace_back(e.v - 1, weight);
		graph[e.v - 1].emplace_back(e.u - 1, std::move(weight));
	}
	auto distance = dijkstra(graph, 0, denominator * BigInt(300LL * n + 1));
	for (int v = 1; v < n; ++v) {
		BigInt a = distance[v], b = denominator;
		while (b) { auto r = a % b; a = std::move(b), b = std::move(r); }
		output.write(distance[v] / a), output.write(' '), output.println(denominator / a);
	}
}