SOJ 703 「SPC #3」布尔立方 ~ Boolean Cube

先来铺垫一个结论。

引理.

模质数 \(q\) 意义下 \(n \times m\) 的秩为 \(r\) 的矩阵的个数为 \[ q^{r(r-1)/2} (r!)_q \binom nr_q \binom mr_q \]

证明：
显然秩为 \(r\) 的 \(r \times m\) 矩阵个数为 \[ \prod_{i=0}^{r-1} (q^m - q^i) = q^{r(r-1)/2} (r!)_q \binom mr_q \]

插入剩下的数的方案数是 \[ [z^{n-r}] \prod_{i=0}^r \frac1{1-q^iz} \]

用类似 \(q\)-二项式定理证明的方法可知其等于 \(\binom nr_q\)。

---

来看这题。
考虑斯特林容斥，这样去掉两个限制的 \(i\times j\times H\) 的立方体的个数就是 \[ \prod_{k=1}^H q^{r_k(r_k-1)/2} (r_k!)_q \binom i{r_k}_q \binom j{r_k}_q \]

可以提出 \(\prod_{k=1}^H q^{r_k(r_k-1)/2} (r_k!)_q\)，那么除此之外总共就是 \[ \left(\sum_{i=1}^L (-1)^{L-i} {L\brack i} \prod_{k=1}^H \binom i{r_k}_q\right)\left(\sum_{j=1}^W (-1)^{W-j} {W \brack j} \prod_{k=1}^H \binom j{r_k}_q\right) \]

这东西不像能结合起来算的样子，所以考虑算出所有 \[ f_i = \prod_{k=1}^H \binom i{r_k}_q \]

注意到 \[ \frac{f_i}{ f_{i-1} } = \prod_{k=1}^H \frac{1-q^i}{ 1-q^{i-r_k} } \]

那么用 CZT 算出 \[ \prod_{k=1}^H (1-q^{-r_k}z) = \exp\left(-\sum_{i\ge 1} \frac{z^i}i \sum_{k=1}^H q^{-ir_k}\right) \]

之后再做一次 CZT 即可。

接下来算答案。
就是计算 \[ \newcommand\me{ \mathrm e } \sum_{i\ge 0} f_i [t^i] \me^{t\ln(1+z)} = \sum_{i\ge 0} f_i [t^i] (1+z)^t \]

转置一下可以发现就是下降幂转普通幂： \[ \sum_{i\ge 0} f_i [z^i] \me^{t\ln(1+z)} = \sum_{i\ge 0} f_i \binom ti \]

分治 NTT 即可。

代码：

#include <cstdio>
#include <vector>
#include <cstring>
#include <algorithm>

using ll = long long;

using namespace std;

const int mod = 998244353;
inline int norm(int x) {
  return x >= mod ? x - mod : x;
}
inline int reduce(int x) {
  return x < 0 ? x + mod : x;
}
inline int neg(int x) {
  return x ? mod - x : 0;
}
inline void add(int &x, int y) {
  if ((x += y - mod) < 0)
    x += mod;
}
inline void sub(int &x, int y) {
  if ((x -= y) < 0)
    x += mod;
}
inline void fam(int &x, int y, int z) {
  x = (x + (ll)y * z) % mod;
}
inline int qpow(int a, int b) {
  int ret = 1;
  for (; b; b >>= 1)
    (b & 1) && (ret = (ll)ret * a % mod),
    a = (ll)a * a % mod;
  return ret;
}

const int N = 1e5;

namespace Poly {
  const int LG = 18;
  const int N = 1 << LG + 1;
  const int G = 3;
  
  int lg2[N + 5];
  int fac[N + 5], ifac[N + 5], inv[N + 5];
  int rt[N + 5];
  
  inline void init() {
    for (int i = 2; i <= N; ++i)
      lg2[i] = lg2[i >> 1] + 1;
    rt[0] = 1, rt[1 << LG] = qpow(G, (mod - 1) >> LG + 2);
    for (int i = LG; i; --i)
      rt[1 << i - 1] = (ll)rt[1 << i] * rt[1 << i] % mod;
    for (int i = 1; i < N; ++i)
      rt[i] = (ll)rt[i & i - 1] * rt[i & -i] % mod;
    fac[0] = 1;
    for (int i = 1; i <= N; ++i)
      fac[i] = (ll)fac[i - 1] * i % mod;
    ifac[N] = qpow(fac[N], mod - 2);
    for (int i = N; i; --i)
      ifac[i - 1] = (ll)ifac[i] * i % mod;
    for (int i = 1; i <= N; ++i)
      inv[i] = (ll)ifac[i] * fac[i - 1] % mod;
  }
  
