题解：P13350 遗传

题解：P13350 遗传

这个做法不依赖树的随机性，而且是线性的。

题目描述

现在有一种和生物 M 有关的遗传病 I，已知其是显性还是隐性，其受一对常染色体等位基因 \((A,a)\) 控制，该基因的遗传符合孟德尔遗传规律，不考虑基因突变、染色体变异，交叉互换等特殊情况。

已知 \(A\) 的基因频率是 \(p\)，即认为在没有其他条件的影响下（在本题中“其他条件”包括“知道其父母”），一个生物 M 的两个基因各自有 \(p\) 的概率为 \(A\)，\(1-p\) 的概率为 \(a\)，且相互独立。

现在有 \(n\) 个生物 M，第 \(i\) 个生物的父亲为 \(fa_i\)，母亲为 \(mo_i\)，若 \(fa_i=0\) 或 \(mo_i=0\) 则表示父亲或母亲未知。为了简化问题，保证 \(fa_i\) 和 \(mo_i\) 要么均为 \(0\)，要么均不为 \(0\)，且恰有一个生物 M 没有孩子，其余生物 M 恰有一个孩子。

给出 \(n\) 个生物的关系，已知部分生物的患病情况，你需要分别求出每一个生物基因型为 \(AA,Aa,aa\) 的概率，对 \(10^9+7\) 取模。

保证数据合法且保证数据随机生成。

对于 \(100\%\) 的测试数据，保证：\(1\leq fa_i,mo_i\leq n \leq10^5\)，\(t\in\{0,1\}\)，\(1\leq a< b\leq 10^9\)，\(k_i\in\{0,1,2\}\)，保证数据随机生成。

题解

初步转化：把树倒过来，双亲改成儿子，孩子改成父亲，变成二叉树；删掉概率，改算方案数，直接认为一个独立的基因有 \(p\) 个方案是显性、\(1-p\) 个方案是隐性。

这个题的方案数计算非常扭曲，一个点的方案数与其父亲和子树都有关系，我们尝试先计算只关注它的子树填法的方案数。

假设 \(a_u\) 是每个点的显性基因个数，这个是随机变量（未定变量）。令 \(f_{u,i}\) 表示考虑了 \(u\) 子树内所有点的填法，在所有的方案中 \(a_u=i\) 的数量。有一个显然的 DP：假设 \(l,r\) 分别是 \(u\) 的两个儿子，枚举 \(i,j\)，计划用 \(f_{l,i}f_{r,j}\) 转移，再枚举 \(c_l,c_r\) 表示抽到了哪个基因，那么最终 \(a_u=[c_l\leq i]+[c_r\leq j]\)。如果这个 \(a_u\) 合法，我们就将 \(f_{l,i}f_{r,j}\) 加到 \(f_{u,a_u}\) 中。这样我们就能求出根的 DP 值了，真是一个非常正确的算法。

考虑怎么倒推回去，求出每个点全局的答案。这有一点像换根 DP，又毫无关系。我们考虑，在用 \(f_{l,i}f_{r,j}\) 向 \(f_{u,a_u}\) 转移时，我们做了两件事：

• 把 \(f_{l,i}\) 对应的所有方案复制了 \(f_{r,j}\) 次，再补充上 \(a_u\) 的具体值。\(r\) 子树的 \(a\) 值我们不关心。
• \(f_{r,j}\) 是对称的。

所以我们记录一个数组 \(bc_{l,a_u,i}\) 把这些 \(f_{r,j}\) 的和记录下来。做第二次搜索，记录一个长度为 \(3\) 的系数数组 \(coe_i\)，表示现在要把 \(f_{u,i}\) 复制 \(coe_i\) 遍。做完复制之后向下递归到儿子 \(v\)，新的 \(coe'_i\) 就设置成 \(\sum_{j=0}^2coe'_jbc_{v,j,i}\)，可以理解成把要复制的东西展开到儿子 \(v\) 这里，记录 \(a_v\) 为三种值时各要复制多少次。这样一路向下递归就做完了。复杂度 \(O(n+\log P)\)（只需要求 \(O(1)\) 次乘法逆元。另外我不知道能不能刻意构造数据使总方案数是模数的倍数，所以还是需要依赖一下数据的随机性，尤其是那个 \(p=a/b\) 的随机性。）

总结

感觉最重要的一步是放弃概率计算，转向方案数计算，毕竟前者需要条件概率知识，太扭曲了。

代码

