Project Euler 280

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

Project Euler 280

题目

Ant and seeds

A laborious ant walks randomly on a \(5\times5\) grid. The walk starts from the central square. At each step, the ant moves to an adjacent square at random, without leaving the grid; thus there are \(2, 3\) or \(4\) possible moves at each step depending on the ant’s position.

At the start of the walk, a seed is placed on each square of the lower row. When the ant isn’t carrying a seed and reaches a square of the lower row containing a seed, it will start to carry the seed. The ant will drop the seed on the first empty square of the upper row it eventually reaches.

What’s the expected number of steps until all seeds have been dropped in the top row?

Give your answer rounded to \(6\) decimal places.

解决方案

令\(M=2,N=5\)。

令坐标为\((x,y)\)，其中：

  • \(x\in\{0,1,\cdots,N-1\}=\mathbb{Z}_N\)
  • \(y\in\mathbb{Z}_N\)
  • 底行为\(y=0\)，顶行为\(y=N-1\)

接下来用两个\(N\)位掩码（或等价的集合）来描述“尚未完成的列”：

  • \(D\subseteq\mathbb{Z}_N\)：当前底行仍需触发的列集合（即底行仍“有事要做”的列）。
  • \(U\subseteq\mathbb{Z}_N\)：当前顶行仍未填满的列集合（即顶行仍为空的列）。

定义期望：\(E_{x,y}(U,D)\)表示从状态\((x,y,U,D)\)出发，直到全部种子都放到顶行所需的期望步数。

初始状态为：\((x,y)=(M,M)\)，且\(U=D=\mathbb{Z}_N\)。

显然，若直接把所有状态都铺开，状态数至少达到\(N^2\cdot(2^N+N\cdot 2^{N-1})\)，规模较大。为减轻计算负担，当宏状态\((U,D)\)固定时，未知量仅有\(N^2\)个，即\(\{E_{x,y}(U,D)\mid x,y\in\{0,1,2,3,4\}\}.\)因此可将整体系统按宏状态\((U,D)\)进行分块：

  • 对于普通格，其方程只会引用同一个\((U,D)\)下的相邻位置期望，因此在这\(N^2\)个未知量内闭合。
  • 对于事件格，其方程形如\(E_{x,0}(U,D)=E_{x,N-1}(D\setminus{x},U),\)右端属于另一个宏状态\((U',D')=(D\setminus{x},U)\)，从而形成不同块之间的连接。

也就是说，我们对任意一种状态都可以进行一个小规模的线性系统（\(N\times N\)个未知数）的计算，而不必把全部状态一次性纳入系统。

最终，在真实过程里存在两类触发：

  • 不携带种子：到达底行\((x,0)\)且该格有种子时，会拾起；
  • 携带种子：第一次到达顶行\((x,N-1)\)且该格为空时，会放下。

为避免额外引入“携带/不携带”的布尔状态，我们使用如下等价变换统一事件：