  struct poly {
    
    vector<int> a;
    inline poly(int x = 0) {
      if (x)
        a.push_back(x);
    }
    inline poly(const vector<int> &o) {
      a = o;
    }
    inline poly(const poly &o) {
      a = o.a;
    }
    
    inline int size() const {
      return a.size();
    }
    inline bool empty() const {
      return a.empty();
    }
    inline void resize(int x) {
      a.resize(x);
    }
    inline int operator[](int x) const {
      if (x < 0 || x >= size())
        return 0;
      return a[x];
    }
    inline void clear() {
      vector<int>().swap(a);
    }
    inline poly modxn(int n) const {
      if (a.empty())
        return poly();
      n = min(n, size());
      return poly(vector<int>(a.begin(), a.begin() + n));
    }
    inline poly rever() const {
      return poly(vector<int>(a.rbegin(), a.rend()));
    }
    
    inline void dif() {
      int n = size();
      for (int i = 0, len = n >> 1; len; ++i, len >>= 1)
        for (int j = 0, *w = rt; j < n; j += len << 1, ++w)
          for (int k = j; k < j + len; ++k) {
            int R = (ll)*w * a[k + len] % mod;
            a[k + len] = reduce(a[k] - R),
            add(a[k], R);
          }
    }
    inline void dit() {
      int n = size();
      for (int i = 0, len = 1; len < n; ++i, len <<= 1)
        for (int j = 0, *w = rt; j < n; j += len << 1, ++w)
          for (int k = j; k < j + len; ++k) {
            int R = norm(a[k] + a[k + len]);
            a[k + len] = (ll)*w * (a[k] - a[k + len] + mod) % mod,
            a[k] = R;
          }
      reverse(a.begin() + 1, a.end());
      for (int i = 0; i < n; ++i)
        a[i] = (ll)a[i] * inv[n] % mod;
    }
    inline void ntt(int type = 1) {
      type == 1 ? dif() : dit();
    }
    
    friend inline poly operator+(const poly &a, const poly &b) {
      vector<int> ret(max(a.size(), b.size()));
      for (int i = 0; i < ret.size(); ++i)
        ret[i] = norm(a[i] + b[i]);
      return poly(ret);
    }
    friend inline poly operator-(const poly &a, const poly &b) {
      vector<int> ret(max(a.size(), b.size()));
      for (int i = 0; i < ret.size(); ++i)
        ret[i] = reduce(a[i] - b[i]);
      return poly(ret);
    }
    friend inline poly operator*(poly a, poly b) {
      if (a.empty() || b.empty())
        return poly();
      if (a.size() < 40 || b.size() < 40) {
        if (a.size() > b.size())
          swap(a, b);
        poly ret;
        ret.resize(a.size() + b.size() - 1);
        for (int i = 0; i < ret.size(); ++i)
          for (int j = 0; j <= i && j < a.size(); ++j)
            ret.a[i] = (ret[i] + (ll)a[j] * b[i - j]) % mod;
        return ret;
      }
      int lim = 1, tot = a.size() + b.size() - 1;
      for (; lim < tot; lim <<= 1);
      a.resize(lim), b.resize(lim);
      a.ntt(), b.ntt();
      for (int i = 0; i < lim; ++i)
        a.a[i] = (ll)a[i] * b[i] % mod;
      a.ntt(-1), a.resize(tot);
      return a;
    }
    
    poly &operator+=(const poly &o) {
      resize(max(size(), o.size()));
      for (int i = 0; i < o.size(); ++i)
        add(a[i], o[i]);
      return *this;
    }
    poly &operator-=(const poly &o) {
      resize(max(size(), o.size()));
      for (int i = 0; i < o.size(); ++i)
        sub(a[i], o[i]);
      return *this;
    }
    poly &operator*=(poly o) {
      return (*this) = (*this) * o;
    }
    
