代码

#include <algorithm>
#include <cstring>

typedef unsigned int uint;
typedef long long unsigned int uint64;

constexpr uint Max_size = 1 << 21 | 5;
constexpr uint g = 3, Mod = 998244353;

inline uint norm_2(const uint x)
{
	return x < Mod * 2 ? x : x - Mod * 2;
}

inline uint norm(const uint x)
{
	return x < Mod ? x : x - Mod;
}

inline uint Mod_mult(const uint x, const uint y) // 32-bit computing
{
	uint res = 0;
	__asm__(
		"mul\t%[y]\n\t"
		"div\t%[Mod]"
		: "=d"(res)
		: "a"(x), [y]"r"(y), [Mod]"r"(Mod)
	);
	return res;
}

struct Z
{
	uint v;
	Z() { }
	Z(const uint _v) : v(_v) { }
};

inline Z operator+(const Z &x1, const Z &x2) { return x1.v + x2.v < Mod ? x1.v + x2.v : x1.v + x2.v - Mod; }
inline Z operator-(const Z &x1, const Z &x2) { return x1.v >= x2.v ? x1.v - x2.v : x1.v + Mod - x2.v; }
inline Z operator*(const Z &x1, const Z &x2) { return Mod_mult(x1.v, x2.v); }
inline Z &operator*=(Z &x1, const Z &x2) { x1.v = Mod_mult(x1.v, x2.v); return x1; }

Z Power(Z Base, int Exp)
{
	Z res = 1;
	for (; Exp; Base *= Base, Exp >>= 1)
		if (Exp & 1)
			res *= Base;
	return res;
}

inline uint mf(uint x)
{
	return (static_cast<uint64>(x) << 32) / Mod;
}

int size;
uint w[Max_size], w_[Max_size];

inline void init(int n)
{
	for (size = 2; size < n; size <<= 1)
		;
	Z pr = Power(g, (Mod - 1) / size);
	size >>= 1;
	w[size] = 1, w_[size] = (static_cast<uint64>(w[size]) << 32) / Mod;
	for (int i = 1; i < size; ++i)
		w[size + i] = (w[size + i - 1] * pr).v, w_[size + i] = (static_cast<uint64>(w[size + i]) << 32) / Mod;
	for (int i = size - 1; i; --i)
		w[i] = w[i * 2], w_[i] = w_[i * 2];
	size <<= 1;
}

inline uint mult_w_2(const uint x, const int i)
{
	uint q = static_cast<uint64>(x) * w_[i] >> 32;
	return x * w[i] - q * Mod;
}

inline uint mult_w(const uint x, const int i)
{
	return norm(mult_w_2(x, i));
}

inline void DFT_fr_2(Z _A[], const int L)
{
	uint *A = reinterpret_cast<uint *>(_A);
	for (int d = L >> 1; d; d >>= 1)
		for (int i = 0; i != L; i += d << 1)
			for (int j = 0; j != d; ++j)
			{
				uint x = norm_2(A[i + j] + A[i + d + j]);
				uint y = mult_w_2(A[i + j] + Mod * 2 - A[i + d + j], d + j);
				A[i + j] = x, A[i + d + j] = y;
			}
}

inline void IDFT_fr(Z _A[], const int L)
{
	uint *A = reinterpret_cast<uint *>(_A);
	for (int d = 1; d != L; d <<= 1)
		for (int i = 0; i != L; i += d << 1)
			for (int j = 0; j != d; ++j)
			{
				uint x = norm_2(A[i + j]);
				uint t = mult_w_2(A[i + d + j], d + j);
				A[i + j] = x + t, A[i + d + j] = x + Mod * 2 - t;
			}
	std::reverse(A + 1, A + L);
	if (L == 2)
		A[0] = norm_2(A[0]), A[1] = norm_2(A[1]);
	int k = __builtin_ctz(L);
	for (int i = 0; i != L; ++i)
	{
		uint64 m = -A[i] & (L - 1);
		A[i] = norm((A[i] + m * Mod) >> k);
	}
}
Z A[Max_size], B[Max_size];

void poly_multiply(uint _A[], int N, uint _B[], int M, uint _C[])
{
	memcpy(reinterpret_cast<uint *>(A), _A, (N + 1) * sizeof(uint));
	memcpy(reinterpret_cast<uint *>(B), _B, (M + 1) * sizeof(uint));
	
	int L;
	for (L = 2; L <= N + M; L <<= 1)
		;
	init(L);
	
	DFT_fr_2(A, L), DFT_fr_2(B, L);
	for (int i = 0; i != L; ++i)
		A[i] *= B[i];
	IDFT_fr(A, L);
	
	memcpy(_C, reinterpret_cast<uint *>(A), (N + M + 1) * sizeof(uint));
}

评测结果