Project Euler 762

Project Euler 762

题目

Amoebas in a 2D grid

Consider a two dimensional grid of squares. The grid has 4 rows but infinitely many columns. An amoeba in square \((x, y)\) can divide itself into two amoebas to occupy the squares \((x+1,y)\) and \((x+1,(y+1) \bmod 4)\), provided these squares are empty.

The following diagrams show two cases of an amoeba placed in square A of each grid. When it divides, it is replaced with two amoebas, one at each of the squares marked with B:

Originally there is only one amoeba in the square \((0, 0)\). After \(N\) divisions there will be \(N+1\) amoebas arranged in the grid. An 分布 may be reached in several different ways but it is only counted once. Let \(C(N)\) be the number of different possible 分布s after \(N\) divisions.

For example, \(C(2) = 2\), \(C(10) = 1301\), \(C(20)=5895236\) and the last nine digits of \(C(100)\) are \(125923036\).

Find \(C(100\,000)\), enter the last nine digits as your answer.

解决方案

如果直接把每一步选哪个变形虫去分裂当作状态来搜索,会出现严重的重复计数:因为同一个最终占据图形,往往可以由不同的分裂顺序得到。

因此,我们不按操作过程计数,而是直接按最终图形计数。为此可以采用一种按列从左到右的规范化扫描方法:扫描到某一列时,只保留这条扫描边界(前沿)对右侧尚未处理部分的影响,而不再关心左侧部分是按什么顺序生成出来的。

在这种表示下,需要记录的不是格子里有几个变形虫,而是当前前沿某个格子从左侧接收到多少路贡献(也可以理解为有多少条分裂链在这里汇合)。由于分裂规则决定了一个格子只能由左边相邻两格产生,而本题又只有 \(4\) 行并按模 \(4\) 循环,因此一个前沿格子的左侧贡献数最多只能是 \(2\)

这就意味着,前沿状态中的关键参数只需要区分三种情况:\(0,1,2\)。于是状态空间被大幅压缩,最终可以建立一个规模很小、可直接计算的递推模型。

定义 \(G(k,m)\),其中\(k\)表示后续还需放置(或等价地还需解释)的最终占据格子数;\(m\)表示当前前沿的重叠层数(只可能是 \(0,1,2\))。

则原题答案满足:\(C(N) = G(N+1, 0)\)因为经过 \(N\) 次分裂后最终有 \(N+1\) 个变形虫(即 \(N+1\) 个最终占据格子)。边界条件统一约定为:\(k < 0\)\(G(k,m)=0\)\(m \ge 3\)\(G(k,m)=0\)

将前沿局部结构按是否重叠分类,得到三条递推。为简写,记:\(A_k = G(k,0),B_k = G(k,1),D_k = G(k,2).\)

当前前沿没有待处理重叠(即 \(m=0\))时,考虑最左端尚未解释的局部结构。按规范化扫描后的局部拼接情况分类,可以得到三部分贡献。

第一类是前沿仍然保持无重叠并继续向右推进。这个局部有两种对称的选择,因此给出 \(2A_{k-1}\)。其中下标变成 \(k-1\),表示该局部会确定(消耗)\(1\) 个最终占据格子;第二类是局部结构会把前沿状态从 \(m=0\) 转入 \(m=1\),对应贡献 \(B_k\)。这一步只改变前沿类型,不额外减少剩余格子数,所以仍然是下标 \(k\);第三类是一个只在起始阶段出现的最小特例:当 \(k=2\) 时,存在一个不能被前两类覆盖的原始局部构型,因此需要额外加上 \(1\)。因此有

\[ A_k = 2A_{k-1} + B_k + \mathbf{1}\{k=2\}. \]

当前前沿处于单层重叠状态(\(m=1\))时,同样看最左端局部结构,并按它对后续前沿的影响分类。

有一种情况会把当前重叠消解掉,使状态回到无重叠,对应 \(A_{k-3}\);这里下标减 \(3\),表示该局部净确定了 \(3\) 个最终占据格子;有两种对称的情况会保持单层重叠并继续向右推进,因此贡献是 \(2B_{k-2}\);还有一种情况会使重叠加深到两层,对应 \(D_{k-1}\);最后还会出现一类由 \(4\) 行环状结构(\(y\)\(\bmod 4\) 循环)带来的回绕拼接,它仍然落在单层重叠状态,对应 \(B_{k-4}\)。这一项正是本题与非环状模型相比的额外项。综上得到

