Project Euler 380
Project Euler 380
题目
Amazing Mazes!
An \(m\times n\) maze is an \(m\times n\) rectangular grid with walls placed between grid cells such that there is exactly one path from the top-left square to any other square.
The following are examples of a \(9\times12\) maze and a \(15\times20\) maze:
Let \(C(m,n)\) be the number of distinct \(m\times n\) mazes. Mazes which can be formed by rotation and reflection from another maze are considered distinct.
It can be verified that \(C(1,1) = 1, C(2,2) = 4, C(3,4) = 2415,\) and \(C(9,12) = 2.5720e46\) (in scientific notation rounded to \(5\) significant digits).
Find \(C(100,500)\) and write your answer in scientific notation rounded to \(5\) significant digits.
When giving your answer, use a lowercase \(e\) to separate mantissa and exponent.
E.g. if the answer is \(1234567891011\) then the answer format would be \(1.2346e12\).
解决方案
把每个格子看成一个顶点，相邻格子之间连一条边（上下左右）。这得到一个图 \(G_{m,n}\)，它正是 \(m\times n\) 的网格图（grid graph），顶点数为 \(|V|=mn\)。
迷宫“打通通道”的选择，等价于在 \(G_{m,n}\) 的边集中选出一部分边作为“可走通道”。题目要求“从左上到任意点恰有一条路径”，这意味着：
- 图必须连通（否则到不了某些格子）；
- 图必须无环（若存在环，则会出现两条不同路径）。
因此可走通道构成一棵包含全部 \(mn\) 个顶点的树，即 \(G_{m,n}\) 的一棵生成树（spanning tree）。
反之，任意一棵生成树都给出了一个迷宫：把树上的边打通，其余边放墙，则从左上到任意点的路径唯一。
所以\(C(m,n)=\tau(G_{m,n}),\)其中 \(\tau(G)\) 表示图 \(G\) 的生成树数量。
设 \(G\) 的拉普拉斯矩阵为 \(L=D-A\)（\(D\) 为度数对角矩阵，\(A\) 为邻接矩阵）。矩阵树定理给出：
\[ \tau(G)=\det(L^{\ast}), \]
其中 \(L^*\) 是删去 \(L\) 的任意一行一列得到的余子式（对连通图任取都相同）。
另一方面，因为 \(L\) 是实对称半正定矩阵，特征值可写为\(0=\lambda_0<\lambda_1\le\cdots\le \lambda_{|V|-1}.\)对连通图，有经典等价形式 \[ \tau(G)=\frac{1}{|V|}\prod_{k=1}^{|V|-1}\lambda_k. \]
因此对网格图 \(G_{m,n}\)，只要能写出拉普拉斯的全部非零特征值，就能得到闭式乘积。
令 \(P_m\) 为 \(m\) 个点的路径图，就是一条长度为 \(m\) 的“链”，顶点为 \(1,2,\dots,m\)，边为\((1,2),(2,3),\dots,(m−1,m)\)。令\(L_{G}\)表示\(G\)的拉普拉斯矩阵。
网格图满足\(G_{m,n}=P_m\square P_n\) （其中\(\square\)是图论中的笛卡尔积）。
拉普拉斯矩阵在图的笛卡尔积下满足Kronecker 和结构：
\[ L_{G_{m,n}} = L_{P_m}\otimes I_n + I_m\otimes L_{P_n}. \]
其中，运算符\(\otimes\)表示Kronecker 积。由于 \(L_{P_m}\otimes I_n\) 与 \(I_m\otimes L_{P_n}\) 彼此可交换，且二者都是实对称矩阵（因此都可正交对角化），它们可以同时对角化。于是 \(L_{G_{m,n}}\) 的特征值恰好由两边特征值两两相加得到：
若 \(L_{P_m}\) 的特征值为 \(\{\alpha_i\}_{i=0}^{m-1}\)，\(L_{P_n}\) 的特征值为 \(\{\beta_j\}_{j=0}^{n-1}\)，则\(\lambda_{i,j}=\alpha_i+\beta_j, 0\le i\le m-1,\ 0\le j\le n-1\)构成 \(L_{G_{m,n}}\) 的全部特征值（共 \(mn\) 个，按重数计）。因此，问题进一步归结为：求路径图 \(P_m\) 的拉普拉斯矩阵的特征值。
\(P_m\) 的拉普拉斯矩阵\(L_{P_m}\)是一个三对角矩阵，满足：
- 对角元：\(L_{11}=L_{mm}=1\)，且 \(L_{kk}=2\)（\(2\le k\le m-1\)）；
- 相邻元：\(L_{k,k+1}=L_{k+1,k}=-1\)；
- 其余为 \(0\)。
设 \(x=(x_1,\dots,x_m)^T\) 为特征向量，对应特征值 \(\alpha\)，即\(Lx=\alpha x.\)逐分量写出方程：
- 端点 \(k=1:x_1-x_2=\alpha x_1\Longleftrightarrow x_2=(1-\alpha)x_1.\)
- 内部点 \(2\le k\le m-1:-x_{k-1}+2x_k-x_{k+1}=\alpha x_k\Longleftrightarrow x_{k+1}-(2-\alpha)x_k+x_{k-1}=0.\)
- 端点 \(k=m:x_m-x_{m-1}=\alpha x_m\Longleftrightarrow x_{m-1}=(1-\alpha)x_m.\)
为了统一记号，引入两个“虚拟”分量 \(x_0,x_{m+1}\)，规定\(x_0=x_1, x_{m+1}=x_m.\)这样做的原因是：如果把内部递推式也写到 \(k=1\)，会得到\(x_2-(2-\alpha)x_1+x_0=0\Longrightarrow x_2=(2-\alpha)x_1-x_0=(1-\alpha)x_1,\)正好等价于端点方程 \(x_2=(1-\alpha)x_1\)（因为 \(x_0=x_1\)）。同理，在 \(k=m\) 处也会匹配 \(x_{m+1}=x_m\)。
因此，我们可以把所有 \(k=1,2,\dots,m\) 的方程统一写成：\(x_{k+1}-(2-\alpha)x_k+x_{k-1}=0, k=1,2,\dots,m,\)并带上边界条件\(x_0=x_1, x_{m+1}=x_m.\)
由\(x_{k+1}-(2-\alpha)x_k+x_{k-1}=0, k=1,2,\dots,m,\)这是一个常系数二阶齐次线性递推。按标准方法先尝试指数型解\(x_k=r^k (r\neq 0).\)
代入递推并约去公共因子 \(r^{k-1}\) 得到特征方程\(r^2-(2-\alpha)r+1=0.\)将其除以 \(r\) 可化为对称形式\(r+\dfrac{1}{r}=2-\alpha\)
由于 \(L_{P_m}\) 是实对称矩阵，\(\alpha\) 为实数；并且 \(L_{P_m}\) 是半正定矩阵，因此 \(\alpha\ge 0\)。
对路径图还可知谱落在 \([0,4]\) 内，于是 \(2-\alpha\in[-2,2]\)，因此可以用余弦作参数化：\(2-\alpha=2\cos\theta, \theta\in[0,\pi].\)这等价于\(\alpha=2-2\cos\theta.\)将 \(2-\alpha=2\cos\theta\) 代回到上面，得到\(r+\dfrac{1}{r}=2\cos\theta,\)从而特征方程的两根为\(r_1=e^{i\theta}, r_2=e^{-i\theta}.\)
因此递推的一般实解可以写成三角形式\(x_k=C\cos(k\theta)+D\sin(k\theta).\)
接下来利用边界条件 \(x_0=x_1\)。由上条方程得\(x_0=C, x_1=C\cos\theta+D\sin\theta,\)条件 \(x_0=x_1\) 等价于\(C=C\cos\theta+D\sin\theta\Longrightarrow D=C\cdot\dfrac{1-\cos\theta}{\sin\theta}=C\tan\dfrac{\theta}{2}.\)
将 \(D\) 消去并用恒等式\(\cos(k\theta)+\tan\frac{\theta}{2}\sin(k\theta)=\dfrac{\cos\left(k\theta-\frac{\theta}{2}\right)}{\cos\frac{\theta}{2}}\)得到\(x_k=C'\cos\left(\left(k-\dfrac12\right)\theta\right)\)
最后利用另一端边界条件 \(x_{m+1}=x_m\)。由上一条的方程有\(\cos\left(\left(m+\dfrac12\right)\theta\right)=\cos\left(\left(m-\dfrac12\right)\theta\right).\)两边相减并用恒等式 \(\cos u-\cos v=-2\sin\dfrac{u+v}{2}\sin\dfrac{u-v}{2}\)，得到\(-2\sin(m\theta)\sin\dfrac{\theta}{2}=0.\)因此要么 \(\theta=0\)（对应 \(\alpha=0\)），要么\(\sin(m\theta)=0\Longleftrightarrow \theta=\dfrac{\pi i}{m}, i=0,1,\dots,m-1.\)
代入\(\alpha=2-2\cos\theta.\)即得路径图 \(P_m\) 的拉普拉斯特征值为
\[ \alpha_i=2-2\cos\left(\frac{\pi i}{m}\right), i=0,1,\dots,m-1. \]
同理 \[ \beta_j = 2-2\cos\left(\frac{\pi j}{n}\right), j=0,1,\dots,n-1. \]
于是网格图的特征值为
\[ \lambda_{i,j} = \left(2-2\cos\frac{\pi i}{m}\right)+\left(2-2\cos\frac{\pi j}{n}\right) =4-2\cos\frac{\pi i}{m}-2\cos\frac{\pi j}{n}. \]
其中唯一的零特征值对应 \((i,j)=(0,0)\)。
代入生成树乘积公式 \(|V|=mn\)，得到闭式：\(\displaystyle{C(m,n)=\frac{1}{mn}\prod_{\substack{0\le i\le m-1\\0\le j\le n-1\\(i,j)\ne(0,0)}}\left(4-2\cos\frac{\pi i}{m}-2\cos\frac{\pi j}{n}\right).}\)
也可以用恒等式 \(2-2\cos\theta=4\sin^2(\theta/2)\) 写成
\[ C(m,n) =\frac{1}{mn}\prod_{\substack{0\le i\le m-1\\0\le j\le n-1\\(i,j)\ne(0,0)}} \left(4\sin^2\frac{\pi i}{2m}+4\sin^2\frac{\pi j}{2n}\right), \]
由于\(C(m,n)\) 极其巨大，直接乘会溢出。取常用对数：\(\displaystyle{L=\log_{10}C(m,n)=\sum_{(i,j)\ne(0,0)}\log_{10}\lambda_{i,j}-\log_{10}(mn).}\)
然后令\(E=\lfloor L\rfloor, M=10^{L-E},\)则\(C(m,n)=M\times 10^E.\)即可完成计算。
代码
1 | import math |