代码

// by newbiechd
#include <algorithm>
#include <cassert>
#include <cstdio>
#include <cstring>
#include <iostream>
#include <vector>

// Delete the debugging information!
#define debug(x) std::cerr << #x << " = " << (x) << std::endl

const int N_MAX = 2097153, mod = 998244353;
int power(int x, int y) {
  int o = 1;
  while (y) {
    if (y & 1)
      o = 1ll * o * x % mod;
    x = 1ll * x * x % mod, y >>= 1;
  }
  return o;
}

unsigned rev[N_MAX], unitRoot[N_MAX];
int getLen(int n) { return 1 << (32 - __builtin_clz(n)); }
int prepare(int n) {
  int len = getLen(n), i;
  for (i = 1; i < len; ++i)
    rev[i] = (rev[i >> 1] >> 1) | (i & 1 ? len >> 1 : 0);
  return len;
}

typedef std::vector<int> Poly;
void swap(int &x, int &y) { x ^= y, y ^= x, x ^= y; }
unsigned reduce(int x) { return x + ((x >> 31) & mod); }
void FFT(Poly &f, int len) {
  static unsigned long long tmp[N_MAX];
  for (int i = 0; i < len; ++i)
    tmp[i] = reduce(f[rev[i]]);
  for (int i = 1; i < len; i <<= 1) {
    for (int j = 0; j < len; j += i << 1)
      for (int k = 0; k < i; ++k) {
        unsigned x = tmp[i | j | k] * unitRoot[i | k] % mod;
        tmp[i | j | k] = tmp[j | k] + mod - x, tmp[j | k] += x;
      }
  }
  for (int i = 0; i < len; ++i)
    f[i] = tmp[i] % mod;
}

void operator*=(Poly &f, Poly &g) {
  const int n = f.size(), m = g.size(), k = n + m - 1, len = prepare(k);
  f.resize(len), g.resize(len), FFT(f, len), FFT(g, len);
  for (int i = 0; i < len; ++i)
    f[i] = 1ull * f[i] * g[i] % mod;
  FFT(f, len), std::reverse(&f[1], &f[len]), f.resize(k);
  const int invLen = power(len, mod - 2);
  for (int i = 0; i < k; ++i)
    f[i] = 1ll * f[i] * invLen % mod;
}

void poly_multiply(unsigned *a, int n, unsigned *b, int m, unsigned *c) {
  ++n, ++m; // n = n + 1, m = m + 1
  const int len = getLen(n + m - 1);
  for (int i = 1; i < len; i <<= 1) {
    unitRoot[i] = 1;
    unsigned temp = power(3, (mod - 1) / (i << 1));
    for (int j = 1; j < i; ++j)
      unitRoot[i | j] = 1ll * unitRoot[(i | j) - 1] * temp % mod;
  }
  
  Poly f(n), g(m);
  memcpy(&f[0], a, n << 2), memcpy(&g[0], b, m << 2), f *= g;
  memcpy(c, &f[0], (n + m - 1) << 2);
}

评测结果