代码

#include <bits/stdc++.h>
using namespace std;

const int N = 3000000;
struct cpx {double i, j;} w[N];
cpx operator + (cpx a, cpx b) {return (cpx){a.i + b.i, a.j + b.j};}
cpx operator - (cpx a, cpx b) {return (cpx){a.i - b.i, a.j - b.j};}
cpx operator * (cpx a, cpx b) {return (cpx){a.i * b.i - a.j * b.j, a.i * b.j + a.j * b.i};}
cpx operator / (cpx a, cpx b) {double d = b.i * b.i + b.j * b.j; return (cpx){(a.i * b.i + a.j * b.j) / d, (a.j * b.i - a.i * b.j) / d};}
cpx operator ! (cpx a) {return (cpx){a.i, -a.j};}
double pi = acos(-1);
void FFT(cpx A[], int n, int f)
{
    for (int i = 1, j = (n >> 1); i < n; ++ i)
    {
        if (i < j) swap(A[i], A[j]); int k;
        for (k = (n >> 1); j & k; j ^= k, k >>= 1); j ^= k;
    }
    for (int i = 2; i <= n; i <<= 1)
    {
        int d = n / i * f;
        for (int j = 0; j < n; j += i)
        {
            int p = (f > 0 ? 0 : n);
            for (int k = j; k < j + (i >> 1); ++ k)
            {
                cpx u = A[k], v = A[k + (i >> 1)] * w[p];
                A[k] = u + v; A[k + (i >> 1)] = u - v; p += d;
            }
        }
    }
    if (f < 0) for (int i = 0; i < n; ++ i) A[i] = A[i] / (cpx){(double)n, 0};
}
void poly_multiply(unsigned *A, int n, unsigned *B, int m, unsigned *C)
{
    static cpx tmp0[N], tmp1[N];
    static cpx tmp00[N], tmp01[N], tmp10[N], tmp11[N];

    int nn = 1; while (nn <= max(n, m)) nn <<= 1;
    w[0] = (cpx){1, 0}; w[1] = (cpx){cos(2 * pi / nn), sin(2 * pi / nn)};
    for (int i = 2; i <= nn; ++ i) w[i] = w[i - 1] * w[1];

    for (int i = 0; i < nn; ++ i) tmp0[i] = (cpx){0, 0};
    for (int i = 0; i <= n; i += 2) tmp0[i / 2].i = A[i];
    for (int i = 1; i <= n; i += 2) tmp0[i / 2].j = A[i];
    FFT(tmp0, nn, 1);
    for (int i = 0; i < nn; ++ i) tmp00[i] = (tmp0[i] + !(tmp0[i == 0 ? 0 : nn - i])) / (cpx){2, 0};
    for (int i = 0; i < nn; ++ i) tmp01[i] = (tmp0[i] - !(tmp0[i == 0 ? 0 : nn - i])) / (cpx){2, 0} * (cpx){0, -1};

    for (int i = 0; i < nn; ++ i) tmp1[i] = (cpx){0, 0};
    for (int i = 0; i <= m; i += 2) tmp1[i / 2].i = B[i];
    for (int i = 1; i <= m; i += 2) tmp1[i / 2].j = B[i];
    FFT(tmp1, nn, 1);
    for (int i = 0; i < nn; ++ i) tmp10[i] = (tmp1[i] + !(tmp1[i == 0 ? 0 : nn - i])) / (cpx){2, 0};
    for (int i = 0; i < nn; ++ i) tmp11[i] = (tmp1[i] - !(tmp1[i == 0 ? 0 : nn - i])) / (cpx){2, 0} * (cpx){0, -1};

    for (int i = 0; i < nn; ++ i) tmp0[i] = tmp00[i] * tmp10[i] + tmp01[i] * tmp11[i] * w[i];
    for (int i = 0; i < nn; ++ i) tmp1[i] = tmp00[i] * tmp11[i] + tmp01[i] * tmp10[i];

    for (int i = 0; i < nn; ++ i) tmp0[i] = tmp0[i] + tmp1[i] / (cpx){0, -1};
    FFT(tmp0, nn, -1);
    for (int i = 0; i <= n + m; i += 2)
    {
        C[i] = (tmp0[i / 2].i + 0.5);
        C[i + 1] = (tmp0[i / 2].j + 0.5);
    }
}

评测结果