代码

#if IN_LOCAL
import local;
#endif // IN_LOCAL
#include<vector>
#include<cstring>
#include<iostream>
#include<algorithm>
#include<functional>
using namespace std;
namespace ModularPolynomial
{
#define CE
#define SC(a,b) static_cast<a>(b)
	using u = unsigned; using ull = unsigned long long;
	CE const u mod = 998244353, G = 3189, modM1 = mod - 1, modM2 = mod - 2, modM1D2 = modM1 / 2, mod2M2 = modM1 * 2; CE const ull modmod2 = SC(ull, mod) * mod * 2;
	inline CE u Power(u a, u b) { u ans = 1; for (; b; a = SC(ull, a) * a % mod, b /= 2) if (b % 2) ans = ans * SC(ull, a) % mod; return ans; }
	inline CE u Inv(const u val) { return Power(val, modM2); } inline CE u Div2(const u val) { return (val % 2 ? val + mod : val) / 2; }
	inline CE void UpdateInc(u& pos, const u val) { if ((pos += val) >= mod) pos -= mod; return; } inline CE u UnderflowReduce(const u val) { return SC(int, val) < 0 ? (val + mod) : val; }
	inline CE void UpdateDec(u& pos, const u val) { if (SC(int, pos -= val) < 0) pos += mod; return; } inline CE u CeilLog2(const u val) { return val > 1 ? 32 - __builtin_clz(val - 1) : 1; }
	struct ArrayPointer { void* p; template<class T1, class T2> ArrayPointer(std::vector<T1, T2>& v) :p(v.data()) {} template<class T> ArrayPointer(T* const _) : p(_) {} };
	struct ConstArrayPointer
	{
		const void* p; template<class T1, class T2> inline ConstArrayPointer(const std::vector<T1, T2>& v) :p(v.data()) {} template<class T> inline ConstArrayPointer(const T* const _) : p(_) {}
	};
	inline void MemoryCopy(const ArrayPointer dest, const ConstArrayPointer source, const size_t size) { memcpy(dest.p, source.p, size); return; }
	inline void MemoryClear(ArrayPointer dest, const size_t size) { memset(dest.p, 0, size); return; }
	CE const u inv_G = Inv(G), sqrt_M1 = Power(G, modM1D2 / 2), inv_sqrtM1 = Inv(sqrt_M1); CE const size_t sizeof_u = sizeof(u); std::pair<u, u> gn[33]; std::vector<u> ur;
	struct Poly :public std::vector<u>
	{
		using base = std::vector<u>; using base::base; inline explicit Poly() {} inline explicit Poly(const u size) :base(size) {}
		inline void resize(const u new_size) { base::resize(new_size); return; } inline void resize(const u new_size, const u init) { base::resize(new_size, init); return; }
		inline void reserve(const u new_size) { base::reserve(new_size); return; } inline u size()const { return SC(u, base::size()); }
		inline void expand(const u new_size) { if (new_size > size()) resize(new_size); return; } inline u operator[](const u pos)const { return base::operator[](pos); }
		inline u& operator[](const u pos) { return base::operator[](pos); } inline u* operator+(const u d) { return data() + d; } inline const u* operator+(const u d)const { return data() + d; }
	};
	inline void Prepare(const u maxi_length)
	{
		if (maxi_length == 1) { ur.resize(2, 1); return; }
		const u siz_gn = CeilLog2(maxi_length), siz_ur = 1u << siz_gn; gn[siz_gn].first = Power(G, modM1D2 >> siz_gn); gn[siz_gn].second = Power(inv_G, modM1D2 >> siz_gn);
		for (u i = siz_gn; i > 1; gn[i - 1].first = SC(ull, gn[i].first) * gn[i].first % mod, gn[i - 1].second = SC(ull, gn[i].second) * gn[i].second % mod, --i);
		ur.resize(siz_ur); ur[siz_ur / 2] = 1; const u w = gn[siz_gn - 1].first; for (u i = siz_ur / 2 + 1; i != siz_ur; ur[i] = ur[i - 1] * SC(ull, w) % mod, ++i);
		for (u i = siz_ur / 2 - 1; i; ur[i] = ur[SC(size_t, i * 2)], --i); return;
	}
	inline void DIFNTT(const u length, u* const a)
	{
		for (u i = length, * const end_p = a + length; i != 1;) for (u i2 = i, *p = (i /= 2, a); p != end_p; p += i)
			for (u j = i, t; j != i2; t = *p, UpdateInc(*p, *(p + i)), *(p + i) = SC(ull, t + mod - *(p + i)) * ur[j++] % mod, ++p);
		return;
	}
	inline void DITNTT(const u length, u* const a)
	{
		for (u i = 1, i2 = 2, * const end_p = a + length; i != length; i = i2, i2 *= 2) for (u* p = a; p != end_p; p += i)
			for (u j = i, t; j != i2; *(p + i) = UnderflowReduce(*p - (t = SC(ull, *(p + i)) * ur[j++] % mod)), UpdateInc(*(p++), t)); return;
	}
	inline void DIFNTT(const u length, Poly& a) { DIFNTT(length, a.data()); return; } inline void DITNTT(const u length, Poly& a) { DITNTT(length, a.data()); return; }
	inline void Mul(const u ans_length, Poly& a, Poly& b)
	{
		const u log_l = CeilLog2(a.size() + b.size() - 1), l = 1u << log_l, inv_l = mod - (modM1 >> log_l); a.resize(l); b.resize(l); DIFNTT(l, a); DIFNTT(l, b);
		for (u i = 0; i < l; b[i] = SC(ull, a[i]) * b[i] % mod, ++i); DITNTT(l, b); a.front() = b.front() * SC(ull, inv_l) % mod;
		for (u i = 1, end_i = std::min(ans_length, l); i < end_i; a[i] = b[l - i] * SC(ull, inv_l) % mod, ++i); a.resize(ans_length); return;
	}
	inline void Inv(const u ans_length, Poly& a, Poly& ans)
	{
		static Poly aux1, aux2; const u l = 1u << CeilLog2(ans_length), limit = a.size(); a.expand(l); ans.resize(l); aux1.resize(l); aux2.resize(l); ans.front() = Inv(a.front());
		for (u i = 1, p_gn = 1; i < ans_length; ++p_gn)
		{
			const u i2 = i * 2, inv_i2 = mod - (modM1 >> p_gn), inv_i = mod - (mod2M2 >> p_gn), cur_gn = gn[p_gn].first, inv_cur_gn = gn[p_gn].second;
			if (i2 > limit) { MemoryCopy(aux1, a, limit * sizeof_u); MemoryClear(aux1 + limit, SC(size_t, i2 - limit) * sizeof_u); }
			else MemoryCopy(aux1, a, i2 * sizeof_u); DIFNTT(i2, aux1); MemoryCopy(aux2, ans, i * sizeof_u); MemoryClear(aux2 + i, i * sizeof_u);
			DIFNTT(i2, aux2); for (u j = 0; j < i2; aux2[j] = (aux1[j] * SC(ull, aux2[j]) + modM1) % mod * aux2[j] % mod, ++j); DITNTT(i2, aux2);
			aux1[0] = aux2[0] * SC(ull, inv_i2) % mod; for (u j = 1; j < i2; aux1[j] = aux2[i2 - j] * SC(ull, inv_i2) % mod, ++j);
			for (u j = 0, val = 1; j < i; aux2[j] = (a[j] + a[i | j] * SC(ull, sqrt_M1)) % mod * val % mod, ans[i | j] = ans[j] * SC(ull, val) % mod, ++j, val = val * SC(ull, cur_gn) % mod);
			DIFNTT(i, aux2); DIFNTT(i, ans + i); for (u j = 0; j < i; aux2[j] = (aux2[j] * SC(ull, ans[i | j]) + modM1) % mod * ans[i | j] % mod, ++j);
			DITNTT(i, aux2); ans[i] = Div2((modmod2 - (SC(ull, aux2.front()) * inv_i + aux1.front()) % mod * inv_sqrtM1 - aux1[i]) % mod); const u end_j = std::min(i, ans_length - i);
			for (u j = 1, val = cur_gn, inv_val = SC(ull, inv_sqrtM1) * inv_cur_gn % mod; j < end_j; ++j)
				ans[i | j] = Div2((modmod2 - (aux2[i - j] * SC(ull, inv_i) + SC(ull, aux1[j]) * val) % mod * inv_val - aux1[i | j]) % mod),
				val = val * SC(ull, cur_gn) % mod, inv_val = inv_val * SC(ull, inv_cur_gn) % mod; i = i2;
		}
		ans.resize(ans_length); return;
	}
	inline void Div(const u ans_length, Poly& a, Poly& b)
	{
		if (ans_length == 1) { a.resize(1); a.front() = a.front() * SC(ull, Inv(b.front())) % mod; return; }
		static Poly aux1, aux2; const u log_l = CeilLog2(ans_length), l = 1u << log_l, lD2 = l / 2, inv_l = mod - (modM1 >> log_l); const size_t sizeof_lD2_u = lD2 * sizeof_u; a.resize(l);
		aux1.reserve(l); aux1.resize(lD2); MemoryCopy(aux1, b, sizeof_lD2_u); Inv(lD2, aux1, aux2); b.resize(l); aux2.resize(l); DIFNTT(l, aux2); aux1.resize(l);  MemoryCopy(aux1, a, sizeof_lD2_u);
		MemoryClear(aux1 + lD2, sizeof_lD2_u); DIFNTT(l, aux1); for (u i = 0; i < l; aux1[i] = aux1[i] * SC(ull, aux2[i]) % mod, ++i); DITNTT(l, aux1);
		aux1.front() = aux1.front() * SC(ull, inv_l) % mod; for (u i = 1; i < lD2; aux1[i] = aux1[l - i] * SC(ull, inv_l) % mod, ++i); MemoryCopy(a, aux1, sizeof_lD2_u);
		MemoryClear(aux1 + lD2, sizeof_lD2_u); DIFNTT(l, aux1); DIFNTT(l, b); for (u i = 0; i < l; b[i] = b[i] * SC(ull, aux1[i]) % mod, ++i); DITNTT(l, b); MemoryClear(aux1, sizeof_lD2_u);
		for (u i = lD2; i < l; aux1[i] = (SC(ull, b[l - i]) * inv_l + mod - a[i]) % mod, ++i); DIFNTT(l, aux1); for (u i = 0; i < l; aux1[i] = aux1[i] * SC(ull, aux2[i]) % mod, ++i);
		DITNTT(l, aux1); for (u i = lD2; i < ans_length; a[i] = SC(ull, mod - aux1[l - i]) * inv_l % mod, ++i); a.resize(ans_length); return;
	}
	inline void Der(const u ans_length, const Poly& a, Poly& ans) { ans.resize(ans_length); for (u i = 1; i <= ans_length; ans[i - 1] = a[i] * SC(ull, i) % mod, ++i); return; }
	inline void Int(const u ans_length, const Poly& a, Poly& ans)
	{
		static std::vector<u> inv(2, 1); for (; inv.size() < ans_length; inv.emplace_back(SC(u, inv[mod % SC(u, inv.size())] * SC(ull, mod - mod / SC(u, inv.size())) % mod)));
		ans.resize(ans_length); ans.front() = 0; for (u i = ans_length; i > 1; --i, ans[i] = SC(ull, a[i - 1]) * inv[i] % mod); return;
	}
	inline void Int(const u ans_length, const Poly& a, Poly& ans, const std::function<u(const u)>& inv)
	{
		ans.resize(ans_length); ans.front() = 0; for (u i = ans_length; i > 1; --i, ans[i] = SC(ull, a[i - 1]) * inv(i) % mod); return;
	}
	inline void Ln(const u ans_length, Poly& a)
	{
		if (ans_length == 1) { a.resize(1); a.front() = 0; return; } static Poly aux; a.resize(ans_length); Der(ans_length - 1, a, aux); Div(ans_length - 1, aux, a); Int(ans_length, aux, a); return;
	}
	inline void Ln(const u ans_length, Poly& a, const std::function<u(const u)>& inv)
	{
		if (ans_length == 1) { a.resize(1); a.front() = 0; return; } static Poly aux; a.resize(ans_length);
		Der(ans_length - 1, a, aux); Div(ans_length - 1, aux, a); Int(ans_length, aux, a, inv); return;
	}
	inline void IntDiv(const Poly& a, const Poly& b, Poly& ans)
	{
		static Poly aux; const u sizaM1 = a.size() - 1, sizbM1 = b.size() - 1, l = a.size() - sizbM1, siz_aux = std::min(b.size(), l); ans.resize(l);
		for (u i = 0; i < l; ans[i] = a[sizaM1 - i], ++i); aux.resize(siz_aux); for (u i = 0; i < siz_aux; aux[i] = b[sizbM1 - i], ++i); Div(l, ans, aux); std::reverse(ans.data(), ans + l); return;
	}
	inline void Mod(Poly& a, Poly& b) { static Poly aux; const u l = b.size() - 1; IntDiv(a, b, aux); a.resize(l); Mul(l, b, aux); for (u i = 0; i < l; UpdateDec(a[i], b[i]), ++i); return; }
	inline void DivMod(Poly& a, Poly& b)
	{
		static Poly aux1, aux2; IntDiv(a, b, aux1); aux2.resize(aux1.size()); MemoryCopy(aux2, aux1, aux1.size() * sizeof_u); a.swap(aux1); const u l = b.size() - 1;
		b.resize(l); Mul(l, b, aux2); for (u i = 0; i < l; b[i] = UnderflowReduce(aux1[i] - b[i]), ++i); return;
	}
#undef SC
#undef CE
}
using ModularPolynomial::Poly;
void poly_multiply(unsigned *a, int n, unsigned *b, int m, unsigned *c)
{
	ModularPolynomial::Prepare(n + m + 1);
	Poly f(a, a + n + 1), g(b, b + m + 1);
	Mul(n + m + 1, f, g);
	for (unsigned i = 0; i < n + m + 1; c[i] = f[i], ++i);
	return;
}

评测结果