代码

#include <algorithm>
#include <math.h>
using namespace std;
namespace { // radix_4_fft
struct Complex {
  typedef double value_t;
  value_t real, imag;
  Complex() = default;
  Complex(value_t a, value_t b) : real(a), imag(b) {}
  Complex(const Complex &rhs) : real(rhs.real), imag(rhs.imag) {}
  ~Complex() = default;
  Complex conj() const { return {real, -imag}; }
  Complex &operator=(const Complex &rhs) { return real = rhs.real, imag = rhs.imag, *this; }
  Complex &operator+=(const Complex &rhs) { return real += rhs.real, imag += rhs.imag, *this; }
  Complex &operator-=(const Complex &rhs) { return real -= rhs.real, imag -= rhs.imag, *this; }
  Complex &operator*=(const Complex &rhs) {
    value_t tmp1 = real * rhs.real - imag * rhs.imag, tmp2 = real * rhs.imag + imag * rhs.real;
    return real = tmp1, imag = tmp2, *this;
  }
  Complex &operator/=(const Complex &rhs) {
    value_t tmp1 =
                (real * rhs.real + imag * rhs.imag) / (rhs.real * rhs.real + rhs.imag * rhs.imag),
            tmp2 =
                (imag * rhs.real - real * rhs.imag) / (rhs.real * rhs.real + rhs.imag * rhs.imag);
    return real = tmp1, imag = tmp2, *this;
  }
  friend Complex operator+(const Complex &lhs, const Complex &rhs) { return Complex(lhs) += rhs; }
  friend Complex operator-(const Complex &lhs, const Complex &rhs) { return Complex(lhs) -= rhs; }
  friend Complex operator*(const Complex &lhs, const Complex &rhs) { return Complex(lhs) *= rhs; }
  friend Complex operator/(const Complex &lhs, const Complex &rhs) { return Complex(lhs) /= rhs; }
};
const Complex::value_t PI = acos(-1.0);
inline int getlog(int n) {
  int lgn = 0;
  while (n >>= 1) ++lgn;
  return lgn;
}
inline void revbin(int n, Complex x[]) {
  for (int i = 0, j = 0; i < n; ++i) {
    if (i < j) std::swap(x[i], x[j]);
    for (int k = n >> 1; (j ^= k) < k; k >>= 1) {
    }
  }
}
inline void dit_fft_core(int n, Complex x[]) {
  int lgn = getlog(n);
  if (lgn & 1) {
    for (int j = 0; j < n; j += 2) {
      Complex A = x[j], B = x[j + 1];
      x[j] = A + B, x[j + 1] = A - B;
    }
  }
  for (int i = 4 << (lgn & 1); i <= n; i <<= 2) {
    int l = i >> 2;
    Complex w = {cos(PI / (l << 1)), sin(-PI / (l << 1))};
    for (int j = 0; j < n; j += i) {
      Complex o1 = {1, 0}, o2 = o1, o3 = o1;
      for (int k = 0; k < l; ++k) {
        Complex A = x[j + k], B = o1 * x[j + k + 2 * l], C = o2 * x[j + k + l],
                D = o3 * x[j + k + 3 * l];
        x[j + k] = A + B + C + D;
        x[j + k + 2 * l] = A - B + C - D;
        std::swap(B.real, B.imag), std::swap(D.real, D.imag);
        B.real = -B.real, D.real = -D.real;
        x[j + k + l] = A - B - C + D;
        x[j + k + 3 * l] = A + B - C - D;
        o1 *= w, o2 = o1 * o1, o3 = o1 * o2;
      }
    }
  }
}
inline void dif_fft_core(int n, Complex x[]) {
  int lgn = getlog(n);
  for (int i = n; i >= 4 << (lgn & 1); i >>= 2) {
    int l = i >> 2;
    Complex w = {cos(PI / (l << 1)), sin(-PI / (l << 1))};
    for (int j = 0; j < n; j += i) {
      Complex o1 = {1, 0}, o2 = o1, o3 = o1;
      for (int k = 0; k < l; ++k) {
        Complex A = x[j + k], B = x[j + k + l], C = x[j + k + 2 * l], D = x[j + k + 3 * l];
        x[j + k] = A + B + C + D;
        x[j + k + l] = o2 * (A - B + C - D);
        std::swap(B.real, B.imag), std::swap(D.real, D.imag);
        B.real = -B.real, D.real = -D.real;
        x[j + k + 2 * l] = o1 * (A - B - C + D);
        x[j + k + 3 * l] = o3 * (A + B - C - D);
        o1 *= w, o2 = o1 * o1, o3 = o1 * o2;
      }
    }
  }
  if (lgn & 1) {
    for (int j = 0; j < n; j += 2) {
      Complex A = x[j], B = x[j + 1];
      x[j] = A + B, x[j + 1] = A - B;
    }
  }
}
inline void dft(int n, Complex x[]) { dif_fft_core(n, x); }
inline void idft(int n, Complex x[]) {
  for (int i = 0; i < n; ++i) x[i].imag = -x[i].imag;
  dit_fft_core(n, x);
  Complex::value_t c = static_cast<Complex::value_t>(1) / (n << 1);
  for (int i = 0; i < n; ++i) x[i].imag *= -c;
}
} // namespace
Complex ans[1U << 21U];
void poly_multiply(unsigned *a, int n, unsigned *b, int m, unsigned *c) {
  for (int i = 0; i <= n || i <= m; ++i) ans[i] = {(double)a[i], (double)b[i]};
  int len = 1;
  while (len < n + m + 1) len <<= 1;
  dft(len, ans);
  for (int i = 0; i < len; ++i) ans[i] *= ans[i];
  idft(len, ans);
  for (int i = 0; i <= n + m; ++i) c[i] = ans[i].imag + 0.5;
}

评测结果