Project Euler 690

发表于 分类于 阅读次数： 本文字数： 4.1k 阅读时长 ≈ 4 分钟

Project Euler 690

题目

Tom and Jerry

Tom (the cat) and Jerry (the mouse) are playing on a simple graph \(G\).

Every vertex of \(G\) is a mousehole, and every edge of \(G\) is a tunnel connecting two mouseholes.

Originally, Jerry is hiding in one of the mouseholes.

Every morning, Tom can check one (and only one) of the mouseholes. If Jerry happens to be hiding there then Tom catches Jerry and the game is over.

Every evening, if the game continues, Jerry moves to a mousehole which is adjacent (i.e. connected by a tunnel, if there is one available) to his current hiding place. The next morning Tom checks again and the game continues like this.

Let us call a graph \(G\) a Tom graph, if our super-smart Tom, who knows the configuration of the graph but does not know the location of Jerry, can guarantee to catch Jerry in finitely many days. For example consider all graphs on 3 nodes:

For graphs \(1\) and \(2\), Tom will catch Jerry in at most three days. For graph \(3\) Tom can check the middle connection on two consecutive days and hence guarantee to catch Jerry in at most two days. These three graphs are therefore Tom Graphs. However, graph \(4\) is not a Tom Graph because the game could potentially continue forever.

Let \(T(n)\) be the number of different Tom graphs with \(n\) vertices. Two graphs are considered the same if there is a bijection \(f\) between their vertices, such that \((v,w)\) is an edge if and only if \((f(v),f(w))\) is an edge.

We have \(T(3) = 3, T(7) = 37, T(10) = 328\) and \(T(20) = 1416269\).

Find \(T(2019)\) giving your answer modulo \(1\,000\,000\,007\).

解决方案

1. 规则重述与目标

给定简单图 \(G\)。Jerry 起始在某个顶点；每天早晨 Tom 选择且仅选择一个顶点检查，若 Jerry 正在该点则被抓；若未抓到，则当天晚上 Jerry 必须沿一条边移动到相邻顶点。Tom 只知道图结构，不知道 Jerry 的位置，但要设计一套检查序列保证有限天内必抓到。

称 \(G\) 为 Tom graph 当且仅当 Tom 有必胜策略。题目要求按同构计数的 \(T(n)\)：有 \(n\) 个顶点的 Tom graphs 数量（无标号），求 \(T(2019)\bmod 10^9+7\)。

若 \(G\) 含有一个环 \(C\)，Jerry 可以一直限制在该环上移动。环上每个顶点度为 2，在已知 Tom 当天检查点 \(c_t\) 的前提下，Jerry 总能从当前顶点选择一个邻点 \(\neq c_t\) 作为下一晚位置，从而永远不被抓到；因此只要图含环，Tom 就无法保证有限步必胜。论文记 \(m(G)\) 为Tom 保证抓到 Jerry 所需的最少检查次数（若不存在则记为 \(\infty\)），更一般地，若 \(G\) 含环则 \(m(G)=\infty\)并推出若 \(m(G)<\infty\) 则 \(G\) 必为树（连通情形）在论文中已给出严格论证。因此：Tom graph 一定是森林。

论文引入一个关键禁用子树 \(T^\ast\)：三条长度为 3 的路径共享同一端点（三臂各 3 条边的树，形式化定义：存在一个顶点 \(r\)（中心点存在三条路径\(r-u_1-v_1-w_1, r-u_2-v_2-w_2, r-u_3-v_3-w_3\)，它们都是长度为 \(3\) 的路径）。并证明：

  • 若一棵树包含 \(T^\ast\) 作为子图，则 Jerry 可对任何 Tom 检查序列构造躲避序列，可无限生存。
  • 若树是 \(T^\ast\)-free（不含 \(T^\ast\) 子图），则 Tom 存在必胜策略（且还能给出最优抓捕时间）。
  • 更关键的是一个等价结构刻画：一棵树 \(T\) 是 \(T^\ast\)-free 当且仅当存在一条路径 \(P\)，使得树上任意顶点到 \(P\) 的距离 \(\le 2\)。而这就是lobster 树的定义。

