Project Euler 654
Project Euler 654
题目
Neighbourly Constraints
Let \(T(n, m)\) be the number of \(m\)-tuples of positive integers such that the sum of any two neighbouring elements of the tuple is \(\le n\).
For example, \(T(3, 4)=8\), via the following eight \(4\)-tuples:
\(\begin{aligned} &(1, 1, 1, 1)\\ &(1, 1, 1, 2)\\ &(1, 1, 2, 1)\\ &(1, 2, 1, 1)\\ &(1, 2, 1, 2)\\ &(2, 1, 1, 1)\\ &(2, 1, 1, 2)\\ &(2, 1, 2, 1) \end{aligned}\)
You are also given that \(T(5, 5)=246\), \(T(10, 10^{2}) \equiv 862820094 \pmod{1\,000\,000\,007}\) and \(T(10^2, 10) \equiv 782136797 \pmod{1\,000\,000\,007}\).
Find \(T(5000, 10^{12}) \bmod 1\,000\,000\,007\).
解决方案
令\(n=5000,m=10^{12}\)。考虑长度为 \(m\) 的正整数序列\((a_1,a_2,\dots,a_m), a_i\ge 1,\)满足相邻约束\(a_i+a_{i+1}\le n (1\le i<m).\)定义 \(T(n,m)\) 为满足约束的序列个数；并约定空序列\(T(n,0)=1.\)
注意到任意 \(a_i\le n-1\)，因此状态空间只需\(\{1,2,\dots,n-1\}.\)设\(d=n-1.\)
令\(f_t(y)=\#\{(a_1,\dots,a_t): a_t=y,a_i+a_{i+1}\le n\},1\le y\le d.\)由约束 \(x+y\le n\) 可得关于\(f\)的状态转移：
\[ f_t(y)= \left \{\begin{aligned} &1 & & \text{if}\quad t=1 \\ &\sum_{x=1}^{n-y} f_{t-1}(x) & & \text{else} \end{aligned}\right. \]
总数为\(\displaystyle{T(n,t)=\sum_{y=1}^{d} f_t(y).}\)因此有\(T(n,1)=d.\)
我们把向量记为\(\mathbf f_t=(f_t(1),f_t(2),\dots,f_t(d))^T.\)定义矩阵 \(M_n\in\{0,1\}^{d\times d}\)：
\[ (M_n)_{y,x}= \begin{cases} 1,& x+y\le n,\\ 0,& x+y>n. \end{cases} \]
那么转移就是\(\mathbf f_{t+1}=M_n\mathbf f_t,\mathbf f_1=\mathbf 1.\)并且\(T(n,t)=\mathbf 1^T \mathbf f_t=\mathbf 1^T M_n^{t-1}\mathbf 1.\)
设 \(M_n\) 的特征多项式为\(P_n(x)=\det(xI_d-M_n)= x^d + c_{d-1}x^{d-1}+\cdots+c_0.\)由 Caylay-Camilton定理 定理，可得\(P_n(M_n)=0,\)即\(M_n^d + c_{d-1}M_n^{d-1}+\cdots+c_0I=0.\) 两侧左乘 \(\mathbf 1^T M_n^{t-1}\)，右乘 \(\mathbf 1\)，得到标量递推：
\[ T(n,t+d)+c_{d-1}T(n,t+d-1)+\cdots+c_0T(n,t)=0(t\ge 0). \]
结论： 只要拿到 \(P_n(x)\) 的系数，再得到前 \(d\) 项\(T(n,0),T(n,1),\dots,T(n,d-1),\)就能用线性递推快速求第 \(m\) 项的方法把 \(m=M\) 做出来。
设\(A_t=T(n,t)(t\ge 0),\)那么有
\[ A_{t+d}+c_{d-1}A_{t+d-1}+\cdots+c_1A_{t+1}+c_0A_t=0(t\ge 0). \]
接下来关键在于：如何在 \(m\) 极大（例如 \(m=10^{12}\)）而 \(n\) 也很大（例如 \(n=5000\)）时，仍能高效求出 \(A_m\)。
通常若直接计算初值 \(A_0,A_1,\dots,A_{d-1}\) 会花费 \(O(n^2)\)，这是不可接受的；我们接下来用一种只用 \(O(n)\) 构造递推 + 初值的方式，并用快速幂在 \(O(\log m)\) 次多项式运算内完成求值。
在多项式环 \(\mathbb{F}[x]/(P_n(x))\) 中（\(\mathbb{F}\) 是模数域），由于\(x^d\equiv-(c_{d-1}x^{d-1}+\cdots+c_1x+c_0)\pmod{P_n(x)},\)因此任何 \(x^k\) 都能唯一化简成次数 \(<d\) 的多项式。
令\(R_m(x)\equiv x^{m+d-1}\pmod{P_n(x)}\)，且\(R_m\)的次数小于\(d\)，并写成\(R_m(x)=r_0+r_1x+\cdots+r_{d-1}x^{d-1}.\)
另一方面，由于 \(c_0\ne 0\)（事实上在偶数 \(n\) 时 \(c_0=\pm 1\)），可以在 \(x=0\) 处对\(\dfrac{1}{P_n(x)}\)展开得到：\(\displaystyle{\frac{1}{P_n(x)}=\sum_{k\ge 0}u_kx^k,}\)其中系数 \({u_k}\) 完全由恒等式 \(\displaystyle{P_n(x)\cdot \sum_{k\ge 0} u_kx^k\equiv 1}\) 决定。
此时可以写成\(x^{m+d-1}\)对\(P_n(x)\)做除法的式子：\(\dfrac{x^{m+d-1}}{P_n(x)}=Q_m(x)+\dfrac{R_m(x)}{P_n(x)}\)，其中 \(Q_m(x)\) 为某个多项式，不影响低阶系数，因此取 \(x^{d-1}\) 项系数有
\[ \left[x^{d-1}\right]\frac{x^{m+d-1}}{P_n(x)}=\left[x^{d-1}\right]\frac{R_m(x)}{P_n(x)}. \]
其中符号\(\left[x^t\right]P(x)\) 表示取幂级数\(P(x)\)中\(x^t\)那一项的系数。
而因为\(\displaystyle{\frac{R_m(x)}{P_n(x)}=\Bigl(\sum_{i=0}^{d-1}r_ix^i\Bigr)\Bigl(\sum_{k\ge0}u_kx^k\Bigr),}\)所以\(\displaystyle{\left[x^{d-1}\right]\frac{R_m(x)}{P_n(x)}=\sum_{i=0}^{d-1}r_iu_{d-1-i}.}\)
于是得到一个纯代数的求值公式：\(\displaystyle{A_m=\sum_{i=0}^{d-1}r_iu_{d-1-i}.}\)
现在问题被分解为三件事：
- 构造特征多项式 \(P_n(x)\)；
- 构造 \(u_0,u_1,\dots,u_{d-1}\)；
- 计算 \(R_m(x)\equiv x^{m+d-1}\pmod{P_n(x)}\)。
若 1、2 都能在 \(O(n)\) 完成，而 3 用幂模多项式做 \(O(\log m)\) 次乘法，则总体可达 \(O(n\log m)\) 级别。
接下来求解特征多项式 \(P_n(x)\) 的递推。以下只讨论 \(n\) 为偶数（也是目标规模 \(n=N\) 的情形）。令\(n=2s, d=n-1=2s-1.\)
对于矩阵\(M_n\)，可以对 \(\det(xI_d-M_n)\) 做一系列按反对角结构配对的消元，从而把行列式化为一个只依赖于 \(n\) 的二阶递推。最终得到：
当 \(n\) 为偶数时，\(P_n(x)\) 满足
\[ P_n(x)=(2x^2-1)P_{n-2}(x)-x^4P_{n-4}(x), \]
并且初值为\(P_2(x)=x-1, P_4(x)=x^3-2x^2-x+1.\)这意味着 \(P_n(x)\) 可以在 \(O(n)\) 次常数倍 + 平移运算中构造出来。
可见，当\(A\)的前面一些值确定后，我们还可以反向逆推到负数索引。即把 负指标初值 定义为\(A_{-k}=u_k(k\ge 0).\)
那么 \({A_t}_{t\in\mathbb Z}\) 就被同一个阶为 \(d\) 的递推唯一延拓到所有整数指标。\({A_t}_{t\in\mathbb Z}\) 同时满足：
- 在非负方向，它与原问题的 \(T(n,t)\) 一致；
- 在负方向，它由 \(1/P_n(x)\) 的幂级数系数给出；
- 正负两侧被同一条阶为 \(d\) 的递推统一连接起来（并且因为 \(c_0\) 可逆而具有唯一性）。
接下来解释在 \(k\le d-1=n-2\) 的范围内，\(u_k\) 的取值不再依赖 \(n\)，而是稳定到某个固定的模式。
对矩阵 \(xI_d-M_n\) 做配对消元（按其反对角的规则结构逐层约化）时，会自然产生一个逐层截断的有理函数列：
\[ F_0(z)=1, F_{j}(z)=z+\frac{1}{F_{j-1}(-z)}(j\ge 1). \]
为了研究 \(z=\infty\) 附近（等价于 \(x=0\) 附近）的展开，把变量换成 \(z=1/x\)，并引入在 \(x=0\) 处正规化的幂级数
\[ G_j(x)=\frac{1}{xF_j(1/x)}. \]
由定义立刻有 \(G_j(0)=1\)，且由递推式直接推出\(G_j(x)=\dfrac{1}{1-x^2G_{j-1}(-x)}.\)
注意右端每次出现的都是 \(x^2\)，因此每迭代一层，最低可能影响的次数至少上升 \(2\)。于是：要确定 \(G_j(x)\) 的展开到 \(x^{2k}\)，只需要 \(j\ge k\)，更深层的递推不会再影响这些低阶项。
因此当 \(n\) 足够大时，低阶系数稳定下来，极限幂级数\(\displaystyle{G(x)=\lim_{j\to\infty}G_j(x)}\)在每个有限阶上都存在，并满足同样的函数方程\(G(x)=\dfrac{1}{1-x^2G(-x)}.\)
把上式中的 \(x\) 替换为 \(-x\) 得到\(G(-x)=\dfrac{1}{1-x^2G(x)}.\)
将两式相乘并比较唯一性可知 \(G(x)=G(-x)\)，即 \(G\) 是偶函数。代回得到一元二次方程\(G(x)=\dfrac{1}{1-x^2G(x)}\Longleftrightarrow x^2G(x)^2-G(x)+1=0.\)
取满足 \(G(0)=1\) 的分支，得到\(\displaystyle{G(x)=\frac{1-\sqrt{1-4x^2}}{2x^2}=\sum_{k\ge 0}C_kx^{2k},}\)其中\(C_k=\dfrac{1}{k+1}\dbinom{2k}{k}\)是 Catalan 数，并等价满足一阶递推\(C_0=1,C_{k+1}=\dfrac{2(2k+1)}{k+2}C_k.\)
现在回到 \(u_k\)。由于 \(A_{-k}=u_k\) 是从 \(1/P_n(x)\) 的展开得到的，而低阶部分稳定到极限结构，最终得到在 \(k\le d-1=n-2\) 的范围内，有\(u_0=1,u_{2k}=0(k\ge 1),u_{2k+1}=C_k\ (k\ge 0).\)
因此\((u_0,u_1,u_2,u_3,\dots)=(1,1,0,1,0,2,0,5,0,14,\dots),\)并且只需 \(O(n)\) 的时间就能计算出前 \(d\) 项\(u_0,u_1,\dots,u_{d-1}.\)
在后续的求和公式中需要按从高到低的顺序使用这些负指标初值，因此将其倒置，定义\(v_i=u_{d-1-i}(0\le i\le d-1).\)
若把模多项式计算得到的余式写为\(\displaystyle{x^{m+d-1}\equiv \sum_{i=0}^{d-1} r_i x^i \pmod{P_n(x)},}\)则线性泛函 \(A_{t}=\bigl[x^{d-1}\bigr]\dfrac{x^{t+d-1}}{P_n(x)}\) 与“负指标初值”完全兼容，从而最终答案可写成
\[ A_m=\sum_{i=0}^{d-1}r_iv_i. \]
这一式把大指数 \(m\) 的困难全部压缩进了 \(r_i\)（多项式环幂运算），而把初值部分压缩进了一个几乎与 \(n\) 无关、并且呈 Catalan 结构的向量 \((v_i)\)。
剩下的任务是计算\(R_m(x)\equiv x^{m+d-1}\pmod{P_n(x)}.\)
这是标准的指数取模多项式问题：在商环 \(\mathbb{F}[x]/(P_n(x))\) 内用二进制快速幂即可。即不断进行\(F(x)\leftarrow F(x)^2\bmod P_n(x),\)若当前二进制位为 \(1\) 再乘一次 \(x\)：\(F(x)\leftarrow xF(x)\bmod P_n(x).\)
因为 \(\deg P_n=d=O(n)\)，整个过程只需要 \(O(\log m)\) 次多项式乘法 + 取模。
综上所述，令 \(n\) 为偶数，\(d=n-1\)，那么构造：
- 特征多项式 \(P_n(x)\)（次数 \(d\)，首项系数为 \(1\)）；
- 逆级数前缀\(\displaystyle{\frac{1}{P_n(x)}=\sum_{k\ge0}u_kx^k,}\)并取\(v_i=u_{d-1-i}(0\le i<d),\)其中\(u_0=1, u_{2k}=0, u_{2k+1}=C_k;\)
- 计算\(\displaystyle{R_m(x)\equiv x^{m+d-1}\pmod{P_n(x)}, R_m(x)=\sum_{i=0}^{d-1}r_ix^i.}\)
则最终
\[ T(n,m)=A_m=\sum_{i=0}^{d-1}r_iv_i. \]
在最后介绍一下\(P_n(x)\)递推式的计算过程。
我们现在用差分把 \(M_n\) 变成稀疏矩阵，考虑对行列式做两类不改变行列式值的初等变换：
- 对所有 \(1\le i\le d-1\)，做行变换\(R_i\leftarrow R_i-R_{i+1}.\)
- 对所有 \(1\le j\le d-1\)，做列变换\(C_j\leftarrow C_j-C_{j+1}.\)
这两组操作都等价于左乘/右乘一个对角线全为 \(1\) 的单位三角矩阵，因此行列式保持不变。
令\(A_n(x)=xI_d-M_n,\)做完上述行差分与列差分后得到新矩阵 \(B_n(x)\)，满足\(\det(B_n(x))=\det(A_n(x))=P_n(x).\)现在只需分析 \(B_n(x)\) 的形状。
记示性函数\(\chi(i,j)=\mathbf 1\{i+j\le n\}.\)对它做行差分与列差分等价于做离散二阶差分：\(\Delta_i\Delta_j\chi(i,j)=\chi(i,j)-\chi(i+1,j)-\chi(i,j+1)+\chi(i+1,j+1).\)
由于集合 \({(i,j):i+j\le n}\) 是一个由反对角线切出来的阶梯区域，它的二阶差分只可能出现在边界附近。直接分类即可得到：
- 当 \(i+j\le n-1\) 时，四项全为 \(1\)，所以差分为 \(0\)；
- 当 \(i+j\ge n+1\) 时，四项全为 \(0\)，所以差分为 \(0\)；
- 只有 \(i+j=n\) 或 \(i+j=n-1\) 时会出现跳变。
更精确地，有
\[ \Delta_i\Delta_j\chi(i,j)= \begin{cases} 1,& i+j=n-1,\\ -1,& i+j=n,\\ 0,& \text{otherwise}. \end{cases} \]
因此，经过“行差分+列差分”以后，\(-M_n\) 的贡献只会落在两条紧邻的反对角线上（\(i+j=n-1\) 与 \(i+j=n\)），矩阵一下子变得非常稀疏。
同理，\(xI_d\) 在做差分后会产生一个带宽为 \(2\)的固定带状结构（只出现在 \(|i-j|\le 1\) 的附近），所以 \(B_n(x)\) 变成一个带宽不超过 \(2\) 的稀疏矩阵：每一行只有常数个非零元。
接下来做一个纯置换：把行与列同时按下面顺序重排（同一个置换作用到行与列，因此行列式仍不变）：\((d,1),(d-1,2),(d-2,3),\ldots,(m+1,m-1),m.\)
也就是把指标从两端向中间拉链式排列。在这个顺序下，矩阵 \(B_n(x)\) 会变成一个块三对角矩阵 \(K_m(x)\)：对角线上是 \(2\times 2\) 的小块，相邻块之间只有固定的耦合，除此之外全为 \(0\)。更具体地，它的左上角 \(4\times4\) 主块形如
\[ \begin{pmatrix} x & -1 & -x & 0\\ -1& 2x & 1 & -x\\ -x& 1 & 2x & -1\\ 0 & -x & -1 & 2x \end{pmatrix}, \]
把它按 \(2\times2\) 分块写成\(\begin{pmatrix}B_0 & C\\C^\ast & B_0\end{pmatrix},\)其中\(B_0=\begin{pmatrix} x & -1\\ -1& 2x \end{pmatrix}, C= \begin{pmatrix} -x & 0\\ 1 & -x \end{pmatrix}, C^\ast= \begin{pmatrix} -x & 1\\ 0 & -x \end{pmatrix}.\)
而整个 \(K_m(x)\) 由同样的块模式向右下不断延伸（只在最末端因为维度是奇数会出现一个 \(1\times 1\) 的收尾块，但这不会影响下面的二阶递推结构）。
因为这些变换都不改变行列式，所以\(P_{2m}(x)=\det(K_m(x)).\)
令\(D_m=\det(K_m(x))=P_{2m}(x).\) 由于 \(K_m(x)\) 是块三对角的，并且第一块只通过 \(C,C^\ast\) 与后面的部分相连，按块展开（等价于对前两行列做消元/拉普拉斯展开）会出现两类贡献：
- 完全落在 \(B_0\) 内部的选取，其系数就是 \(\det(B_0)\)，剩下的余子式正是 \(K_{m-1}(x)\)，贡献为 \(\det(B_0)D_{m-1}\)；
- 通过 \(C,C^\ast\) “跨到下一块” 的选取，其系数是 \(\det(C)\det(C^\ast)\)，剩下的余子式变成 \(K_{m-2}(x)\)，贡献为 \(-\det(C)\det(C^\ast)D_{m-2}\)（符号来自行列式展开的交错性）。
因此得到\(D_m=\det(B_0)D_{m-1}-\det(C)\det(C^\ast)D_{m-2}.\)现在计算这些 \(2\times2\) 行列式即可：
\[ \det(B_0)=x\cdot 2x-(-1)\cdot(-1)=2x^2-1,\det(C)=(-x)(-x)-0\cdot 1=x^2,\det(C^\ast)=x^2. \]
代回得到\(D_m=(2x^2-1)D_{m-1}-x^4D_{m-2}.\)
把 \(D_m=P_{2m}(x)\) 换回 \(n=2m\) 的记号，就得到
\[ P_n(x)=(2x^2-1)P_{n-2}(x)-x^4P_{n-4}(x), \]
- \(n=2\) 时，\(d=1\)，\(M_2=[1]\)，所以\(P_2(x)=\det([x]-[1])=x-1.\)
- \(n=4\) 时，\(d=3\)， \(M_4=\begin{pmatrix}1&1&1\\1&1&0\\1&0&0\end{pmatrix}, xI_3-M_4= \begin{pmatrix}x-1&-1&-1\\-1&x-1&0\\-1&0&x\end{pmatrix}.\)直接计算行列式得到\(P_4(x)=x^3-2x^2-x+1.\)
代码
1 | import numpy as np |