    poly deriv() const {
      if (empty())
        return poly();
      vector<int> ret(size() - 1);
      for (int i = 0; i < size() - 1; ++i)
        ret[i] = (ll)(i + 1) * a[i + 1] % mod;
      return poly(ret);
    }
    poly integ() const {
      if (empty())
        return poly();
      vector<int> ret(size() + 1);
      for (int i = 0; i < size(); ++i)
        ret[i + 1] = (ll)a[i] * inv[i + 1] % mod;
      return poly(ret);
    }
    inline poly inver(int m) const {
      poly ret(qpow(a[0], mod - 2)), f, g;
      for (int k = 1; k < m;) {
        k <<= 1, f.resize(k), g.resize(k);
        for (int i = 0; i < k; ++i)
          f.a[i] = operator[](i), g.a[i] = ret[i];
        f.ntt(), g.ntt();
        for (int i = 0; i < k; ++i)
          f.a[i] = (ll)f[i] * g[i] % mod;
        f.ntt(-1);
        for (int i = 0; i < (k >> 1); ++i)
          f.a[i] = 0;
        f.ntt();
        for (int i = 0; i < k; ++i)
          f.a[i] = (ll)f[i] * g[i] % mod;
        f.ntt(-1);
        ret.resize(k);
        for (int i = (k >> 1); i < k; ++i)
          ret.a[i] = neg(f[i]);
      }
      return ret.modxn(m);
    }
    inline pair<poly, poly> div(poly o) const {
      if (size() < o.size())
        return make_pair(poly(), *this);
      poly f, g;
      f = (rever().modxn(size() - o.size() + 1) * o.rever().inver(size() - o.size() + 1))
          .modxn(size() - o.size() + 1).rever();
      g = (modxn(o.size() - 1) - o.modxn(o.size() - 1) * f.modxn(o.size() - 1)).modxn(o.size() - 1);
      return make_pair(f, g);
    }
    inline poly log(int m) const {
      return (deriv() * inver(m)).integ().modxn(m);
    }
    inline poly exp(int m) const {
      poly ret(1), iv, it, d = deriv(), itd, itd0, t1;
      if (m < 70) {
        ret.resize(m);
        for (int i = 1; i < m; ++i) {
          for (int j = 1; j <= i; ++j)
            ret.a[i] = (ret[i] + (ll)j * operator[](j) % mod * ret[i - j]) % mod;
          ret.a[i] = (ll)ret[i] * inv[i] % mod;
        }
        return ret;
      }
      for (int k = 1; k < m;) {
        k <<= 1;
        it.resize(k >> 1);
        for (int i = 0; i < (k >> 1); ++i)
          it.a[i] = ret[i];
        itd = it.deriv(), itd.resize(k >> 1);
        iv = ret.inver(k >> 1), iv.resize(k >> 1);
        it.ntt(), itd.ntt(), iv.ntt();
        for (int i = 0; i < (k >> 1); ++i)
          it.a[i] = (ll)it[i] * iv[i] % mod,
          itd.a[i] = (ll)itd[i] * iv[i] % mod;
        it.ntt(-1), itd.ntt(-1), sub(it.a[0], 1);
        for (int i = 0; i < k - 1; ++i)
          sub(itd.a[i % (k >> 1)], d[i]);
        itd0.resize((k >> 1) - 1);
        for (int i = 0; i < (k >> 1) - 1; ++i)
          itd0.a[i] = d[i];
        itd0 = (itd0 * it).modxn((k >> 1) - 1);
        t1.resize(k - 1);
        for (int i = (k >> 1) - 1; i < k - 1; ++i)
          t1.a[i] = itd[(i + (k >> 1)) % (k >> 1)];
        for (int i = k >> 1; i < k - 1; ++i)
          sub(t1.a[i], itd0[i - (k >> 1)]);
        t1 = t1.integ();
        for (int i = 0; i < (k >> 1); ++i)
          t1.a[i] = t1[i + (k >> 1)];
        for (int i = (k >> 1); i < k; ++i)
          t1.a[i] = 0;
        t1.resize(k >> 1), t1 = (t1 * ret).modxn(k >> 1), t1.resize(k);
        for (int i = (k >> 1); i < k; ++i)
          t1.a[i] = t1[i - (k >> 1)];
        for (int i = 0; i < (k >> 1); ++i)
          t1.a[i] = 0;
        ret -= t1;
      }
      return ret.modxn(m);
    }
    inline poly sqrt(int m) const {
      poly ret(1), f, g;
      for (int k = 1; k < m;) {
        k <<= 1;
        f = ret, f.resize(k >> 1);
        f.ntt();
        for (int i = 0; i < (k >> 1); ++i)
          f.a[i] = (ll)f[i] * f[i] % mod;
        f.ntt(-1);
        for (int i = 0; i < k; ++i)
          sub(f.a[i % (k >> 1)], operator[](i));
        g = (2 * ret).inver(k >> 1), f = (f * g).modxn(k >> 1), f.resize(k);
        for (int i = (k >> 1); i < k; ++i)
          f.a[i] = f[i - (k >> 1)];
        for (int i = 0; i < (k >> 1); ++i)
          f.a[i] = 0;
        ret -= f;
      }
      return ret.modxn(m);
    }
    inline poly pow(int m, int k1, int k2 = -1) const {
      if (empty())
        return poly();
      if (k2 == -1)
        k2 = k1;
      int t = 0;
      for (; t < size() && !a[t]; ++t);
      if ((ll)t * k1 >= m)
        return poly();
      poly ret;
      ret.resize(m);
      int u = qpow(a[t], mod - 2), v = qpow(a[t], k2);
      for (int i = 0; i < m - t * k1; ++i)
        ret.a[i] = (ll)operator[](i + t) * u % mod;
      ret = ret.log(m - t * k1);
      for (int i = 0; i < ret.size(); ++i)
        ret.a[i] = (ll)ret[i] * k1 % mod;
      ret = ret.exp(m - t * k1), t *= k1, ret.resize(m);
      for (int i = m - 1; i >= t; --i)
        ret.a[i] = (ll)ret[i - t] * v % mod;
      for (int i = 0; i < t; ++i)
        ret.a[i] = 0;
      return ret;
    }
  };
}
using Poly::fac;
using Poly::ifac;
using Poly::inv;
using Poly::init;
using Poly::poly;

