代码

#include <math.h>

template<typename _Tp>
void swap(_Tp& x1, _Tp& x2)
{
	_Tp tmp(x1);
	x1 = x2, x2 = tmp;
}

const double Pi(3.1415926535897932384626433832795);

const int Max_L((1 << 21) + 5);

struct Complex
{
    double R, I;
    Complex(const double &_R = 0.0, const double &_I = 0.0) : R(_R), I(_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*=(Complex &C1, const Complex &C2)
{
    return C1 = C1 * C2;
}

inline Complex operator~(Complex &C)
{
    return Complex(C.R, -C.I);
}

int bitrev[Max_L];
Complex w[Max_L];

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

void DFT(Complex A[], int L)
{
    for (int i(0); i != L; ++i)
        if (i < bitrev[i])
            swap(A[i], A[bitrev[i]]);
    for (int d(1), f(L >> 1); d != L; d <<= 1, f >>= 1)
        for (int i(0); i != L; i += d << 1)
        {
            Complex *l(A + i), *r(A + i + d), *cur(w);
            for (int j(0); j != d; ++l, ++r, cur += f, ++j)
            {
                Complex tmp(*r * *cur);
                *r = *l - tmp, *l = *l + tmp;
            }
        }
}

void IDFT(Complex A[], int L)
{
    for (int i(0); i != L; ++i)
        if (i < bitrev[i])
            swap(A[i], A[bitrev[i]]);
    for (int d(1), f(L >> 1); d != L; d <<= 1, f >>= 1)
        for (int i(0); i != L; i += d << 1)
        {
            Complex *l(A + i), *r(A + i + d), *cur(w + L);
            for (int j(0); j != d; ++l, ++r, cur -= f, ++j)
            {
                Complex tmp(*r * *cur);
                *r = *l - tmp, *l = *l + tmp;
            }
        }
    for (int i(0); i != L; ++i)
        A[i] *= 1.0 / L;
}

void poly_multiply(unsigned *a, int n, unsigned *b, int m, unsigned *c)
{
	int L;
	for (L = 1; L <= n + m; L <<= 1)
		;
	init(L);
	static Complex A[Max_L], B[Max_L];
	for (int i(0); i <= n; ++i)
		A[i] = a[i];
	for (int i(0); i <= m; ++i)
		B[i] = b[i];
	DFT(A, L), DFT(B, L);
	for (int i(0); i != L; ++i)
		A[i] *= B[i];
	IDFT(A, L);
	for (int i(0); i <= n + m; ++i)
		c[i] = A[i].R + 0.5;
}

评测结果