WBC004 E - Link 解説

2021年3月18日に開かれたウインドベル(@Twi_Stamp)主催の競プロ有志コンテストWind Bel Contest 004 – Good Bye 2A Day 2の解説です。

問題文

時刻 $0$ にノード数 $N_0$ , リンク(無向辺)数 $M_0$ の単純・連結とは限らない(多重辺や自己ループを含む可能性がある)ネットワークを準備する。ノードには $1$ から $N_0$ の番号が振られているが、リンクは区別されない。このネットワークは、以下の制約を満たすネットワークを重複無く列挙し、その中から一様ランダムに $1$ つ選んだものである。ただし、 $2$ つのネットワークが異なるとはある $i,j$ が存在して、ノード $i$ とノード $j$ を結ぶリンクの本数が異なることとする。

次数が $1$ 以上のノードのうち番号が一番大きいものをノード $K$ として、任意の $i (1\leq i\leq K−1)$ についてノード $i$ とノード $i+1$ を結ぶリンクが存在する。

この後、時刻 $0.5,1.5,2.5,\cdots$ にひとつずつ、 $2$ 本のリンクを持ったノードを追加していく。(これらのノードの番号は順に $N_0+1,N_0+2,\cdots$ とする。)

すなわち、時刻 $t\in Z\geq 0$ におけるネットワークの全ノード数と総リンク数はそれぞれ $N(t)=N_0+t,M(t)=M_0+2t$ である。

時刻 $t+0.5$ に追加されるノードの持つ $2$ 本のリンクがノード $i (1\leq i\leq N(t))$ と繋がる確率をそれぞれ $P_1(i),P_2(i)$ とする。

$P_1(i)=\frac{1}{N(t)}$
時刻 $t$ におけるノード $i$ の次数を $k_i$ として、 $P_2(i)=\frac{k_i}{\sum^{N(t)} _ {j=1} k_j}$

時刻 $t=L$ にネットワークが連結である確率を求めよ。

制約

$M_0 = 1 (200$ 点 $)$
$L = 0 (400$ 点 $)$
$2 \leq N_0 \leq 100$
$1 \leq M_0 \leq 10$
$0 \leq L \leq 10^{4}$
入力はすべて整数

解説

問題の整理

この問題、問題文の読解がとても大変ですよね（この問題は元ネタからこんな調子だったので仕方ないところがあるのですが）。こういうときは、絵を描きながら１文ずつ読んでいくと理解の助けになると思います。グラフ問題は特に絵を書くことの恩恵を受けやすいです。この問題の内容を絵で表すとこんな感じになります。まずは問題文前半で説明されている初期状態です。

f:id:kefeme:20210320140937p:plain

要するに、頂点数 $N_0$ , 辺数 $M_0$ で、連結なカタマリが1つだけあるグラフなのですね。ただし、図では明示できていない条件として連結なカタマリの中では隣り合う数字が辺でつながれていなければなりません（つまり、1-2, 2-3, 3-4, ...という辺がないといけない）。

次に、毎時発生する遷移の様子です。

f:id:kefeme:20210320143639p:plain

毎時、腕を $2$ 本持った新しい頂点が飛んできます。それぞれの腕がどの頂点とつながるかについては、問題文のゴツい数式に気圧されそうですが、落ち着いて読んでみると単に確率を (場合の数) / (総数) という形で表しているに過ぎず、片方の腕は「接続先の次数に比例した確率」、もう片方の腕は「等確率」で各頂点につながるということが読み取れます。

そして、以上のルールをもとに我々が最終的に答えなければいけないものが、初期状態としてありえるものからランダムに $\mathbf{1}$ つピックアップし、 $\mathbf{L}$ 回の遷移を行ったあと、グラフ全体が連結になっている確率です。難解な問題文からシンプル?な問題設定が露わになりました。

全体の方針

初期状態と遷移の仕方が定められているあたりがとてもDPっぽいですね。ある時刻の状態（ここでは状態＝確率分布）から次の時刻の状態が一意に定まるので、時間を添字に持つDPを用いることを考えます。また、制約が $L \leq 10^{4}$ なので、単なる時間の1次元DPにとどまらず添字を増やしていく余裕はあるということを頭の隅に置いておきます。察するに、おそらく連結成分の大きさを表す $K$ も添え字に持ってDPするのでしょうね。

まずは遷移を掘り下げていきます。我々にとって興味があるのは「時刻 $L$ で全体が連結かどうか」なので、「連結」という視点を念頭に置いておきます。

遷移の確率を考える

