代码

#include <math.h>
using namespace std;

typedef double real;

struct Complex {
  real x, y;
  Complex operator + (const Complex &b) const { return (Complex){x + b.x, y + b.y}; }
  Complex operator - (const Complex &b) const { return (Complex){x - b.x, y - b.y}; }
  Complex operator * (const Complex &b) const
  { return (Complex){x * b.x - y * b.y, x * b.y + y * b.x}; }
  Complex MulR(real c) const { return (Complex){x * c, y * c}; }
  Complex MulI(real c) const { return (Complex){-y * c, x * c}; }
  Complex Conj(void) const { return (Complex){x, -y}; }
};

const int maxn = 1 << 21;
const real pi = acos(-1);

Complex root[maxn >> 1];
int rev[maxn];

void Init(int n) {
  rev[0] = 0;
  for (int i = 1; i < n; ++i) {
    rev[i] = i & 1 ? rev[i - 1] | n >> 1 : rev[i >> 1] >> 1;
  }
  for (int i = 0; i < n >> 1; ++i) {
    real alpha = 2 * pi * i / n;
    root[i] = (Complex){cos(alpha), sin(alpha)};
  }
}

void Dft(Complex *a, int n) {
  for (int i = 0; i < n; ++i) {
    if (i < rev[i]) {
      swap(a[i], a[rev[i]]);
    }
  }
  Complex *end = root + (n >> 1);
  for (int m = 1, step = n >> 1; m < n; m <<= 1, step >>= 1) {
    for (int i = 0; i < n; i += m << 1) {
      for (Complex *l = a + i, *r = l + m, *w = root; w < end; w += step) {
        Complex t = *r * *w;
        *r = *l - t;
        *l = *l + t;
        ++l; ++r;
      }
    }
  }
}

void Idft(Complex *a, int n) {
  reverse(a + 1, a + n);
  Dft(a, n);
  real inv = 1.0 / n;
  for (int i = 0; i < n; ++i) {
    a[i] = a[i].MulR(inv);
  }
}

void Ddft(Complex *a, Complex *b, int n) {
  for (int i = 0; i < n; ++i) {
    a[i].y = b[i].x;
  }
  Dft(a, n);
  for (int i = 0; i < n; ++i) {
    b[i] = a[i ? n - i : 0].Conj();
  }
  for (int i = 0; i < n; ++i) {
    Complex t = (a[i] + b[i]).MulR(0.5);
    b[i] = (a[i] - b[i]).MulI(-0.5);
    a[i] = t;
  }
}

void Didft(Complex *a, Complex *b, int n) {
  for (int i = 0; i < n; ++i) {
    a[i] = a[i] + b[i].MulI(1);
  }
  Idft(a, n);
  for (int i = 0; i < n; ++i) {
    b[i] = (Complex){a[i].y, 0};
    a[i].y = 0;
  }
}

void Transform(double &u, double &v, double &x, double &y) {
  double a = (u * u - v * v - x * x + y * y) * 0.25;
  u = x = (u * v + x * y) * 0.5; v = -a; y = a;
}

void Convolution(Complex *a, int n) {
  Init(n);
  Dft(a, n);
  for (int i = 0; i << 1 <= n; ++i) {
    int j = i ? n - i : 0;
    Transform(a[i].x, a[i].y, a[j].x, a[j].y);
  }
  Idft(a, n);
}

Complex A[maxn];

void poly_multiply(unsigned *a, int n, unsigned *b, int m, unsigned *c)
{
  for (int i=0;i<=n+m;++i) A[i]=(Complex){(real)a[i],(real)b[i]};
  Convolution(A,1<<21);
  for (int i=0;i<=n+m;++i) c[i]=round(A[i].x);
}

评测结果