Project Euler 390

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

Project Euler 390

题目

Triangles with non rational sides and integral area

Consider the triangle with sides \(\sqrt 5,\sqrt {65}\) and \(\sqrt {68}\). It can be shown that this triangle has area \(9\).

\(S(n)\) is the sum of the areas of all triangles with sides \(\sqrt{1+b^2}\), \(\sqrt {1+c^2}\) and \(\sqrt{b^2+c^2}\,\) (for positive integers \(b\) and \(c\)) that have an integral area not exceeding \(n\).

The example triangle has \(b=2\) and \(c=8\).

\(S(10^6)=18018206\).

Find \(S(10^{10})\).

解决方案

Project Euler 390 题解

记三边\(a=\sqrt{1+b^2},d=\sqrt{1+c^2},e=\sqrt{b^2+c^2},\)面积为 \(S\)。使用海伦公式的等价形式（对任意三角形成立）：

\[ 16S^2 =2(a^2d^2+d^2e^2+e^2a^2)-(a^4+d^4+e^4). \]

代入\(a^2=1+b^2, d^2=1+c^2, e^2=b^2+c^2\)，逐项计算：

  • \(a^2d^2=(1+b^2)(1+c^2)=1+b^2+c^2+b^2c^2\)
  • \(d^2e^2=(1+c^2)(b^2+c^2)=b^2+c^2+b^2c^2+c^4\)
  • \(e^2a^2=(b^2+c^2)(1+b^2)=b^2+c^2+b^4+b^2c^2\)

因此\(a^2d^2+d^2e^2+e^2a^2=1+3b^2+3c^2+3b^2c^2+b^4+c^4.\)再算，得到

  • \(a^4=(1+b^2)^2=1+2b^2+b^4,\)
  • \(d^4=1+2c^2+c^4,\)
  • \(e^4=(b^2+c^2)^2=b^4+2b^2c^2+c^4,\)

所以\(a^4+d^4+e^4=2+2b^2+2c^2+2b^4+2c^4+2b^2c^2.\)代回，得到：

\[ \begin{aligned} 16S^2 &=2(1+3b^2+3c^2+3b^2c^2+b^4+c^4)-(2+2b^2+2c^2+2b^4+2c^4+2b^2c^2)\\ &=4b^2+4c^2+4b^2c^2. \end{aligned} \]

得到核心恒等式：

\[ 4S^2=b^2c^2+b^2+c^2. \]

令\(n=2S\)，可见\(n\in \mathbb{Z}\)，则上式等价为：

\[ n^2=b^2c^2+b^2+c^2. \]

面积为整数并且不超过 \(N\) 等价于\(n\le 2N.\)

因为 \(n=2S\) 为偶数，所以 \(n^2\equiv 0\pmod 4\)。观察\(n^2=b^2c^2+b^2+c^2 \pmod 4.\)若 \(b\) 为奇数，则 \(b^2\equiv 1\pmod4\)：

  • 若 \(c\) 为偶数，则 \(c^2\equiv 0\pmod4\)，右侧 \(\equiv 0+1+0\equiv 1\pmod4\)，矛盾；
  • 若 \(c\) 为奇数，则右侧 \(\equiv 1+1+1\equiv 3\pmod4\)，矛盾。

同理若 \(c\) 为奇数也矛盾，因此 \(b,c\)必为偶数。

由\(n^2=b^2c^2+b^2+c^2\)整理得到一个对称形式：\(n^2=(b^2+1)(c^2+1)-1\Longrightarrow n^2+1=(b^2+1)(c^2+1).\)更关键的是，把它看成对固定 \(b\) 的二元丢番图：\(n^2-(b^2+1)c^2=b^2.\)记\(D=b^2+1,\)则方程为

\[ n^2-Dc^2=b^2.\tag{1} \]

这就是典型的广义佩尔方程（右边不是 \(1\) 而是 \(b^2\)）。