「接続先の次数に比例した確率でくっつく」というのが明らかに意味深な問題設定なので、その意味するところを考えてみましょう。

もったいぶらずに言ってしまうと、「次数で結合先を選ぶ腕」は次数が $0$ の頂点、つまり孤立した頂点にはつながらないということです。逆に、すでに連結になっているカタマリにしかつながらないとも言えます。

となると、遷移の結果を左右する権利を握っているのは「等確率で全頂点につながるもう片方の腕」です。これが連結成分内の頂点にくっついた場合は何も起きないに等しいのですが、その一方、孤立した頂点にくっついた場合は連結成分に新たに1つ頂点が加わることになります。

それぞれの事象が発生する確率は、連結した頂点数と孤立した頂点数を数えてそれぞれ全頂点数で割ることで $\frac{N-N _ {0}+K}{N}, \frac{N _ {0}-K}{N}$ と求まります。時刻 $t$ における全頂点数を $N$ とおいています。また、 $K$ は上で述べたように、今現在の連結成分の大きさを表しています（問題文の $K$ とは異なるものであることに注意してください）。

f:id:kefeme:20210320165307p:plain

初期状態を計算する

さぁ数え上げです！私はここで一生つまずいていました……

最初に載せた画像を再掲します。

f:id:kefeme:20210320140937p:plain

この条件を満たす辺の張り方の数を連結成分の大きさ $K$ ごとに数え上げたいです。なお、頂点は区別して辺は区別しません。

まず、1-2, 2-3, 3-4, ...の辺は絶対に必要なので先に張ってしまいましょう。これで、 $K-1$ 本の辺を消費します（植木算注意！）。

残った $M _ {0} - K + 1$ 本の辺を多重辺自己ループなんでもありで張っていきます。この際、下手な数え方をすると重複の処理で爆死するので気をつけてください。

上手な数え方としては、張れる辺 (1-1, 1-2, 1-3, ...) を列挙し、それぞれに $M _ {0} - K + 1$ 本の辺を分配するというやり方があります。これは数え上げ典型問題であるところの「箱にボール分ける問題（箱区別してボール区別しないで空箱許すバージョン）」に帰着できます。知ってる人には「重複組合せ」「 $_{n} \mathrm{H} _{r}$ 」で通じるやつです。

f:id:kefeme:20210320182106p:plain

これをこの問題に応用すると、まず張れる辺の種類は $K$ 頂点から重複を許して $2$ 頂点を選んでその間に辺を張るため $_{K+2-1} \mathrm{C} _{2} = \frac{1}{2} K (K + 1)$ 種類となり、その中で $M _ {0} - K + 1$ 本の辺を分配するので、求める場合の数は

$_{K(K-1)/2 + M _ {0}} \mathrm{C} _{M _ {0} - K + 1}$

となります。

DPに落とし込む

DP配列の定義

$dp[t][k] =$ 時刻 $t$ 時点で連結成分内の頂点 $1\sim N _ {0}$ が $k$ であるような確率

初期値

$dp[0][k] = _{k(k-1)/2 + M _ {0}} \mathrm{C} _{M _ {0} - k + 1}$

(確率に直すためにあとで総和で割ります)

漸化式

$dp[t + 1][k + 1] = dp[t][k] \times \frac{N _ {0}-K}{N _ {0}+t} + dp[t][k + 1] \times \frac{t + k + 1}{N _ {0}+t}$

もらうDPです。左の項は孤立していた頂点が選ばれて連結成分に1つ加わる場合、右の項は連結成分が選ばれて変わらない場合の遷移です。

更新順序

素直に添字が小さい順で問題ありません。

実装例

from math import comb

def solve():
    # input
    n, m, l = map(int, input().split())

    # DP配列の生成
    dp = [[0] * (n + 10) for _ in range(l + 10)]

    # 初期値の設定
    dp[0] = [0] * (n + 10)
    for k in range(min(n + 1, m + 2)): # nかm+1で打ち止め
        dp[0][k] = comb(k * (k - 1) // 2 + m, m - k + 1)
    init_sum = sum(dp[0])
    dp[0] = [x / init_sum for x in dp[0]]

    # DPの更新
    for t in range(l):
        for k in range(n):
            dp[t + 1][k + 1] = dp[t][k] * (n - k) / (n + t) + dp[t][k + 1] * (t + k + 1) / (n + t)
    
    # output
    print(f"{dp[l][n]:.15f}") # 桁数を固定しないと指数表記になってWAを食らうことがある


t = int(input())
for _ in range(t): solve()

ひとこと

数え上げ地力が足りない

くふむ

けふぇめ