代码

#define NDEBUG
#include <algorithm>
#include <cassert>
#include <cstring>
#include <functional>
#include <initializer_list>
#include <iostream>
#include <memory>
#include <random>
#include <vector>

namespace field998244353 {

using i32 = int32_t;
using u32 = uint32_t;
using i64 = int64_t;
using u64 = uint64_t;

static constexpr u32 P = 998244353;

static u32 get_r(u32 mod) {
  u32 iv = mod;
  for (u32 i = 0; i != 4; ++i) iv *= 2 - mod * iv;
  return iv;
}

static const u32 r = get_r(P), r2 = -u64(P) % P;

// from https://min-25.hatenablog.com/entry/2017/08/20/171214
struct MontgomeryModInt32 {
public:
  using i32 = int32_t;
  using u32 = uint32_t;
  using i64 = int64_t;
  using u64 = uint64_t;

  MontgomeryModInt32() = default;
  ~MontgomeryModInt32() = default;
  MontgomeryModInt32(u32 v) : v(reduce(u64(v) * r2)) {}
  MontgomeryModInt32(const MontgomeryModInt32 &rhs) : v(rhs.v) {}
  static u32 reduce(u64 x) {
    u64 y = u32(x >> 32) - u32((u64(u32(x) * r) * P) >> 32);
    return y + (P & -(i32(y) < 0));
  }
  u32 get() const { return reduce(v); }
  MontgomeryModInt32 &operator=(const MontgomeryModInt32 &rhs) { return v = rhs.v, *this; }
  MontgomeryModInt32 operator-() const {
    MontgomeryModInt32 res;
    res.v = (P & -(v != 0)) - v;
    return res;
  }
  MontgomeryModInt32 &operator+=(const MontgomeryModInt32 &rhs) {
    v += rhs.v - P, v += P & -(i32(v) < 0);
    return *this;
  }
  MontgomeryModInt32 &operator-=(const MontgomeryModInt32 &rhs) {
    v -= rhs.v, v += P & -(i32(v) < 0);
    return *this;
  }
  MontgomeryModInt32 &operator*=(const MontgomeryModInt32 &rhs) {
    return v = reduce(u64(v) * rhs.v), *this;
  }
  friend MontgomeryModInt32 operator+(const MontgomeryModInt32 &lhs,
                                      const MontgomeryModInt32 &rhs) {
    return MontgomeryModInt32(lhs) += rhs;
  }
  friend MontgomeryModInt32 operator-(const MontgomeryModInt32 &lhs,
                                      const MontgomeryModInt32 &rhs) {
    return MontgomeryModInt32(lhs) -= rhs;
  }
  friend MontgomeryModInt32 operator*(const MontgomeryModInt32 &lhs,
                                      const MontgomeryModInt32 &rhs) {
    return MontgomeryModInt32(lhs) *= rhs;
  }
  friend std::istream &operator>>(std::istream &is, MontgomeryModInt32 &rhs) {
    return is >> rhs.v, rhs.v = reduce(u64(rhs.v) * r2), is;
  }
  friend std::ostream &operator<<(std::ostream &os, const MontgomeryModInt32 &rhs) {
    return os << rhs.get();
  }
  MontgomeryModInt32 pow(i64) const;

private:
  u32 v;
};

MontgomeryModInt32 MontgomeryModInt32::pow(i64 y) const {
  if ((y %= P - 1) < 0) y += P - 1; // phi(P) = P - 1
  MontgomeryModInt32 res(1), x(*this);
  for (; y != 0; y >>= 1, x *= x)
    if (y & 1) res *= x;
  return res;
}

namespace number_theoretic_transform {

u32 get_len(u32 n) { // if n=0, boom
  --n;
  n |= n >> 1, n |= n >> 2, n |= n >> 4, n |= n >> 8, n |= n >> 16;
  return ++n;
}

static MontgomeryModInt32 ROOT[1 << 21], IROOT[1 << 21];

void init_root_of_unity(u32 n) {
  static u32 lim = 1;
  static MontgomeryModInt32 G(3);
  if (lim < n) {
    const u32 l = n >> 1;
    MontgomeryModInt32 g = G.pow((P - 1) / n) /*, ig = G.pow(-i32((P - 1) / n))*/;
    ROOT[l] /*= IROOT[l]*/ = 1;
    for (u32 i = l + 1; i < n; ++i) ROOT[i] = ROOT[i - 1] * g /*, IROOT[i] = IROOT[i - 1] * ig*/;
    for (u32 i = l - 1; i >= lim; --i) ROOT[i] = ROOT[i << 1] /*, IROOT[i] = IROOT[i << 1]*/;
    lim = n;
  }
}

void radix_2_dif_ntt(u32 n, MontgomeryModInt32 x[], MontgomeryModInt32 root[]) {
  for (u32 i = n; i >= 4; i >>= 1)
    for (u32 j = 0, l = i >> 1; j != n; j += i)
      for (u32 k = 0; k != l; ++k) {
        MontgomeryModInt32 u = x[j + k], v = x[j + k + l];
        x[j + k] = u + v, x[j + k + l] = (u - v) * root[l + k];
      }
  if (n > 1)
    for (u32 i = 0; i != n; i += 2) {
      MontgomeryModInt32 u = x[i], v = x[i + 1];
      x[i] = u + v, x[i + 1] = u - v;
    }
}

void radix_2_dit_ntt(u32 n, MontgomeryModInt32 x[], MontgomeryModInt32 root[]) {
  if (n > 1)
    for (u32 i = 0; i != n; i += 2) {
      MontgomeryModInt32 u = x[i], v = x[i + 1];
      x[i] = u + v, x[i + 1] = u - v;
    }
  for (u32 i = 4; i <= n; i <<= 1)
    for (u32 j = 0, l = i >> 1; j != n; j += i)
      for (u32 k = 0; k != l; ++k) {
        MontgomeryModInt32 u = x[j + k], v = root[l + k] * x[j + k + l];
        x[j + k] = u + v, x[j + k + l] = u - v;
      }
}

void dft(u32 n, MontgomeryModInt32 x[]) { init_root_of_unity(n), radix_2_dif_ntt(n, x, ROOT); }

void idft(u32 n, MontgomeryModInt32 x[]) {
  radix_2_dit_ntt(n, x, ROOT), std::reverse(x + 1, x + n);
  MontgomeryModInt32 iv(P - (P - 1) / n);
  for (u32 i = 0; i != n; ++i) x[i] *= iv;
}

} // namespace number_theoretic_transform

using namespace number_theoretic_transform;

} // namespace field998244353

using namespace field998244353;

void poly_multiply(unsigned *a, int n, unsigned *b, int m, unsigned *c) {
  static MontgomeryModInt32 d[1 << 21], e[1 << 21];
  for (int i = 0; i <= n; ++i) d[i] = a[i];
  for (int i = 0; i <= m; ++i) e[i] = b[i];
  int len = get_len(n + m + 1);
  dft(len, d), dft(len, e);
  for (int i = 0; i != len; ++i) d[i] *= e[i];
  idft(len, d);
  for (int i = 0; i <= n + m; ++i) c[i] = d[i].get();
}

评测结果