グラフが連結である確率（3）

グラフ プログラミング

m点で辺の数がm-1で連結のとき、グラフは木になる。ループがなくて連結のグラフを木と呼ぶ。
m点で辺の数がm-1のグラフの数を数える。辺の種類はm(m-1)/2ある。ここからm-1辺を選ぶから、



m(m-1)/2Cm-1

そのうち、連結なグラフの数をT(m)で表す。前々回、T(4) = 16、T(5) = 125であることを確認した。一般にはどうなるだろう。


m点のうちの1点を取り出す。これは他の点と連結している。1点かもしれないし、2点かもしれない。1点だとすると、この点を取り除くとこのグラフは連結である。だから、1点で繋がっている場合は、T(m-1)あることになる。2点なら連結のグラフが2つあることになる。2つのグラフは点の数が[ m - 2, 1 ], [ m - 3, 2 ], ...となる。要するに、m-1の分割を列挙する必要がある。それをPythonで再帰を使って書く。



def divide(n, l = 0):

    if n == 1:

        return [ [ 1 ] ]

    elif l == 0:

        l = n

    elif n < l:

        l = n

    

    b = [ ]

    if l < n:

        mx = l + 1

    else:

        mx = l

    for i in range(1, mx):

        for d in divide(n - i, min(i, l)):

            b.append([ i ] + d)

    if l == n:

        b.append([ n ])

    return b

divide(4)とすると、



[[1, 1, 1, 1], [2, 1, 1], [2, 2], [3, 1], [4]]

と返ってくる。例えば、[2, 2]を考えると、4つの中から2つを選ぶから、4C2、残りの2つから2つを選ぶから、2C2。ただし、2つのグループが入れ替わっても同じだから、4C22C2。また考えている点とそれぞれのグループと結ぶ辺は、それぞれ2通りずつある。また、それぞれのグループの中で連結だからT(2)通りある。結局、



4C22C2/2!2*2*T(2)*T(2) = 12

12通りある。このように計算してすべての分割を足し合わせる。
コードは下に示すが、このうちconvertは[2, 2]→[0, 0, 2]と変換する。2が2つありほかは0だからこのようになる。重複を考えて、[0!, 0!, 2!]で割るためのものである。
2時間くらい回したが、T(69)までしか出なかった。



T(2) = 1 p(2) = 1/1 = 1.0

T(3) = 3 p(3) = 3/3 = 1.0

T(4) = 16 p(4) = 16/20 = 0.8

T(5) = 125 p(5) = 125/210 = 0.595238095238

T(6) = 1296 p(6) = 1296/3003 = 0.431568431568

T(7) = 16807 p(7) = 16807/54264 = 0.309726522188

T(8) = 262144 p(8) = 262144/1184040 = 0.221397925746

T(9) = 4782969 p(9) = 4782969/30260340 = 0.158060649682

T(10) = 100000000 p(10) = 100000000/86163135 = 0.112846039347

p(11) = 0.0806173186438

p(12) = 0.0576467330605

p(13) = 0.0412645757851

p(14) = 0.0295699238707

p(15) = 0.0212122549758

 ...

p(67) = 1.29847366758e-009

p(68) = 9.48684441114e-010

p(69) = 6.93190168084e-010




def divide(n, l = 0):

    if n == 1:

        return [ [ 1 ] ]

    elif l == 0:

        l = n

    elif n < l:

        l = n

    

    b = [ ]

    if l < n:

        mx = l + 1

    else:

        mx = l

    for i in range(1, mx):

        for d in divide(n - i, min(i, l)):

            b.append([ i ] + d)

    if l == n:

        b.append([ n ])

    return b
def calcTree(T, m):

    N = 0

    d = divide(m - 1)

    for e in d:

        N += calcN(e)

    T.append(N)
def calcN(d):

    a = convert(d)

    n = fac(sum(d))

    

    for e in d:

        n *= e * T[e]

        n /= fac(e)

    for e in a:

        n /= fac(e)

    return n
def sum(a):

    n = 0

    for e in a:

        n += e

    return n
def convert(d):

    a = [ ]

    for e in d:

        if e >= len(a):

            for i in range(len(a), e):

                a.append(0)

            a.append(1)

        else:

            a[e] += 1

    return a
def fac(n):

    p = 1

    for i in range(1, n + 1):

        p *= i

    return p
def allN(m):

    n = m * (m - 1) / 2

    return fac(n) / fac(m - 1) / fac(n - m + 1)
n = 1000

T = [ 1, 1 ]

for m in range(2, n + 1):

    calcTree(T, m)

    a = allN(m)

    print m, T[m], a, float(T[m]) / a