代码

#include<algorithm>
#include<cmath>

constexpr double pi = acos(-1.);

int lim, l, to[1<<21|1];

struct Complex {
  double x, y;

  Complex() {}

  Complex(double X, double Y): x(X), y(Y) {}

  Complex operator+(const Complex &rhs) const {
    return Complex(x + rhs.x, y + rhs.y);
  }

  Complex operator-(const Complex &rhs) const {
    return Complex(x - rhs.x, y - rhs.y);
  }

  Complex operator*(const Complex &rhs) const {
    return Complex(x * rhs.x - y * rhs.y, x * rhs.y + y * rhs.x);
  }
} F[1<<21|1];

void FFT(Complex *a, int type) {
  for (int i = 0; i <= lim; ++i)
    if (i < to[i]) std::swap(a[i], a[to[i]]);
  Complex rt, w, t1, t2;
  for (int m = 1; m < lim; m <<= 1) {
    rt = Complex(cos(pi / m), type * sin(pi / m));
    for (int j = 0; j < lim; j += m << 1) {
      w = Complex(1, 0);
      for (int k = 0; k < m; ++k, w = w * rt) {
        t1 = a[j | k], t2 = w * a[j | k | m];
        a[j | k] = t1 + t2, a[j | k | m] = t1 - t2;
      }
    }
  }
  if (type == -1)
    for (int i = 0; i <= lim; ++i)
      a[i].x /= lim, a[i].y /= lim;
}

void poly_multiply(unsigned *a, int N, unsigned *b, int M, unsigned *c) {
  for (lim = 1, l = -1; lim <= N + M; lim <<= 1, ++l);
  for (int i = 0; i <= lim; ++i)
    to[i] = (to[i >> 1] >> 1) | ((i & 1) << l);
  for (int i = 0, t; i <= N; ++i)
    F[i].x = a[i];
  for (int i = 0, t; i <= M; ++i)
    F[i].y = b[i];
  FFT(F, 1);
  for (int i = 0; i <= lim; ++i)
    F[i] = F[i] * F[i];
  FFT(F, -1);
  for (int i = 0; i <= N + M; ++i)
    c[i] = F[i].y / 2 + 0.5;
}

评测结果