#include <bits/stdc++.h>
using namespace std;
#ifdef LOCAL
#define debug(...) fprintf(stderr, __VA_ARGS__)
#else
#define endl "\n"
#define debug(...) void(0)
#endif
using LL = long long;
template <unsigned umod>
struct modint {/*{{{*/
  static constexpr int mod = umod;
  unsigned v;
  modint() = default;
  template <class T, enable_if_t<is_integral<T>::value, int> = 0>
    modint(const T& y) : v((unsigned)(y % mod + (is_signed<T>() && y < 0 ? mod : 0))) {}
  modint& operator+=(const modint& rhs) { v += rhs.v; if (v >= umod) v -= umod; return *this; }
  modint& operator-=(const modint& rhs) { v -= rhs.v; if (v >= umod) v += umod; return *this; }
  modint& operator*=(const modint& rhs) { v = (unsigned)(1ull * v * rhs.v % umod); return *this; }
  modint& operator/=(const modint& rhs) { assert(rhs.v); return *this *= qpow(rhs, mod - 2); }
  friend modint operator+(modint lhs, const modint& rhs) { return lhs += rhs; }
  friend modint operator-(modint lhs, const modint& rhs) { return lhs -= rhs; }
  friend modint operator*(modint lhs, const modint& rhs) { return lhs *= rhs; }
  friend modint operator/(modint lhs, const modint& rhs) { return lhs /= rhs; }
  template <class T> friend modint qpow(modint a, T b) {
    modint r = 1;
    for (assert(b >= 0); b; b >>= 1, a *= a) if (b & 1) r *= a;
    return r;
  }
  friend int raw(const modint& self) { return self.v; }
  friend ostream& operator<<(ostream& os, const modint& self) { return os << raw(self); }
  explicit operator bool() const { return v != 0; }
  modint operator-() const { return modint(0) - *this; }
  bool operator==(const modint& rhs) const { return v == rhs.v; }
  bool operator!=(const modint& rhs) const { return v != rhs.v; }
};/*}}}*/
using mint = modint<1000000007>;
constexpr int N = 1e5 + 10;
int n, t, ch[N][2], kl[N];
mint gpr;
bool vis[N];
mint f[N][3], bc[N][3][3];
bool check(int u, int c) {
  if (kl[u] == 2) return true;
  if (c == 0) return (t == 1) ^ (kl[u] == 0);
  else return (t == 0) ^ (kl[u] == 0);
}
void dfs(int u) {
  if (!ch[u][0]) {
    mint p = gpr, q = 1 - p;
    f[u][0] = q * q, f[u][1] = 2 * p * q, f[u][2] = p * p;
    for (int i = 0; i < 3; i++) if (!check(u, i)) f[u][i] = 0;
  } else {
    dfs(ch[u][0]);
    dfs(ch[u][1]);
    int l = ch[u][0], r = ch[u][1];
    for (int i = 0; i < 3; i++) {
      for (int cl = 1; cl < 3; cl++) {
        for (int j = 0; j < 3; j++) {
          for (int cr = 1; cr < 3; cr++) {
            int cu = (cl <= i) + (cr <= j);
            if (check(u, cu)) {
              bc[l][cu][i] += f[r][j];
              bc[r][cu][j] += f[l][i];
              f[u][cu] += f[l][i] * f[r][j];
            }
          }
        }
      }
    }
  }
}
int rt = 0;
void dfs2(int u, array<mint, 3> arr) {
  if (u != rt) {
    for (int i = 0; i < 3; i++) f[u][i] *= arr[i];
  }
  if (ch[u][0]) {
    for (int v : ch[u]) {
      array<mint, 3> nxt = {};
      for (int i = 0; i < 3; i++) {
        for (int j = 0; j < 3; j++) nxt[j] += arr[i] * bc[v][i][j];
      }
      dfs2(v, nxt);
    }
  }
}
int main() {
#ifndef LOCAL
  cin.tie(nullptr)->sync_with_stdio(false);  
#endif
  int x, y;
  cin >> n >> t >> x >> y, gpr = mint(x) / y;
  for (int i = 1; i <= n; i++) {
    cin >> kl[i] >> ch[i][0] >> ch[i][1];
    vis[ch[i][0]] = vis[ch[i][1]] = true;
  }
  for (int i = 1; i <= n; i++) if (!vis[i]) rt = i;
  dfs(rt);
  dfs2(rt, {1, 1, 1});
  for (int i = 1; i <= n; i++) {
    mint s = 1 / (f[i][0] + f[i][1] + f[i][2]);
    for (int j = 0; j < 3; j++) cout << f[i][2 - j] * s << " \n"[j == 2];
  }
  return 0;
}