Lilliput Steps

小さな一歩から着実に. 数学やプログラミングのことを書きます.

AOJ 2345 : Network Reliability

AOJ 動的計画法 グラフ 数学

問題文 : Network Reliability | Aizu Online Judge

概要

無向グラフG = (V, E)と整数 p が与えられる．Gの各枝が独立にpパーセントの確率で消失するとき，グラフが連結である確率を求めよ．

制約

  • 1 \leq |V| \leq 14
  • 0 \leq |E| \leq 100

コメント

制約からbitDPを思い浮かべるところまではすぐなんだけど，部分問題に落とすところで非常に苦戦してしまったのでメモ．

解法

f(S) := 部分グラフ G' = (S, E(S)) が連結である確率とおく．求めたいのはf(V)である．
しかしこのままではあまりいい性質がないので，f(V)と一致する別の量を考えてみる．

グラフが連結であるとは，ある頂点v \in Vを固定したときに，vからV\setminus vの任意の頂点に到達可能，すなわちvの属する連結成分の頂点集合がVに一致することである．
すなわち，

g(S) := 部分グラフ G' = (S, E(S)) で固定された頂点vの連結成分の頂点集合がSになる確率

とおくと，f(V) = g(V)となる．

この言い換えを行うことで，

\begin{align*} f(V) = g(V) &= 1 - \mathrm{P}(G \text{が非連結})\\ &= 1 - \sum_{V' \subset V} \mathrm{P}(v \text{の連結成分の頂点集合が} V' \text{と一致する})\\ &= 1 - \sum_{V' \subset V} g(V') \times \mathrm{P}(V' \text{と} V \setminus V' \text{を繋ぐ枝がすべて消失している})\\ &= 1 - \sum_{V' \subset V} g(V') \times (p/100)^{\Delta_{V', V \setminus V'}} \end{align*}

と変形できる．ここで{\Delta_{V', V \setminus V'}}は頂点集合V'と頂点集合V \setminus V'の間の枝の総数である．
よって，各状態について部分集合全体の和をとればよい．nビットの2進数で，kビットが1,n-kビットが0であるものは\begin{pmatrix}n \\ k\end{pmatrix}通りあり，そのビットを集合と見たときに部分集合は2^k通りあるので，状態遷移の回数は

\displaystyle\sum_{i=0}^n \begin{pmatrix} n \\ i\end{pmatrix} 2^i = (1 + 2)^n = 3^n

となる．よって\mathrm{O}(3^{|V|})時間でこの問題が解けた．

コード

#include <bits/stdc++.h>
 
using namespace std;
 
int n, m;
int cnt[14][14];
int cntS[1 << 14];
double memo[1 << 14], p;
 
// P(V : univ and equiv_class(0) == V) = 1 - P(V : univ and equiv_class(0) != V)
//                  = 1 - Σ_{V' : proper subset of V} P(V : univ and equiv_class(0) == V')
//                  = 1 - Σ_{V' : proper subset of V} P(V \ V' and V' is disconnected and V' : univ and equiv_class(0) == V')
double calc(int bit)
{
    if (bit == 1) return 1;
    if (memo[bit] >= 0) return memo[bit];
 
    double ret = 1;
    for (int i = bit; i; i = (bit & (i - 1))) {
        if (i == bit || i % 2 == 0) continue;
        ret -= pow(p, cntS[bit] - cntS[i] - cntS[bit & ~i]) * calc(i);
    }
    return memo[bit] = ret;
}
 
int main()
{
    scanf("%d %d %lf", &n, &m, &p);
    p /= 100;
 
    for (int i = 0; i < m; i++) {
        int a, b;
        scanf("%d %d", &a, &b);
        cnt[a - 1][b - 1]++;
        cnt[b - 1][a - 1]++;
    }
 
    for (int i = 0; i < (1 << n); i++)
        memo[i] = -1;
 
    for (int i = 0; i < (1 << n); i++) {
        for (int j = 0; j < n; j++) {
            for (int k = j; k < n; k++) {
                if (((i >> j) & 1) && ((i >> k) & 1)) cntS[i] += cnt[j][k];
            }
        }
    }
 
    printf("%.10f\n", calc((1 << n) - 1));
 
    return 0;
}