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

全ての場合について試すのは無理があるので、前々回と同じようにまとめて考える。

その前に、Pythonで補助的なクラスArrayGenを作る。
例えば、階乗の計算はこのように行っていたが、

def fac(n):

    p = 1

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

        p *= i

    return p

これだと、n!ならn-1回の乗算を行わなければならない。しかし、(n-1)!の計算を記憶していれば1階の乗算で済む。このように計算結果を記憶しておき、新たな計算をしなければならないときは漸化式に従って計算する。

class ArrayGen:

    def __init__(self, f, initial):

        self.f = f

        self.a = [ ]

        if type(initial) == list:

            self.a.extend(initial)

        else:

            self.a.append(initial)

    

    def val(self, n):

        if n < len(self.a):

            return self.a[n]

        else:

            v = self.f(self, n)

            self.a.append(v)

            return v
fac = ArrayGen(lambda o, n: n * o.val(n - 1), 1)

for i in range(101):

    print fac.val(i)

ArrayGenの引数は、漸化式になる関数と、初期値である。初期値は0のときの値、あるいは初期値が複数ある場合は配列にすればよい。例えば、フィボナッチ数列なら次のようである。

class ArrayGen:

    ...
fib = ArrayGen(lambda o, n: o.val(n - 1) + o.val(n - 2), [ 0, 1 ])

for n in range(10):

    print fib.val(n)

さて、今度は連結でないグラフを考え、全てのグラフの数から連結でないグラフの数を引いて、連結であるグラフの数を求める。連結でないグラフは、複数の連結グラフに分かれる。n点の分け方を列挙する。

def divide(o, n):

    a = [ [ n ] ]

    for i in range(n - 1, 0, -1):

        for d in o.val(n - i):

            if d[0] <= i:

                a.append([ i ] + d)

    return a
div = ArrayGen(divide, 0)

for n in range(1, 5):

    print div.val(n)

さらに、この分かれたグラフに辺を分配する。ただし、前回見たように、連結でありうる辺の数は上限・下限が決まっている。

def distribute(nmin, nmax, n):

    nlen = len(nmin)

    if nlen == 1:

        if nmin[0] <= n and n <= nmax[0]:

            yield [ n ]

    else:

        res = n

        for a in nmin:

            res -= a

        for i in range(nmin[0], nmax[0] + 1):

            for d in distribute(nmin[1:], nmax[1:], n - i):

                yield [ i ] + d
def dist_edge(a, n):

    nmin = [ i - 1 for i in a ]

    nmax = [ i * (i - 1) / 2 for i in a ]

    for distr in distribute(nmin, nmax, n):

        yield distr

まず、点の数の配列と辺の数をdist_edgeに渡している。そうすると辺の数の上限・下限が決まるので、これをdistributeに渡す。この中で、再帰的に分配をしている。
分配が決まれば、個々の連結のグラフの場合の数は求まっているとすれば、このパターンの全体のグラフの数が求まる。
これを実装したところ、n=22まで正確に求められた。

同じ点と辺の数の比でも点の数が増すにつれて連結である確率が低くなることが分かる。

class ArrayGen:

    def __init__(self, f, initial):

        self.f = f

        self.a = [ ]

        if type(initial) == list:

            self.a.extend(initial)

        else:

            self.a.append(initial)

    

    def val(self, n):

        if n < len(self.a):

            return self.a[n]

        else:

            return self.f(self, n)
def divide(o, n):

    a = [ [ n ] ]

    for i in range(n - 1, 0, -1):

        for d in o.val(n - i):

            if d[0] <= i:

                a.append([ i ] + d)

    o.a.append(a)

    return o.a[n]
def distribute(nmin, nmax, n):

    nlen = len(nmin)

    if nlen == 1:

        if nmin[0] <= n and n <= nmax[0]:

            yield [ n ]

    else:

        res = n

        for a in nmin:

            res -= a

        for i in range(nmin[0], nmax[0] + 1):

            for d in distribute(nmin[1:], nmax[1:], n - i):

                yield [ i ] + d
def dist_edge(a, n):

    nmin = [ i - 1 for i in a ]

    nmax = [ i * (i - 1) / 2 for i in a ]

    for distr in distribute(nmin, nmax, n):

        yield distr
div = ArrayGen(divide, 0)

fac = ArrayGen(lambda o, n: n * o.val(n - 1), 1)
def calcGraph(m, n):

    N = 0

    for d in div.val(m):        # distribution of vertices

        if len(d) == 1:

            continue

        N += calcNumVertexDist(d, n)

    return combination(m * (m - 1) / 2, n) - N
def calcNumVertexDist(d, n):

    N = 0

    for de in dist_edge(d, n):  # distribution of edges

        N += calcNumEdgeDist(d, de)

    

    a = convert(d)

    N *= fac.val(sum(d))

    for e in d:

        N /= fac.val(e)

    for e in a:

        N /= fac.val(e)

    return N
def calcNumEdgeDist(d, e):

    n = 1

    for i in range(len(d)):

        n *= G[d[i]][e[i]]

    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 combination(m, n):

    N = 1

    if n * 2 > m:

        n = m - n

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

        N *= m - i + 1

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

        N /= i

    

    return N
G = [ [ ], [ 1 ] ]

for m in range(2, 30):

    G.append([ ])

    for n in range(m * (m - 1) / 2 + 1):

        allN = combination(m * (m - 1) / 2, n)

        if n < m - 1:

            N = 0

        elif m * (m - 3) / 2 + 1 < n:

            N = allN

        else:

            N = calcGraph(m, n)

        G[m].append(N)

        print "p( %d , %d ) = %d / %d = %.8f" % (m, n, N, allN, float(N) / allN)