\[ B_k = A_{k-3} + 2B_{k-2} + D_{k-1} + B_{k-4}. \]

当前前沿处于双层重叠状态(\(m=2\))时,由于本题中一个前沿格子的左侧贡献数最多为 \(2\),不会再出现新的有效状态 \(m=3\),因此这里只剩两类转移。

一种情况会把双层重叠部分消解为单层重叠,对应 \(B_{k-4}\);另一种情况保持双层重叠并向右推进,这种局部有两种对称选择,因此贡献 \(2D_{k-3}\)。于是递推在这里截断为:

\[ D_k = B_{k-4} + 2D_{k-3}. \]

接下来把 \(A_k,B_k,D_k\) 消元,得到只含 \(A_k\) 的递推。我们使用生成函数完成。定义普通生成函数:\(\displaystyle A(x)=\sum_{k\ge0} A_k x^k,B(x)=\sum_{k\ge0} B_k x^k,D(x)=\sum_{k\ge0} D_k x^k\)利用负下标项视为 \(0\) 的约定,将三条递推逐项求和,可得:

  • \(A_k = 2A_{k-1} + B_k + \mathbf{1}\{k=2\}\)可得\((1-2x)A(x) - B(x) = x^2;\)
  • \(B_k = A_{k-3} + 2B_{k-2} + D_{k-1} + B_{k-4}\)可得\(-x^3A(x) + (1-2x^2-x^4)B(x) - xD(x) = 0;\)
  • \(D_k = B_{k-4} + 2D_{k-3}\)可得\(-x^4B(x) + (1-2x^3)D(x) = 0.\)

\(B(x), D(x)\) 消去后,可得到 \(A(x)\) 是一个有理函数,其分母为:\(1 - 2x - 2x^2 + x^3 + 3x^4 + 5x^5 - 4x^6 + 2x^7 - 4x^8\)

因此 \(A_k\) 满足 8 阶线性递推。再利用 \(C(n)=A_{n+1}\),得到同样的递推(仅下标整体平移):

\[ C(n) = 2C(n-1) + 2C(n-2) - C(n-3) - 3C(n-4) - 5C(n-5) + 4C(n-6) - 2C(n-7) + 4C(n-8) \]

从上面的状态递推直接小范围计算可得:\(C(0,\dots,9)=[1,1,2,4,9,20,46,105,243].\)定义长度为 \(8\) 的列向量状态:

\[ S(n) = [C(n), C(n-1), C(n-2), C(n-3), C(n-4), C(n-5), C(n-6), C(n-7)]^T \]

则对 \(n \ge 9\),有 \(S(n)=M S(n-1)\),其中转移矩阵 \(M\) 如下:

\[ M= \begin{bmatrix} 2 & 2 & -1 & -3 & -5 & 4 & -2 & 4 \\ 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 \end{bmatrix} \]

于是\(S(8) = [C(8),C(7),...,C(1)]^T = [243,105,46,20,9,4,2,1]^T,S(N) = M^{(N-8)} S(8)\)

答案就是 \(S(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
N = 100000
MOD = 10 ** 9

def mat_mul(A, B, mod):
return [[sum(A[i][k] * B[k][j] for k in range(len(A[0]))) % mod for j in range(len(B[0]))] for i in range(len(A))]


def mat_pow(a, e):
n = len(a)
r = [[1 if i == j else 0 for j in range(n)] for i in range(n)]
while e:
if e & 1:
r = mat_mul(r, a, MOD)
a = mat_mul(a, a, MOD)
e >>= 1
return r

def C(n):
if n == 0:
return 1
base = [1, 2, 4, 9, 20, 46, 105, 243]
if n <= 8:
return base[n - 1] % MOD
coef = [2, 2, -1, -3, -5, 4, -2, 4]
a = [[0] * 8 for _ in range(8)]
a[0] = [x % MOD for x in coef]
for i in range(1, 8):
a[i][i - 1] = 1
s8 = [[base[7]], [base[6]], [base[5]], [base[4]], [base[3]], [base[2]], [base[1]], [base[0]]]
p = mat_pow(a, n - 8)
return mat_mul(p, s8, MOD)[0][0] % MOD

print(C(N))