CF2029H Message Spread

把集合随时间的变化写出来 \(S_1\to S_2\to\cdots\to S_k\)，考虑把答案拆到每一段处理。即设 \(f_S\) 为 \(\left\{1\right\}\to S\) 的概率，\(g_S\) 为 \(S\) 发生变化的期望步数。那么答案就是 \(\sum f_Sg_S\)。

\(g_S\) 是容易算的，考虑 \(f_S\) 怎么求。

一个暴力的想法是枚举新增的集合 \(T\)，\(f_S\times\dfrac{\prod\limits_{i\in T}(1-\prod\limits_{j\in S}(1-e_{i,j}))\times \prod\limits_{i\in S}\prod\limits_{j\notin S\cup T}(1-e_{i,j})}{1-\prod\limits_{i\in S}\prod\limits_{j\notin S}(1-e_{i,j})}\to f_{S\cup T}\)。但是这个式子和 \(S\) 和 \(T\) 强相关，并不好处理。

实际上难做的就是分子上左边的部分，处理的是 \(T\) 中每个点都和 \(S\) 有连边。考虑容斥，记 \(h_S\) 为当前集合在 \(S\) 之内的概率，那么左边那一块就全部消掉了，即 \(f_S\times\dfrac{ \prod\limits_{i\in S}\prod\limits_{j\notin S\cup T}(1-e_{i,j})-\prod\limits_{i\in S}\prod\limits_{j\notin S}(1-e_{i,j})}{1-\prod\limits_{i\in S}\prod\limits_{j\notin S}(1-e_{i,j})}\to h_{S\cup T}\)。

这个系数实际上就只和 \(S,T,S\cup T\) 有关了，于是我们要解决的就是类似半在线的子集卷积。这个考虑枚举 \(|S|\)，每次把 \(h\) 做一次高维差分，算出 \(f_S\) 再前缀和回去。这个转移直接用子集卷积实现就是 \(2^nn^3\)。

考虑优化。事实上可以发现，和 \(T\) 有关的项可以预处理 FWT 之后的结果，并且每次做子集卷积时，\(S\) 的部分只有当前枚举的 \(\text{popcount}\) 可能不是全 \(0\)，因此只要一次 FWT。此时瓶颈就是对于每一位做 IFWT。

进一步地，我们用到的 \(h\) 也只是当前枚举 \(\text{popcount}\) 上的结果，所以可以把若干次 FWT 之后的 \(h\) 加在一起，枚举到了再进行 IFWT，这样 IFWT 也只会执行 \(n\) 次。总复杂度就做到了 \(2^nn^2\)。

#include <bits/stdc++.h>
using namespace std;
using ll = long long;

const int mod = 998244353;
void Add(int &x, ll y) {
  x = (x + y) % mod;
}
int Pow(int x, int y) {
  int b = x, r = 1;
  for(; y; b = (ll)b * b % mod, y /= 2) {
    if(y & 1) r = (ll)r * b % mod;
  }
  return r;
}

const int kN = 22, kS = (1 << 21) + 5;
int n, m;
int e[kN][kN];
int prd[kS], inv[kS];
int f[kS], g[kS], h[kS], hh[kN][kS];
int A[kS], B[kS], C[kS];
int fwt_B[kN][kS];
vector<int> vec[kN];

void FWT_OR(int n, int *f) {
  for(int i = 0; i < n; i++) {
    for(int j = 0; j < (1 << n); j += (2 << i)) {
      int p = j, q = j + (1 << i);
      for(; p < j + (1 << i); p++, q++) {
        if((f[q] += f[p]) >= mod) f[q] -= mod;
      }
    }
  }
}
void IFWT_OR(int n, int *f) {
  for(int i = 0; i < n; i++) {
    for(int j = 0; j < (1 << n); j += (2 << i)) {
      int p = j, q = j + (1 << i);
      for(; p < j + (1 << i); p++, q++) {
        if((f[q] -= f[p]) < 0) f[q] += mod;
      }
    }
  }
}

void Convolution(int n, int pc, int *A, int *B, int res[][kS]) {
  static int fwt_A[kS];
  memcpy(fwt_A, A, (1 << n) * sizeof(int));
  FWT_OR(n, fwt_A);
  for(int i = 0; i + pc <= n; i++) {
    for(int j = 0; j < (1 << n); j++) {
      Add(res[i + pc][j], (ll)fwt_A[j] * fwt_B[i][j]);
    }
  }
}

int main() {
  // freopen("1.in", "r", stdin);
  // freopen("1.out", "w", stdout);
  ios::sync_with_stdio(0), cin.tie(0);
  cin >> n >> m;
  fill_n(prd, 1 << n, 1);
  for(int i = 1; i <= m; i++) {
    int u, v, p, q;
    cin >> u >> v >> p >> q;
    u--, v--;
    p = (ll)p * Pow(q, mod - 2) % mod;
    e[u][v] = e[v][u] = p;
    prd[(1 << u) | (1 << v)] = mod + 1 - p;
  }
  for(int i = 0; i < n; i++) {
    for(int j = 0; j < (1 << n); j += (2 << i)) {
      int p = j, q = j + (1 << i);
      for(; p < j + (1 << i); p++, q++) {
        prd[q] = (ll)prd[q] * prd[p] % mod;
      }
    }
  }
  for(int i = 0; i < (1 << n); i++) {
    inv[i] = Pow(prd[i], mod - 2);
  }
  int U = (1 << n) - 1;
  for(int i = 1; i < (1 << n); i += 2) {
    g[i] = Pow(mod + 1 - (ll)prd[U] * inv[i] % mod * inv[U ^ i] % mod, mod - 2);
  }
  for(int i = 1; i < (1 << n); i += 2) {
    h[i] = 1;
    vec[__builtin_popcount(i)].push_back(i);
  }
  for(int i = 1; i < (1 << n - 1); i++) {
    fwt_B[__builtin_popcount(i)][i] = prd[i << 1];
  }
  for(int i = 0; i < n; i++) {
    FWT_OR(n - 1, fwt_B[i]);
  }
  int res = 0;
  for(int pc = 1; pc < n; pc++) {
    IFWT_OR(n - 1, hh[pc - 1]);
    for(int s : vec[pc]) {
      Add(h[s], (ll)hh[pc - 1][s >> 1] * inv[s]);
    }
    IFWT_OR(n, h);
    for(int s : vec[pc]) {
      int no_out = (ll)prd[U] * inv[s] % mod * inv[U ^ s] % mod;
      Add(res, (ll)h[s] * g[s]);
      f[s] = (ll)h[s] * Pow(mod + 1 - no_out, mod - 2) % mod;
      Add(h[s], (ll)f[s] * (mod - no_out));
    }
    FWT_OR(n, h);
    memset(A, 0, sizeof(A));
    for(int s : vec[pc]) {
      A[s >> 1] = (ll)f[s] * prd[U] % mod * inv[U ^ s] % mod;
      Add(h[s], (ll)A[s >> 1] * inv[s]);
    }
    Convolution(n - 1, pc - 1, A, B, hh);
  }
  cout << res << "\n";
  return 0;
}