代码

#include <cstdio>
#include <cctype>
#include <cmath>
#include <algorithm>

constexpr double Pi = 3.1415926535897932384626433832795;
constexpr int Max_size = 1 << 21 | 5;

struct Complex { double R, I; };

inline Complex operator+(const Complex &C1, const Complex &C2) { return Complex{C1.R + C2.R, C1.I + C2.I}; }
inline Complex operator-(const Complex &C1, const Complex &C2) { return Complex{C1.R - C2.R, C1.I - C2.I}; }
inline Complex operator*(const Complex &C1, const Complex &C2) { return Complex{C1.R * C2.R - C1.I * C2.I, C1.R * C2.I + C1.I * C2.R}; }
inline Complex operator~(const Complex &C) { return Complex{C.R, -C.I}; }
inline Complex &operator*=(Complex &C1, const Complex &C2) { return C1 = C1 * C2; }

int size;
int bitrev[Max_size];
Complex w[Max_size];

inline void init(int n)
{
	for (size = 1; size < n; size <<= 1)
		;
	for (int i = 0; i != size; ++i)
		bitrev[i] = bitrev[i >> 1] >> 1 | ((i & 1) * (size >> 1));
	size >>= 1;
	Complex pr = {cos(Pi / size), sin(Pi / size)};
	w[size] = {1, 0};
	for (int i = 1; i != size; ++i)
		w[size + i] = w[size + i - 1] * pr;
	for (int i = size - 1; i; --i)
		w[i] = w[i << 1];
	size <<= 1;
}

inline void DFT(Complex A[], int L)
{
	int shift = __builtin_ctz(size) - __builtin_ctz(L);
	for (int i = 0; i != L; ++i)
		if (i < bitrev[i] >> shift)
			std::swap(A[i], A[bitrev[i] >> shift]);
	for (int d = 1; d != L; d <<= 1)
		for (int i = 0; i != L; i += d << 1)
			for (int j = 0; j != d; ++j)
			{
				Complex tmp = A[i + d + j] * w[d + j];
				A[i + d + j] = A[i + j] - tmp, A[i + j] = A[i + j] + tmp;
			}
}

inline void Real_DFT(Complex A[], int L)
{
	if (L == 1)
		return;
	static Complex tmp[Max_size];
	L >>= 1;
	for (int i = 0; i != L; ++i)
		tmp[i] = {A[i << 1].R, A[i << 1 | 1].R};
	DFT(tmp, L);
	tmp[L] = tmp[0];
	for (int i = 0; i != L; ++i)
	{
		Complex t1 = (tmp[i] + ~tmp[L - i]) * Complex{0.5, 0}, t2 = (tmp[i] - ~tmp[L - i]) * Complex{0, -0.5} * w[L + i];
		A[i] = t1 + t2, A[L + i] = t1 - t2;
	}
}

inline void Real_IDFT(Complex A[], int L)
{
	if (L == 1)
		return;
	static Complex tmp[Max_size];
	A[L] = A[0];
	std::reverse(A + 1, A + L);
	L >>= 1;
	for (int i = 0; i != L; ++i)
	{
		Complex t1 = (A[i] + A[L + i]) * Complex{0.5, 0}, t2 = (A[i] - A[L + i]) * Complex{0, 0.5} * w[L + i];
		tmp[i] = t1 + t2;
	}
	DFT(tmp, L);
	for (int i = 0; i != L; ++i)
		A[i << 1] = {tmp[i].R / L, 0}, A[i << 1 | 1] = {tmp[i].I / L, 0};
}

int N, M, L;
Complex A[Max_size], B[Max_size];

void poly_multiply(unsigned *a, int n, unsigned *b, int m, unsigned *c)
{
	N = n, M = m;
	for (int i = 0; i <= N; ++i)
		A[i].R = a[i];
	for (int i = 0; i <= M; ++i)
		B[i].R = b[i];
	
	for (L = 2; L <= N + M; L <<= 1)
		;
	init(L);
	
	Real_DFT(A, L), Real_DFT(B, L);
	for (int i = 0; i != L; ++i)
		A[i] *= B[i];
	Real_IDFT(A, L);
	
	for (int i = 0; i <= N + M; ++i)
		c[i] = A[i].R + 0.5;
}

评测结果