[SDOI 2014] 重建(矩阵树定理推广) | 错题本

    技术2025-10-25  10

    文章目录

    题目分析代码

    题目

    [SDOI2014] 重建

    分析

    问题要求 ∑ T ( ∏ e ∈ T p e ∏ e ∉ T ( 1 − p e ) ) \sum_{T} \left( \prod_{e \in T} p_e \prod_{e \notin T} (1 - p_e)\right) TeTpee/T(1pe) 由于矩阵树只能处理 e ∈ T e \in T eT,所以把式子变成 ∑ T ( ∏ e ∈ T p e ∏ e ( 1 − p e ) ∏ e ∈ T ( 1 − p e ) ) = ∏ e ( 1 − p e ) ∑ T ( ∏ e ∈ T p e 1 − p e ) \begin{aligned} & \sum_{T} \left( \prod_{e \in T} p_e \frac{\prod_{e} (1 - p_e)}{\prod_{e \in T} (1 - p_e)} \right) \\ =& \prod_{e} (1 - p_e) \sum_{T} \left( \prod_{e \in T} \frac{p_e }{1 - p_e} \right)\end{aligned} =T(eTpeeT(1pe)e(1pe))e(1pe)T(eT1pepe) 于是用矩阵树定理推广计算 ∑ T ( ∏ e ∈ T p e 1 − p e ) \sum_{T} \left( \prod_{e \in T} \frac{p_e }{1 - p_e} \right) T(eT1pepe),最后乘上 ∏ e ( 1 − p e ) \prod_{e} (1 - p_e) e(1pe) 即为答案。 注意 p e = 1 p_e = 1 pe=1 会除零而出问题,给 p e p_e pe 减个 ϵ \epsilon ϵ 即可。


    矩阵树定理推广:计算所有生成树边权乘积的和,即 ∑ T ( ∏ e ∈ T x e ) \sum_{T} \left( \prod_{e \in T} x_e \right) T(eTxe) 其中 x e x_e xe 为一个有理数。

    x e = p e t x_e = \dfrac{p_e}{t} xe=tpe,其中 p e , t ∈ Z p_e, t \in \Z pe,tZ,那么 ∑ T ( ∏ e ∈ T x e ) = 1 t n − 1 ∑ T ( ∏ e ∈ T p e ) \begin{aligned} & \sum_{T} \left( \prod_{e \in T} x_e \right) \\ =& \frac{1}{t^{n - 1}} \sum_{T} \left( \prod_{e \in T} p_e \right)\end{aligned} =T(eTxe)tn11T(eTpe) 当边权 p p p为整数的时候很好计算:把 u , v u, v u,v 之间连上 p p p 条边,然后计算生成树的个数就是原图所有生成树边权乘积和。

    但是 t n − 1 t^{n - 1} tn1 很显然精度不够,所以我们用 行列式的可乘性 把 t n − 1 t^{n - 1} tn1 分配给行列式中每个元素,于是我们有了一个新的基尔霍夫矩阵:度为连出去所有边的边权和,边权为原边权的相反数。然后正常矩阵树定理得到的答案就是所求的了。

    代码

    #include <bits/stdc++.h> const int MAXN = 50; const double eps = 1e-5; int N; double Mat[MAXN + 5][MAXN + 5]; void AddEdge(int u, int v, double w) { Mat[u][u] += w, Mat[v][v] += w; Mat[u][v] -= w, Mat[v][u] -= w; } inline double Abs(const double &x) { return x < 0 ? -x : x; } inline int Comp(const double &x, const double &y) { if (Abs(x - y) < eps) return 0; return x > y ? 1 : -1; } double Det(int n) { n--; double ret = 1; for (int i = 1; i <= n; i++) { int p = i; for (int j = i; j <= n; j++) if (Comp(Mat[p][i], Mat[j][i]) < 0) p = j; if (p != i) std::swap(Mat[i], Mat[p]); if (!Comp(Mat[i][i], 0)) return 0; ret *= Mat[i][i]; for (int j = i + 1; j <= n; j++) { double tmp = Mat[j][i] / Mat[i][i]; for (int k = i; k <= n; k++) Mat[j][k] -= Mat[i][k] * tmp; } } return Abs(ret); } int main() { scanf("%d", &N); double All = 1.0; for (int i = 1; i <= N; i++) for (int j = 1; j <= N; j++) { double p; scanf("%lf", &p); if (i < j) { if (Comp(p, 1) == 0) p -= eps; AddEdge(i, j, p / (1 - p)); All *= (1 - p); } } printf("%.5f", Det(N) * All); }
    Processed: 0.010, SQL: 9