因此得到结论：对任意连通简单图 \(G\)，\(G\) 是 Tom graph 当且仅当 \(G\) 是 lobster 树。再加上 1 个点、2 个点的退化情形也都是可保证抓到的（显然或两天内抓到），于是连通 Tom graphs 的计数就是无标号lobster 树的计数。

森林中 Jerry 永远不会跨越连通分量；若某分量不是 Tom graph，则 Jerry 起始在该分量即可永远躲避；反之若每个分量都是 Tom graph，则 Tom 可以把每个分量的一套必胜检查序列按任意顺序串联执行：若 Jerry 在某分量中，他终会在对应那段检查中被捕，其他分量的检查只是空耗但不影响必胜性。

故：对任意简单图 \(G\)，\(G\) 是 Tom graph 当且仅当 \(G\) 的每个连通分量都是 lobster 树（即 \(G\) 是 lobster 森林）。于是问题转化为按同构计数的无标号lobster 森林有多少个。

设\(a_m\)表示大小为 \(m\) 的连通Tom graphs（即lobster 树）同构类数量；\(T(n)\)：大小为 \(n\) 的 Tom graphs（lobster 森林）同构类数量。对固定的 \(m\)，若森林里恰有 \(k\) 个大小为 \(m\) 的连通分量，那么这些分量只是一组 可重元素的多重集合：从 \(a_m\) 种类型里取 \(k\) 个（允许重复、无序），数量是\(\dbinom{a_m+k-1}{k}.\)对应生成函数贡献为\(\displaystyle\sum_{k\ge 0}\binom{a_m+k-1}{k}x^{mk}=\frac{1}{(1-x^m)^{a_m}}.\)将所有 \(m\) 独立相乘，得到总生成函数\(\displaystyle\sum_{n\ge 0}T(n)x^n=\prod_{m\ge 1}\frac{1}{(1-x^m)^{a_m}}.\)这正是 Euler transform 的形式。

进一步为了计算系数，可对上式取对数并做形式幂级数运算。令\(\displaystyle F(x)=\sum_{n\ge 0}T(n)x^n, T(0)=1.\) 则\(\displaystyle \log F(x)=\sum_{m\ge 1}a_m\sum_{k\ge 1}\frac{x^{mk}}{k}=\sum_{n\ge 1}\frac{1}{n}\left(\sum_{d\mid n}da_d\right)x^n.\)

