代码

#include <math.h>
#define N 1048576
#define pi 3.1415926535897932384626
#define op cplx operator
struct cplx
{
	double c, s;
	cplx(double _c = 0, double _s = 0) : c(_c), s(_s) {}
	op + (const cplx &w) const { return cplx(c + w.c, s + w.s); }
	op - (const cplx &w) const { return cplx(c - w.c, s - w.s); }
	op *(const cplx &w) const { return cplx(c * w.c - s * w.s, s * w.c + c * w.s); }
	op - () const { return cplx(-c, -s); }
} A[N], B[N], C[N], E[N];
typedef cplx *const cc;
int rev[N];
inline void dft(cc X)
{
	for (register int i = 1; i < N; i <<= 1)
	{
		for (register int j = 0; j < N; j += i << 1)
		{
			cc f = X + j, g = f + i, e = E + i;
			for (register int k = 0; k < i; ++k)
			{
				cplx x = f[k], y = g[k] * e[k];
				f[k] = x + y;
				g[k] = x - y;
			}
		}
	}
}
inline void cal(unsigned *a, int n, cc A)
{
	for (register int i = 0; i <= n; i += 2)
		A[rev[i >> 1]] = cplx(a[i], i == n ? 0 : a[i + 1]);
	dft(A);
	A[N] = A[0];
}
inline void trans(cplx &a0, cplx &a1, const cplx &v0, const cplx &v1)
{
	a0 = cplx(v0.c + v1.c, v0.s - v1.s);
	a1 = cplx(-v0.s - v1.s, v0.c - v1.c);
}
void poly_multiply(unsigned *a, int n, unsigned *b, int m, unsigned *c)
{
	for (register int i = 1; i < N; ++i)
		rev[i] = rev[i >> 1] >> 1 | (i & 1) << 19;
	E[1] = cplx(1, 0);
	for (register int i = 2; i < N; i <<= 1)
	{
		cc e0 = E + i / 2, e1 = E + i;
		cplx w(cos(pi / i), sin(pi / i));
		for (register int j = 0; j < i; j += 2)
			e1[j] = e0[j >> 1], e1[j + 1] = e1[j] * w;
	}
	cal(a, n, A);
	cal(b, m, B);
	cplx _i(0, 1);
	for (register int i = 0; i < N; ++i)
	{
		cplx A0, A1, B0, B1;
		trans(A0, A1, A[i], A[N - i]);
		trans(B0, B1, B[i], B[N - i]);
		cplx e = i < N / 2 ? E[N / 2 + i] : -E[i];
		C[rev[i]] = A0 * B0 + A1 * B1 * e - (A0 * B1 + A1 * B0) * _i;
	}
	for (register int i = 1; i < N; ++i)
		E[i].s *= -1;
	dft(C);
	for (register int i = 0; i <= n + m; i += 2)
	{
		cplx x = C[i >> 1] - C[N - 1];
		c[i] = int(x.c / (N * 4) + .5);
		if (i == n + m)
			break;
		c[i + 1] = int(x.s / (N * 4) + .5);
	}
}

评测结果