OEIS A001349 を計算する Python プログラム

count-unlabeled-connected-graph.py

from fractions import Fraction
import math

# 定数項のみ val 、それ以外は 0
def constant_fps(n, val) :
    res = [Fraction(0) for _ in range(n)]
    res[0] = val
    return res

# 全部ゼロ
def zero_fps(n) :
    return constant_fps(n, Fraction(0))

# 多項式の掛け算
def mult_poly(f, g) :
    n = len(f)
    m = len(g)
    res = zero_fps(n+m-1)
    for i in range(n) :
        for j in range(m) :
            res[i+j] += f[i] * g[j]
    return res

# 加算
# f,g の長さは一致していること
def add_fps(f, g) :
    n = len(f)
    return [f[i] + g[i] for i in range(n)]

# 定数倍
def times_fps(f, t) :
    n = len(f)
    return [f[i] * t for i in range(n)]

# 微分
def derivative_fps(f) :
    n = len(f)
    return [f[i+1] * Fraction(i+1) for i in range(n-1)]

# derivative_fps の逆変換、ただし res[0] = 0
def integral_fps(f) :
    n = len(f)
    return [Fraction(0)] + [f[i-1] * Fraction(1,i) for i in range(1,n+1)]

# exp(f) = 1 / 0! + f / 1! + f^2 / 2! + f^3 / 3! + ...
def exp_fps(f) :
    n = len(f)
    res = constant_fps(n, Fraction(1))
    for j in reversed(range(1,n+1)) :
        res = list(mult_poly(res, f)[0:n])
        res = times_fps(res, Fraction(1,j))
        res = add_fps(res, constant_fps(n,Fraction(1)))
    return res

# fg = 1 となる g を返す
# f[0] != 0 でなければならない
def inv_fps(f) :
    n = len(f)
    f0 = f[0]
    g = times_fps(f, Fraction(1) / f0)
    res = constant_fps(n, Fraction(1))
    for i in range(1, n) :
        for j in range(i) :
            res[i] -= res[j] * g[i-j]
    res = times_fps(res, Fraction(1) / f0)
    return res

# exp_fps の逆変換
#  log f = ∫ f'/f
def log_fps(f) :
    n = len(f)
    return integral_fps(mult_poly(derivative_fps(f), inv_fps(f))[0:n-1])

# 重さ k のアイテムが f[k] 種類ある
#     各種類のアイテムは無数にある
# res[n] : 重さの総和が n となるようなアイテムの組み合わせは何通りあるか？
def euler_transform_fps(f) :
    n = len(f)
    log_ans = zero_fps(n)
    for k in range(1,n) :
        t = 1
        while t*k < n :
            log_ans[t*k] += f[k] * Fraction(1,t)
            t += 1
    return exp_fps(log_ans)

# euler_transform_fps の逆変換
def inverse_euler_transform_fps(f) :
    n = len(f)
    g = log_fps(f)
    ans = zero_fps(n)
    for k in range(1,n) :
        ans[k] = g[k]
        t = 1
        while t*k < n :
            g[t*k] -= ans[k] * Fraction(1,t)
            t += 1
    return ans

# res[n] : 総和が n となる正整数の組み合わせを全列挙
def list_partitions(n) :
    memo = [[] for _ in range(n+1)]
    memo[0].append([])
    for m in range(1, n+1) :
        for k in range(m) :
            for pre in memo[k] :
                if len(pre) <= m-k :
                    tmp = [1] * (m-k)
                    for i in range(len(pre)) :
                        tmp[i] += pre[i]
                    memo[m].append(tmp)
    return memo

# OEIS A000088
# res[n] : n 頂点のラベル無し無向グラフの個数
def count_unlabeled_graphs_fps(n) :
    p = list_partitions(n-1)
    res = zero_fps(n)
    for k in range(n) :
        # k! 通りある置換すべてについてその置換で変化しない
        #   辺の張り方を数えて足し上げ、 k! で割る

        # 置換のサイクルのサイズのパターンを全パターン
        for partition in p[k] :
            # 辺の張り方
            g2 = 0
            for a in range(len(partition)) :
                for b in range(a) :
                    g2 += math.gcd(partition[a], partition[b])
            for cy in partition :
                if cy >= 2 :
                    g2 += cy // 2
            numer = 2 ** g2

            # サイクルの配置の仕方の個数の 1/k! 倍、の逆数
            denom = 1
            counter = [0] * (k+1)
            for cy in partition :
                denom *= cy
                counter[cy] += 1
            for z in range(k+1) :
                denom *= math.factorial(counter[z])

            res[k] += Fraction(numer, denom)
    return res

# OEIS A001349
# res[n] : n 頂点のラベル無し "連結" 無向グラフの個数
def count_unlabeled_connected_graphs_fps(n) :
    res = inverse_euler_transform_fps(count_unlabeled_graphs_fps(n))
    res[0] = Fraction(1)
    return res

# 40 ならすぐに答えが出るはず

# f = count_unlabeled_graphs_fps(15)
# print(list(a.numerator for a in f))
# 
# g = count_unlabeled_connected_graphs_fps(15)
# print(list(a.numerator for a in g))

f = count_unlabeled_connected_graphs_fps(45)
for d in range(len(f)) :
    print(d, ":", f[d])