Project Euler 835

题目

Supernatural Triangles

A Pythagorean triangle is called supernatural if two of its three sides are consecutive integers.

Let \(S(N)\) be the sum of the perimeters of all distinct supernatural triangles with perimeters less than or equal to \(N\).

For example, \(S(100) = 258\) and \(S(10000) = 172004\).

Find \(S(10^{10^{10}})\). Give your answer modulo \(1234567891\).

解决方案

令\(n=10^{m}\)，其中\(m=10^{10}\)。

不失一般性，假设直角三角形的三条边分别为\(x,y,z\)，其中满足\(x< y< z\)。那么\(z\)是这个直角三角形的斜边。

那么题目中的超自然直角三角形分两种情况进行讨论。其中一种是\(z-y=1\)，另一种是\(y-x=1\)。

\(z-y=1\)

令\(z=y+1\)，那么\(x^2=(y+1)^2-y^2=2y+1\)。也就是说，当\(2y+1\)是一个平方数时，对应的\(x,y,y+1\)就构成了一个直角三角形。

令\(x=2k-1\)，其中\(k=\{2,3,4,\dots\}\)，那么有\(y=2k^2-2k,z=2k^2-2k+1\)，那么整个三角形的周长为\(x+y+z=2k(2k-1)\)。

令\(f(k)=2k(2k+1)\)，现在的任务是求最大的正整数\(k_0\)使得\(f(k_0)\le n\)。可以解得\(k_0\)为：

\[k_0=\left\lfloor\dfrac{1}{4} (\sqrt{4 n+1}+1)\right\rfloor\]

由于\(n\)是一个偶完全平方数，那么有\(k_0=\left\lfloor\dfrac{1}{4} (\sqrt{4 n+1}+1)\right\rfloor=\left\lfloor\dfrac{1}{4} (\sqrt{4 n}+1)\right\rfloor=\dfrac{\sqrt{n}}{2}\)。

也就是说，\(k_0=5\cdot 10^{m/2-1}\)。

最终这个子问题的答案为\(\displaystyle{S_0(k_0)=\sum_{k=2}^{k_0}f(k)=\dfrac{k_0(k_0+1)(4k_0-1)}{3}-2}\)。

\(y-x=1\)

令\(y=x+1\)，那么有\(x^2+(x+1)^2=z^2\)。化简后配方，得到\((2x+1)^2-2z^2=-1\)，令\(u=2x+1\)，那么得到负佩尔方程\((1)\)：

\[u^2-2z^2=-1\qquad(1)\]

并且，方程\((1)\)的每一对正整数解\((u,z)\)恰好对应着一个超自然直角三角形的周长\(c=u+z\)。

假设负佩尔方程\(u^2-Dz^2=-1\)存在最小特解为\((u_1,z_1)\)。那么，负佩尔方程的通解\((u_k,z_k)\)由以下式子导出。

\[u_{k+1}+z_{k+1}\sqrt{D}=(u_{k}+z_{k}\sqrt{D})(u_1+z_1\sqrt{D})^2\]

容易发现，\(u^2-2z^2=-1\)的最小特解为\((1,1)\)，代入上面的导出递推公式，得到通解递推公式：

\[\left \{\begin{aligned} & u_{k+1}=3u_k+4z_k\\ & z_{k+1}=2u_k+3z_k \end{aligned}\right. \]

令\(c_k=u_k+z_k\)，那么消去上面的递推式\(u_k\)，得到

\[\left \{\begin{aligned} & z_{k+1}=z_k+2s_k&\qquad(2)\\ & s_{k+1}=2z_k+5s_k&\qquad(3) \end{aligned}\right. \]

将\((2)\)代入\((3)\)，消去\(z_k\)，得到

\[s_{k+1}=2z_{k+1}+s_k\qquad(4)\]

由\((3)\)得：

\[s_{k+2}=2z_{k+1}+5s_{k+1}\qquad(5)\]

由\((4),(5)\)消去\(z\)，最终得到

\[s_{k+2}=6s_{k+1}-s_k\qquad(6)\]

为了方便后续计算，我们通过构造函数\(g(k)\)来代替\(s_k\)，其中\(g(k)=s_k\)，并且\(g\)在\(0\)下有意义。

构造\(g(0)=0,g(1)=2,g(k)=6g(k-1)-g(k-2)\)，现在的任务是求最大的正整数\(k_1\)使得\(g(k_1)\le n\)。不难计算得出，\(g(k)\)的通项公式为

\[g(k)=\dfrac{(3+2 \sqrt{2})^k-(3-2 \sqrt{2})^k}{2 \sqrt{2}}\]

假设\(n\)不在序列\(g\)中，当\(k\)非常大时，\(\dfrac{(3-2 \sqrt{2})^k}{2 \sqrt{2}}\)这一项的值非常小，因此不考虑。此时相当于求最大的\(k_1\)使得\(\dfrac{(3+2 \sqrt{2})^{k_1}}{2 \sqrt{2}}\le 10^{m}\)。两边对\(10\)取对数，最终得到

\[k_1=\left\lfloor \dfrac{m+\log_{10}(2\sqrt{2})}{\log_{10}(3+2\sqrt{2})}\right\rfloor\]

最终这个子问题的答案为\(\displaystyle{S_1(k_1)-g(1)=\sum_{k=1}^{k_1}g(k)-g(1)}\)。

接下来的问题是求解\(S_1(k_1)\)。将\(S_1(k)\)和\(g(k)\)写成矩阵形式为：

\[\begin{bmatrix} 0 & 1 & 0\\ -1 & 6 & 0\\ 0 & 1 & 1\\ \end{bmatrix} \begin{bmatrix} g(k-1)\\ g(k)\\ S_1(k-1) \end{bmatrix} = \begin{bmatrix} g(k)\\ g(k+1)\\ S_1(k) \end{bmatrix}\]

那么通过矩阵快速幂，可以以\(O(\log k_1)\)的时间复杂度内求解\(S_1(k_1)\)。

综合两个子问题，最终答案为\(S_0(k_0)+(S_1(k_1)-2)-\mathbf{12}\)，注意最后一项\(12\)是指\((3,4,5)\)的直角三角形，它在两个子问题中重复计算了一次。

OEIS相关信息：\(S_0\)在A002939中找到，\(S_1\)在A001542中找到。

代码

from math import log10, sqrt
from tools import mod_inverse

N = 10 ** 10
mod = 1234567891


def mul(a: list, b: list):
    return [[sum(a[i][k] * b[k][j] for k in range(len(b))) % mod for j in range(len(b[0]))] for i in range(len(a))]


n = 5 * pow(10, N // 2 - 1, mod) % mod
a1 = n * (n + 1) * (4 * n - 1) * mod_inverse(3, mod) % mod
sqt2 = sqrt(2)
m = int((N + log10(2 * sqt2)) / log10(3 + 2 * sqt2))
a = [[0, 2, 0]]
b = [[0, -1, 0], [1, 6, 1], [0, 0, 1]]
while m:
    if m & 1:
        a = mul(a, b)
    b = mul(b, b)
    m >>= 1
a2 = a[0][2]
ans = (a1 + a2 - 16) % mod
print(ans)