inline int binom(int n, int m) {
  return n < m || m < 0 ? 0 : (ll)fac[n] * ifac[m] % mod * ifac[n - m] % mod;
}

inline poly czt(const poly &f, int c, int m) {
  int n = f.size(), ci = qpow(c, mod - 2);
  poly a, b, ret;
  a.resize(n), b.resize(n + m - 1), ret.resize(m);
  vector<int> cpow(n + m - 1), cipow(max(n, m));
  cpow[0] = 1;
  for (int i = 1, pw = 1; i < cpow.size(); ++i)
    cpow[i] = (ll)cpow[i - 1] * pw % mod,
    pw = (ll)pw * c % mod;
  cipow[0] = 1;
  for (int i = 1, pw = 1; i < cipow.size(); ++i)
    cipow[i] = (ll)cipow[i - 1] * pw % mod,
    pw = (ll)pw * ci % mod;
  for (int i = 0; i < n; ++i)
    a.a[i] = (ll)f[i] * cipow[i] % mod;
  for (int i = 0; i < n + m - 1; ++i)
    b.a[n + m - 2 - i] = cpow[i];
  a *= b;
  for (int i = 0; i < m; ++i)
    ret.a[i] = (ll)cipow[i] * a[n + m - 2 - i] % mod;
  return ret;
}

namespace QAnalog {
  const int q = 2;
  int qPow[N + 5];
  int n[N + 5], fac[N + 5], ifac[N + 5];
  inline void init() {
    qPow[0] = 1;
    for (int i = 1; i <= N; ++i)
      qPow[i] = (ll)qPow[i - 1] * q % mod;
    int qi = qpow(reduce(1 - q), mod - 2);
    fac[0] = 1;
    for (int i = 1; i <= N; ++i)
      n[i] = (ll)(1 - qPow[i] + mod) * qi % mod,
      fac[i] = (ll)fac[i - 1] * n[i] % mod;
    ifac[N] = qpow(fac[N], mod - 2);
    for (int i = N; i; --i)
      ifac[i - 1] = (ll)ifac[i] * n[i] % mod;
  }
  inline int binom(int n, int m) {
    return n < m || m < 0 ? 0 : (ll)fac[n] * ifac[m] % mod * ifac[n - m] % mod;
  }
}

int n, H, L[N + 5], W[N + 5], R, lim;
int r[N + 5], coe = 1;
int f[N + 5], g[N + 5], ans;