定义\(\displaystyle c_n=\sum_{d\mid n}da_d.\)对 \(\log F\) 求导得到\(\displaystyle \frac{F'(x)}{F(x)}=\sum_{n\ge 1}c_n x^{n-1}.\)比较系数（将 \(\displaystyle F'(x)=\sum_{n\ge 1}nT(n)x^{n-1}\) 与右侧 \(\displaystyle\left(\sum_{k\ge 0}T(k)x^k\right)\left(\sum_{j\ge 1}c_j x^{j-1}\right)\) 相乘）得递推\(\displaystyle nT(n)=\sum_{k=1}^{n}c_kT(n-k),(n\ge 1),\) 也就是：

\[ T(n)=\frac{1}{n}\sum_{k=1}^{n}c_kT(n-k). \]

无标号lobster 树的序列是 OEIS A130131。其条目给出（等价的）可计算生成函数描述：以分拆生成函数\(\displaystyle P(x)=\prod_{k\ge 1}\frac{1}{1-x^k}\)为基础，定义\(M(x)=xP(x),E(x)=x\Big(P(x)-\dfrac{1}{1-x}\Big),S(x)=\dfrac{1}{1-M(x)}.\)则直径 \(>4\)部分为

\[ B(x)=\frac{E(x)^2S(x)+E(x^2)S(x^2)(1+M(x))}{2}, \]

直径 \(\le 4\)部分为 \[ C(x)=xP(x)-\frac{x^3}{(1-x)^2(1+x)}, \]

因此总生成函数 \[ A(x)=\sum_{n\ge 1}a_n x^n=B(x)+C(x). \]

\(P(x)\) 的系数是分拆数 \(p(n)\)，用Euler 五边形数定理的递推可在 \(O(n\sqrt n)\) 内求得 \(p(0,\dots,N)\)。有了 \(P(x)\) 的截断多项式（到 \(x^{N}\)），就能按上面的代数组合式做若干次多项式乘法与幂级数求逆，得到 \(A(x)\) 的前 \(N\) 项，也即 \(a_1,\dots,a_{N}\)。

代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
import numpy as np

N = 2019
MOD = 1000000007


def partitions(n):
    p = [0] * (n + 1)
    p[0] = 1
    g = []
    sgn = []
    k = 1
    while True:
        g1 = k * (3 * k - 1) // 2
        if g1 > n:
            break
        g.append(g1)
        sgn.append(1 if k & 1 else -1)
        g2 = k * (3 * k + 1) // 2
        if g2 <= n:
            g.append(g2)
            sgn.append(1 if k & 1 else -1)
        k += 1
    for i in range(1, n + 1):
        s = 0
        for gg, ss in zip(g, sgn):
            if gg > i:
                break
            s += ss * p[i - gg]
        p[i] = s % MOD
    return p


def solve(n):
    P = partitions(n)

    def arr(lst):
        a = np.array(lst, dtype=object)
        if a.size < n + 1:
            b = np.zeros(n + 1, dtype=object)
            b[:a.size] = a
            return b
        return a[:n + 1].copy()

    def add(a, b):
        return (a + b) % MOD

    def sub(a, b):
        return (a - b) % MOD

    def mul(a, b):
        return np.convolve(a, b)[:n + 1] % MOD

    def inv_series(a):
        a = arr(a)
        inva0 = pow(int(a[0]), MOD - 2, MOD)
        res = np.zeros(n + 1, dtype=object)
        res[0] = inva0
        for i in range(1, n + 1):
            s = int(np.dot(a[1:i + 1], res[i - 1::-1])) % MOD
            res[i] = (-s * inva0) % MOD
        return res

    def sub_x2(a):
        res = np.zeros(n + 1, dtype=object)
        res[::2] = a[:n // 2 + 1]
        return res

    one = np.zeros(n + 1, dtype=object)
    one[0] = 1

    M = np.zeros(n + 1, dtype=object)
    M[1:] = np.array(P[:-1], dtype=object)

    G = inv_series([1, MOD - 1])
    xG = np.zeros(n + 1, dtype=object)
    xG[1:] = G[:-1]

    E = sub(M, xG)
    S = inv_series(sub(one, M))
    E2 = sub_x2(E)
    S2 = sub_x2(S)

    B = add(mul(mul(E, E), S), mul(mul(E2, S2), add(one, M)))
    B = (B * ((MOD + 1) // 2)) % MOD

    invden = inv_series([1, MOD - 1, MOD - 1, 1])
    x3inv = np.zeros(n + 1, dtype=object)
    x3inv[3:] = invden[:-3]
    A = add(B, sub(M, x3inv))

    invs = [0] * (n + 1)
    invs[1] = 1
    for i in range(2, n + 1):
        invs[i] = MOD - (MOD // i) * invs[MOD % i] % MOD

    c = np.zeros(n + 1, dtype=object)
    for d in range(1, n + 1):
        ad = int(A[d]) % MOD
        if ad == 0:
            continue
        t = (d * ad) % MOD
        for m in range(d, n + 1, d):
            c[m] = (c[m] + t) % MOD

    T = np.zeros(n + 1, dtype=object)
    T[0] = 1
    for i in range(1, n + 1):
        s = int(np.dot(c[1:i + 1], T[i - 1::-1])) % MOD
        T[i] = (s * invs[i]) % MOD

    return int(T[n]) % MOD


print(solve(N))