代码

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

typedef long long int64;

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

int size;
int bitrev[Max_size];
double c[Max_size], s[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;
	c[size] = 1, s[size] = 0;
	double pc = cos(Pi / size), ps = sin(Pi / size);
	for (int i = 1; i != size; ++i)
		c[size + i] = c[size + i - 1] * pc - s[size + i - 1] * ps, s[size + i] = s[size + i - 1] * pc + c[size + i - 1] * ps;
	for (int i = size - 1; i; --i)
		c[i] = c[i << 1], s[i] = s[i << 1];
	size <<= 1;
}

inline void DHT(double 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]);
	if (1 < L)
		for (int i = 0; i != L; i += 2)
		{
			double x = A[i], y = A[i + 1];
			A[i] = x + y, A[i + 1] = x - y;
		}
	if (2 < L)
		for (int i = 0; i != L; i += 4)
		{
			double x = A[i], y = A[i + 2];
			A[i] = x + y, A[i + 2] = x - y;
			x = A[i + 1], y = A[i + 2 + 1];
			A[i + 1] = x + y, A[i + 2 + 1] = x - y;
		}
	if (4 < L)
		for (int i = 0; i != L; i += 8)
		{
			{
				double x = A[i + 5], y = A[i + 7];
				A[i + 5] = c[5] * x + s[5] * y, A[i + 7] = s[5] * x - c[5] * y;
			}
			double x = A[i], y = A[i + 4];
			A[i] = x + y, A[i + 4] = x - y;
			x = A[i + 1], y = A[i + 4 + 1];
			A[i + 1] = x + y, A[i + 4 + 1] = x - y;
			x = A[i + 2], y = A[i + 4 + 2];
			A[i + 2] = x + y, A[i + 4 + 2] = x - y;
			x = A[i + 3], y = A[i + 4 + 3];
			A[i + 3] = x + y, A[i + 4 + 3] = x - y;
		}
	for (int d = 8; d < L; d <<= 1)
		for (int i = 0; i != L; i += d << 1)
		{
			{
				int j = d + 1, k = d + d - 1;
				for (; j + 7 < k - 7; j += 8, k -= 8)
				{
					double x = A[i + j], y = A[i + k];
					A[i + j] = c[j] * x + s[j] * y, A[i + k] = s[j] * x - c[j] * y;
#define calc(offset)\
					x = A[i + j + offset], y = A[i + k - offset];\
					A[i + j + offset] = c[j + offset] * x + s[j + offset] * y, A[i + k - offset] = s[j + offset] * x - c[j + offset] * y;
					calc(1)calc(2)calc(3)calc(4)calc(5)calc(6)calc(7)
#undef calc
				}
				for (; j < k; ++j, --k)
				{
					double x = A[i + j], y = A[i + k];
					A[i + j] = c[j] * x + s[j] * y, A[i + k] = s[j] * x - c[j] * y;
				}
			} 
			for (int j = 0; j != d; j += 8)
			{
				double x = A[i + j], y = A[i + d + j];
				A[i + j] = x + y, A[i + d + j] = x - y;
#define calc(offset)\
				x = A[i + j + offset], y = A[i + d + j + offset];\
				A[i + j + offset] = x + y, A[i + d + j + offset] = x - y;
				calc(1)calc(2)calc(3)calc(4)calc(5)calc(6)calc(7)
#undef calc
			}
		}
}

int N, M, L;
double A[Max_size], B[Max_size], C[Max_size];

void poly_multiply(unsigned *a, int n, unsigned *b, int m, unsigned *c)
{
	N = n, M = m;
	for (int i = 0, x; i <= N; ++i)
		A[i] = a[i];
	for (int i = 0, x; i <= M; ++i)
		B[i] = b[i];
	
	for (L = 2; L <= N + M; L <<= 1)
		;
	init(L);
	
	DHT(A, L), DHT(B, L);
	C[0] = A[0] * B[0] / L;
	for (int i = 1; i <= L >> 1; ++i)
	{
		double x = A[i] + A[L - i], y = A[i] - A[L - i];
		C[i] = (B[i] * x + B[L - i] * y) / 2 / L, C[L - i] = (B[L - i] * x - B[i] * y) / 2 / L;
	}
	DHT(C, L);
	
	for (int i = 0; i <= N + M; ++i)
		c[i] = C[i] + 0.5;
	
	return 0;
}

评测结果