namespace TellegensPrinciple {
  #define ls (p << 1)
  #define rs (ls | 1)
  poly seg[N * 4 + 5];
  int ans[N + 5];
  inline poly mulT(poly a, poly b, int k = -1) {
    if (a.empty() || b.empty())
      return poly();
    int n = a.size(), m = b.size();
    if (k == -1)
        k = n - m + 1;
    if (k <= 0)
        return poly();
    if (k < 40 || b.size() < 40) {
      poly ret;
      ret.resize(k);
      for (int i = 0;i < k;++i)
        for (int j = 0;j < b.size();++j)
          fam(ret.a[i], a[i + j], b[j]);
      return ret;
    }
    reverse(b.a.begin(), b.a.end());
    if (k == n - m + 1) {
      int lim = 1;
      for (; lim < n; lim <<= 1);
      a.resize(lim), b.resize(lim);
      a.ntt(), b.ntt();
      for (int i = 0; i < lim; ++i)
        a.a[i] = (ll)a[i] * b[i] % mod;
      a.ntt(-1);
      for (int i = 0; i < k; ++i)
        a.a[i] = a[m - 1 + i];
      for (int i = k; i < lim; ++i)
        a.a[i] = 0;
      a.resize(n - m + 1);
    } else {
      a *= b, a.resize(n + m - 1);
      for (int i = 0; i < k; ++i)
        a.a[i] = a[m - 1 + i];
      for (int i = k; i < n + m - 1; ++i)
        a.a[i] = 0;
      a.resize(k);
    }
    return a;
  }
  void build(int p, int l, int r) {
    if (l == r) {
      seg[p] = poly({(ll)neg(l - 1) * inv[l] % mod, inv[l]});
      return ;
    }
    int mid = l + r >> 1;
    build(ls, l, mid), build(rs, mid + 1, r);
    seg[p] = seg[ls] * seg[rs];
  }
  void solve(int p, int l, int r, poly f) {
    if (l == r) {
      for (int i = 0; i < f.size(); ++i)
        fam(ans[l], seg[p][i], f[i]);
      ans[l] = (ll)ans[l] * fac[l] % mod;
      return ;
    }
    int mid = l + r >> 1;
    solve(ls, l, mid, f.modxn(mid - l + 2)), solve(rs, mid + 1, r, mulT(f, seg[ls]));
  }
  void solve(poly f) {
    int n = f.size() - 1;
    build(1, 1, n), solve(1, 1, n, f);
  }
  #undef ls
  #undef rs
}

int main() {
  init(), QAnalog::init();
  
  scanf("%d%d", &n, &H);
  for (int i = 1; i <= H; ++i)
    scanf("%d", r + i), R = max(R, r[i]),
    coe = (ll)coe * qpow(2, (ll)r[i] * (r[i] - 1) / 2 % (mod - 1)) % mod * QAnalog::fac[r[i]] % mod;
  f[R] = 1;
  for (int i = 1; i <= H; ++i)
    f[R] = (ll)f[R] * QAnalog::binom(R, r[i]) % mod;
  for (int i = 1; i <= n; ++i)
    scanf("%d%d", L + i, W + i), lim = max({lim, L[i], W[i]});
  poly temp;
  temp.resize(R + 1);
  for (int i = 1; i <= H; ++i)
    add(temp.a[r[i]], 1);
  temp = czt(temp, inv[2], H + 1), temp.a[0] = 0;
  for (int i = 1; i <= H; ++i)
    temp.a[i] = (ll)(mod - inv[i]) * temp[i] % mod;
  temp = temp.exp(H + 1);
  temp = czt(temp, 2, lim + 1);
  for (int i = R + 1; i <= lim; ++i)
    f[i] = (ll)f[i - 1] * qpow(reduce(1 - QAnalog::qPow[i]), H) % mod * qpow(temp[i], mod - 2) % mod;
  TellegensPrinciple::solve(poly(vector<int>(f, f + lim + 1)));
  copy(TellegensPrinciple::ans + 1, TellegensPrinciple::ans + lim + 1, g + 1);
  for (int i = 1; i <= n; ++i) {
    ans = (ll)coe * g[L[i]] % mod * g[W[i]] % mod;
    printf("%d\n", ans);
  }
}