当蚂蚁到达当前底行\((x,0)\)且\(x\in D\)时，立刻触发一次事件，并执行：

  1. 清掉该列：\(D' = D\setminus{x}\)；
  2. 上下翻转位置：\((x,y)\mapsto(x,N-1-y)\)。此处\(y=0\)，所以\((x,0)\mapsto(x,N-1)\)；
  3. 交换上下行任务：\((U,D)\mapsto(D',U)\)。

这样，“拾起后需要去顶行放下”的阶段被翻译成“在翻转后的视角里，仍旧是到达底行触发一次事件”。因此我们只需保留一种事件规则：始终只判断\(y=0\)且\(x\in D\)。

综上，可得到统一的状态转移方程：

\[ E_{x,y}(U,D)= \left\{ \begin{aligned} &0, & & \text{if}\quad U=\varnothing,\\ &E_{x,N-1-y}(D\setminus{x},U), & & \text{else if}\quad y=0\land x\in D,\\ &1+\frac{1}{\deg(x,y)}\sum_{(u,v)\in N(x,y)}E_{u,v}(U,D), & & \text{else}. \end{aligned} \right. \]

其中第一行是终止条件（吸收态）：\(U=\varnothing\Rightarrow E_{x,y}(U,D)=0\)；第二行表示若当前位置满足\(y=0\)且\(x\in D\)，则立即触发事件且不消耗步数；第三行表示否则走一步到相邻格。这里\(\deg(x,y)\)为可走邻居数，\(N(x,y)\)为相邻位置集合。在\(N\times N\)棋盘上，\(\deg(x,y)\)仅有三类取值：

  • 角点：\(\deg(x,y)=2\)；
  • 边但非角：\(\deg(x,y)=3\)；
  • 内部：\(\deg(x,y)=4\)。

最终所求答案即为\[E_{M,M}(\mathbb{Z}_N,\mathbb{Z}_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
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
from functools import cache
import numpy as np

# 网格边长 N（本题 N=5）
N = 5

# 将二维坐标 (x,y) 映射到一维编号 p ∈ [0, N^2)
# 约定：x 为列（0..N-1），y 为行（0..N-1），y=0 为底行，y=N-1 为顶行
idx = [[0 for _ in range(N)] for _ in range(N)]
for i in range(N):
    for j in range(N):
        idx[i][j] = i + j * N


@cache
def dfs(sx, sy, bits_u, bits_d):
    """
    返回期望步数 E(sx, sy; bits_u, bits_d)。

    状态变量含义：
      - (sx, sy): 蚂蚁当前位置。
      - bits_d: “当前底行 y=0 仍需触发的列集合”，用 N-bit 掩码表示。
               第 x 位为 1 表示：当蚂蚁到达 (x,0) 时会立刻发生一次“事件”（拾起或放下）。
      - bits_u: “当前顶行 y=N-1 仍需触发的列集合”，用 N-bit 掩码表示。
               它不会直接在转移方程里被检查，而是在事件发生后与 bits_d 交换角色。

    为什么只检查 y==0 的事件就够了？
      通过“上下翻转 + 交换 bits_u/bits_d”的等价变换，把
        “携带种子去顶行找空位放下”
      变成在翻转后的世界里
        “不携带种子去新底行触发一次事件”
      从而统一用“到达 y==0 且 bits_d 有该列”这一种事件处理。

    本函数的核心是：对固定 (bits_u, bits_d) 建立 N^2 个线性方程，未知量为各格子的期望 E[x,y]，
    解出整张网格的 E，再返回起点 (sx,sy) 的那一项。
    """

    # 线性方程组：M * E = B
    #
    # 未知向量 E 是长度 N^2 的列向量：
    #   E[p] = E[x,y]，其中 p = idx[x][y]
    #
    # M 是系数矩阵 (N^2 × N^2)，B 是常数项列向量 (N^2 × 1)。
    M = np.zeros((N**2, N**2), dtype=float)
    B = np.zeros((N**2, 1), dtype=float)

    for x in range(N):
        for y in range(N):
            p = idx[x][y]

            # 我们总是把方程写成 “E[p] - sum(coeff * E[neighbor]) = constant”
            # 先写下 E[p] 的系数为 1
            M[p, p] = 1.0

            # ----------------------------
            # (1) 事件格的方程：y==0 且 bits_d 的第 x 位为 1
            # ----------------------------
            #
            # 事件含义（统一表述）：
            #   当蚂蚁在“当前底行”到达 (x,0) 且该列仍需触发，
            #   会立刻发生一次“状态切换”，且该切换不消耗步数（没有 +1）。
            #
            # 在原问题里，这个切换可能是：
            #   - 不携带时：拾起底行的种子
            #   - 携带时：在顶行空格放下种子
            #
            # 由于我们对状态做了翻转/交换的等价变换，两类切换都可归一到这里处理。
            if y == 0 and (bits_d & (1 << x)):
                # 触发一次事件后，这一列的“底行目标”被完成，因此清除此位
                new_bits_d = bits_d - (1 << x)

                # bits_u == 0 表示顶行目标已经全完成。
                # 在真实过程里，这时通常不会再触发“拾起/放下”逻辑；
                # 这里沿用传统写法：让该方程变成 E[p]=0（吸收态），通过保持 B[p]=0 来达成。
                if bits_u == 0:
                    continue

                # 事件发生后的等价变换：
                #   - 上下翻转坐标：(x, y) -> (x, N-1-y)
                #     这里 y=0，所以翻转后位置是 (x, N-1)（顶行）
                #   - 交换 bits_u / bits_d 的角色：
                #       新的“底行目标集合” = new_bits_d
                #       新的“顶行目标集合” = bits_u
                #
                # 由于事件不消耗步数，所以直接有：
                #   E(x,0; bits_u, bits_d) = E(x,N-1; new_bits_d, bits_u)
                #
                # 我们把右侧当成一个已知常数（递归调用返回的期望），放入 B[p]：
                B[p, 0] = dfs(x, N - 1 - y, new_bits_d, bits_u)
                continue

            # ----------------------------
            # (2) 普通格的方程：一步随机游走
            # ----------------------------
            #
            # 若当前位置不是事件格，则蚂蚁会走一步到相邻格子（不出界），且等概率选择。
            #
            # 设当前位置可走邻居数为 t，邻居集合为 N(p)，则期望满足：
            #
            #   E[p] = 1 + (1/t) * sum_{q in N(p)} E[q]
            #
            # 移项得到标准线性形式：
            #
            #   E[p] - (1/t) * sum_{q in N(p)} E[q] = 1
            #
            # 这就是我们在 M 与 B 中填写的内容：
            #   - B[p] = 1
            #   - 对每个邻居 q，M[p,q] = -(1/t)
            t = 0
            if x < N - 1:
                t += 1
            if x > 0:
                t += 1
            if y < N - 1:
                t += 1
            if y > 0:
                t += 1

            f = 1.0 / t
            B[p, 0] = 1.0

            if x < N - 1:
                M[p, idx[x + 1][y]] = -f
            if x > 0:
                M[p, idx[x - 1][y]] = -f
            if y < N - 1:
                M[p, idx[x][y + 1]] = -f
            if y > 0:
                M[p, idx[x][y - 1]] = -f

    # 解线性方程组得到 E（长度 N^2 的解向量）
    # np.linalg.solve 会用数值稳定的方式求解 M*E=B（比显式求逆更可靠）
    E = np.linalg.solve(M, B)

    # 返回起点 (sx,sy) 的期望步数
    return E[idx[sx][sy], 0]


# 初始状态：
#   - 蚂蚁从中心格出发 (N//2, N//2)
#   - 底行 N 个种子都存在 => bits_d = (1<<N)-1（二进制全 1）
#   - 顶行 N 个空位都未填 => bits_u = (1<<N)-1（二进制全 1）
ans = dfs(N // 2, N // 2, (1 << N) - 1, (1 << N) - 1)
print(f"{ans:.6f}")