考虑标准佩尔方程\(u^2-Dv^2=1.\)当 \(D=b^2+1\) 时，有一个显式整数解：\(u_1=2b^2+1, v_1=2b.\)可验证：\((2b^2+1)^2-(b^2+1)(2b)^2=4b^4+4b^2+1-(4b^4+4b^2)=1.\)

因此 \(u_1+v_1\sqrt{D}\) 是该二次整数环中的一个单位元。

若 \((n,c)\) 满足 \((1)\)，则在环 \(\mathbb Z[\sqrt D]\) 中\((n+c\sqrt D)(u_1+v_1\sqrt D)\)

仍满足同样的范数条件，从而得到新的解 \((n',c')\)：\(n'=u_1n+v_1Dc=(2b^2+1)n+2b(b^2+1)c,c'=v_1n+u_1c=2bn+(2b^2+1)c.\)即递推公式：

\[ \begin{aligned} n_{k+1}&=(2b^2+1)n_k+2b(b^2+1)c_k,\\ c_{k+1}&=2bn_k+(2b^2+1)c_k. \end{aligned} \tag{2} \]

方程 \((1)\) 显然有平凡解\((n,c)=(b,0),\)因为 \(b^2-D\cdot 0=b^2\)。虽然 \(c=0\) 不对应题目需要的三角形，但它可以作为递推的起点；用 \((2)\) 推一步立刻得到 \(c>0\)，从而产生合法三角形解。

原问题关于 \(b,c\) 对称，只与平方有关。所以交换后，\((b,c,n)\Rightarrow (c,b,n)\) 仍是解；变号后，\((b,c,n)\Rightarrow (b,-c,n)\) 仍是解（因为只出现 \(c^2\)）。实践上可以这样做：

  1. 对每个偶数 \(b_0\)，从种子 \((b,c,n)=(b_0,0,b_0)\) 出发；
  2. 沿着固定 \(b\) 的递推 \((2)\) 不断生成更大的 \((n,c)\)；
  3. 每得到一个合法解（\(n\le 2N\) 且 \(c\ne 0\)），就把它对应的面积 \(n/2\) 加入答案；
  4. 同时把该解当作“换根”，把新的种子加入栈：\((b,c,n)\Rightarrow (|c|,|b|,n), (|c|,-|b|,n)\)用来覆盖其它等价类。

为保证不重计数，只需用集合去重：把三角形用规范化键\((\min(|b|,|c|),\max(|b|,|c|),n)\)表示即可。

从种子 \((n,c)=(b,0)\) 用 \((2)\) 推一步：\(n_1=(2b^2+1)b=2b^3+b,\)所以第一批非退化解的面积为\(S=\dfrac{n_1}{2}=b^3+\dfrac b2>b^3.\)

而递推会使 \(n\) 单调快速增大，因此若 \(b^3>N\)，该分支不可能产生面积 \(\le N\) 的三角形。因此只需\(b\le \lfloor N^{1/3}\rfloor.\)可见枚举非常轻量。

代码

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
N = 10 ** 10

def solve(N: int) -> int:
    lim = 2 * N
    b_max = int(N ** (1 / 3)) + 2
    seen = set()
    ans = 0

    for b0 in range(2, b_max + 1, 2):
        stack = [(b0, 0, b0)]
        while stack:
            b, c, n = stack.pop()
            while True:
                n1 = (2 * b * b + 1) * n + 2 * b * (b * b + 1) * c
                c1 = 2 * b * n + (2 * b * b + 1) * c
                n, c = n1, c1

                if n > lim:
                    break
                if c == 0 or b == 0:
                    continue

                b = abs(b)
                c = abs(c)
                x, y = (b, c) if b <= c else (c, b)
                key = (x, y, n)
                if key in seen:
                    continue
                seen.add(key)

                ans += n // 2
                stack.append((c, b, n))
                stack.append((c, -b, n))

    return ans

print(solve(N))