Project Euler 332

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

Project Euler 332

题目

Spherical triangles

A spherical triangle is a figure formed on the surface of a sphere by three great circular arcs intersecting pairwise in three vertices.

Let \(C(r)\) be the sphere with the centre \((0,0,0)\) and radius r.

Let \(Z(r)\) be the set of points on the surface of \(C(r)\) with integer coordinates.

Let \(T(r)\) be the set of spherical triangles with vertices in \(Z(r)\). Degenerate spherical triangles, formed by three points on the same great arc, are not included in \(T(r)\).

Let \(A(r)\) be the area of the smallest spherical triangle in \(T(r)\).

For example \(A(14)\) is \(3.294040\) rounded to six decimal places.

Find \(\sum \limits_{r = 1}^{50} {A(r)}\). Give your answer rounded to six decimal places.

解决方案

引理：若 \(r\) 为偶数，则 \(Z(r)\) 中任意点 \((x,y,z)\) 的三个坐标都为偶数。

证明：当 \(r\) 为偶数时，\(r^2\equiv 0\pmod 4\)。对任意整数，平方模 \(4\) 只可能是 \(0\) 或 \(1\)。若 \(x\) 为奇数则 \(x^2\equiv 1\pmod 4\)，同理奇数平方都为 \(1\)。 在三个数平方之和中，若存在奇数平方项，则和模 \(4\) 为 \(1,2,\) 或 \(3\)，不可能为 \(0\)。因此必须\(x^2\equiv y^2\equiv z^2\equiv 0\pmod 4,\)即 \(x,y,z\) 全为偶数。证毕。

于是对 \(r=2m,Z(2m)=\{(2x,2y,2z):(x,y,z)\in Z(m)\}.\)

把点除以半径得到单位向量 \((x,y,z)/r\)，可见 \(Z(2m)\) 与 \(Z(m)\) 给出的单位方向集合完全相同。球面三角形的角度（球面边长、角、球面盈余）只取决于单位方向，而面积会额外乘上 \(r^2\)，因此\(A(2m)=4A(m), A(4m)=16A(m),\ldots\)

所以只要计算所有奇数半径的 \(A(r)\)，就能把 \(r,2r,4r,\dots\le 50\) 的贡献用 \(4^k\) 直接加进去。

接下来进行球面三角形面积的向量公式推导。取球面三角形三个顶点向量为 \(a,b,c\in\mathbb R^3\)，满足\(\lVert a\rVert=\lVert b\rVert=\lVert c\rVert=r.\)令单位向量\(u=\dfrac{a}{r}, v=\dfrac{b}{r}, w=\dfrac{c}{r}.\)

球面三角形的面积为\(S=r^2E,\)其中球面盈余（定义为球面曲率造成的角和超额）\(E=\alpha+\beta+\gamma-\pi,\)\(\alpha,\beta,\gamma\) 为球面三角形三个顶点角。

设三条球面边（中心角）分别为\(A=\arccos(v\cdot w), B=\arccos(w\cdot u), C=\arccos(u\cdot v).\)其中 \(u\cdot v\) 表示点积。由球面三角学可得（将 L’Huilier 公式化简到只含 \(\cos A,\cos B,\cos C\) 的形态）：

\[\tan^2\frac E2=\frac{1+2\cos A\cos B\cos C-\cos^2A-\cos^2B-\cos^2C}{(1+\cos A+\cos B+\cos C)^2}.\]

取正根并用 \(\operatorname{atan2}\) 写成更稳定的形式：

\[\frac E2=\operatorname{atan2}(\sqrt{1+2\cos A\cos B\cos C-\cos^2A-\cos^2B-\cos^2C},1+\cos A+\cos B+\cos C).\]

注意到\(\cos A=v\cdot w, \cos B=w\cdot u, \cos C=u\cdot v.\)

考虑单位向量的 Gram 矩阵

\[G=\begin{pmatrix} u\cdot u & u\cdot v & u\cdot w\\ v\cdot u & v\cdot v & v\cdot w\\ w\cdot u & w\cdot v & w\cdot w \end{pmatrix} =\begin{pmatrix} 1 & u\cdot v & u\cdot w\\ u\cdot v & 1 & v\cdot w\\ u\cdot w & v\cdot w & 1 \end{pmatrix}.\]

线性代数恒等式给出\(\det(G)=(u\cdot(v\times w))^2,\)其中 \(v\times w\) 为叉积。把 \(\det(G)\) 展开，恰好得到

\[\det(G)=1+2(u\cdot v)(v\cdot w)(w\cdot u)-(u\cdot v)^2-(v\cdot w)^2-(w\cdot u)^2.\]

这正是上一节分子中的平方项，因此\(\sqrt{1+2\cos A\cos B\cos C-\cos^2A-\cos^2B-\cos^2C}=|u\cdot(v\times w)|.\)

把分母写回点积和：\(1+\cos A+\cos B+\cos C = 1+u\cdot v+v\cdot w+w\cdot u.\)

于是得到非常紧凑且数值稳定的公式：

\[E=2\operatorname{atan2}(|u\cdot(v\times w)|,1+u\cdot v+v\cdot w+w\cdot u).\]

代入 \(u=a/r\) 等并整理比例：

\[|u\cdot(v\times w)| =\left|\frac a r\cdot\left(\frac b r\times\frac c r\right)\right| =\frac{|a\cdot(b\times c)|}{r^3},\]

\[1+u\cdot v+v\cdot w+w\cdot u =1+\frac{a\cdot b+b\cdot c+c\cdot a}{r^2} =\frac{r^2+a\cdot b+b\cdot c+c\cdot a}{r^2}.\]

因此\(E=2\operatorname{atan2}(|a\cdot(b\times c)|,r(r^2+a\cdot b+b\cdot c+c\cdot a)),\)面积即\(S=r^2E.\)

这里 \(a\cdot(b\times c)\) 与各点积都可以用整数运算精确计算。

对于退化情形，三点在同一条大圆弧上 \(\Longleftrightarrow\) 三点与原点共面 \(\Longleftrightarrow\) 向量 \(a,b,c\) 线性相关 \(\Longleftrightarrow a\cdot(b\times c)=0.\)因此只要检查三重积是否为 \(0\)，即可精确排除退化情形。

对固定 \(r\)，由\(E=2\operatorname{atan2}(\Delta,rN),\Delta=|a\cdot(b\times c)|,N=r^2+a\cdot b+b\cdot c+c\cdot a\)可知在 \(\Delta>0\) 且 \(N>0\) 时，\(E\) 随 \(\dfrac{N}{r}\)（等价于 \(\dfrac{N}{\Delta}\)）增大而减小。于是最小面积等价于在所有非退化三角形中最大化\(\dfrac{N}{\Delta}=\dfrac{r^2+a\cdot b+b\cdot c+c\cdot a}{|a\cdot(b\times c)|}.\)比较两个分数可用交叉相乘完全避免浮点误差： \(\dfrac{N_1}{\Delta_1}>\dfrac{N_2}{\Delta_2}\iff N_1\Delta_2>N_2\Delta_1.\)

另外，方程 \(x^2+y^2+z^2=r^2\) 对坐标置换与符号翻转对称（共 \(48\) 种正交变换），面积表达式只依赖点积与三重积的绝对值，因此不妨把第一个顶点 \(a\) 限制在一个规范域内，例如\(0\le x\le y\le z.\)这样可以显著减少枚举量，而不影响最小值结果。

结合第 1 节偶数半径缩放，只需对奇数 \(r\le 50\) 计算一次 \(A(r)\)，并按 \(r,2r,4r,\dots\) 把贡献累加即可。

代码

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
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
from tools import int_sqrt
from math import atan2
from itertools import permutations, product

N = 50


def points_on_sphere(r: int):
    """
    枚举球面 C(r) 上的整点：x^2 + y^2 + z^2 = r^2。

    返回两份列表：
    - base: 仅在规范域 0<=x<=y<=z 内的解 (x,y,z)，用于当作锚点 a
            （利用符号翻转与坐标置换对称性，把 48 倍的等价点折叠起来）
    - full: 将 base 展开符号翻转与坐标置换后的全集 Z(r)，用于枚举 b,c
    """
    r2 = r * r

    # base 只收集非负且排序的三元组 0<=x<=y<=z，避免同一轨道重复出现
    base = []
    for x in range(r + 1):
        x2 = x * x
        for y in range(x, r + 1):
            # rem = r^2 - x^2 - y^2，若 rem<0 则更大的 y 只会更负，可直接 break
            rem = r2 - x2 - y * y
            if rem < 0:
                break

            # z = ceil/floor? 这里要的就是整数平方根（floor），并检查是否恰好平方
            z = int_sqrt(rem)

            # 维持排序 x<=y<=z；若 z<y 则该 (x,y) 不可能得到规范解
            if z < y:
                continue
            if z * z == rem:
                base.append((x, y, z))

    # full：把 base 中每个 (x,y,z) 通过
    #   - 坐标置换（最多 6 种）
    #   - 三个坐标独立符号翻转（2^3=8 种）
    # 展开成全集 Z(r)
    #
    # 注意：当 x,y,z 有重复值时 permutations 会产生重复，需要用 set 去重。
    pts = set()
    for x, y, z in base:
        for px, py, pz in set(permutations((x, y, z), 3)):
            for sx, sy, sz in product((1, -1), repeat=3):
                pts.add((sx * px, sy * py, sz * pz))

    return base, list(pts)


def min_area_for_odd_r(r: int) -> float:
    """
    计算奇数半径 r 的 A(r)。

    几何公式（向量均以原点为球心，|a|=|b|=|c|=r）：

    设
      det = | a · (b × c) |
      num = r^2 + a·b + b·c + c·a

    则该球面三角形对应的球面盈余（也等于立体角）为
      E = 2 * atan2(det, r * num)

    面积为
      S = r^2 * E

    退化情形：
      三点与原点共面  <=>  det = 0
      这正对应“三点在同一条大圆弧上”，题目要求排除。

    由于 r 固定、atan2 对正比值单调递增，最小化 E 等价于最小化 det/(r*num)，
    也等价于最大化 num/det（纯整数比较可避免浮点误差）。
    """
    base, Z = points_on_sphere(r)
    r2 = r * r

    # best_num/best_det 维护当前发现的最大 num/det（对应最小面积）
    # 初始化为 0，表示尚未找到任何合法三角形
    best_num = 0
    best_det = 1

    m = len(Z)

    # 关键剪枝：只让 a 在规范域 base 内枚举（约缩小到 1/48），
    # b,c 仍用全集 Z(r) 枚举，这样不会漏掉任何三角形：
    # 任意三角形可通过对称变换把其中一个顶点搬到规范域。
    for a in base:
        ax, ay, az = a

        # 枚举 b（固定后再枚举 c），dab = a·b 可提前算一次
        for j in range(m):
            bx, by, bz = Z[j]
            dab = ax * bx + ay * by + az * bz  # a·b

            # 枚举 c，要求 j<k 避免重复组合（b,c 无序即可）
            for k in range(j + 1, m):
                cx, cy, cz = Z[k]

                # det = |a · (b × c)|，完全用整数运算：
                # b×c = (by*cz-bz*cy, bz*cx-bx*cz, bx*cy-by*cx)
                det = abs(
                    ax * (by * cz - bz * cy)
                    + ay * (bz * cx - bx * cz)
                    + az * (bx * cy - by * cx)
                )
                if det == 0:
                    # 退化：a,b,c 与原点共面 => 在同一大圆上 => 不计入 T(r)
                    continue

                # num = r^2 + a·b + b·c + c·a（同样是整数）
                # 其中 b·c、c·a 都现场计算
                num = (
                    r2
                    + dab
                    + (bx * cx + by * cy + bz * cz)          # b·c
                    + (cx * ax + cy * ay + cz * az)          # c·a
                )

                # 对于最小面积候选，num 一定应为正；
                # num<=0 对应 atan2(det, r*num) 接近 π/2 或更大，面积很大，不可能最小
                if num <= 0:
                    continue

                # 比较 num/det 是否更大：num/det > best_num/best_det
                # 用交叉相乘避免浮点误差：num*best_det > best_num*det
                if num * best_det > best_num * det:
                    best_num = num
                    best_det = det

    # 找到最佳比值后，计算对应球面盈余 E 与面积 S：
    # E = 2*atan2(det, r*num),  S = r^2*E
    E = 2.0 * atan2(best_det, r * best_num)
    return r2 * E


def solve(R) -> float:
    """
    计算 sum_{r=1..R} A(r)。

    关键优化：只计算奇数半径。

    证明思路（只用模 4）：
      若 r 为偶数，则 r^2 ≡ 0 (mod 4)。
      任意整数平方满足 x^2 ≡ 0 或 1 (mod 4)。
      若 x,y,z 中存在奇数，则对应平方 ≡1 (mod 4)，
      三个平方之和在 (mod 4) 下只能是 0,1,2,3 中的某个。
      想要 x^2+y^2+z^2 ≡ 0 (mod 4) 且只有三个项，唯一可能是三项全 ≡0，
      即 x,y,z 全为偶数。

    因此当 r 为偶数时，Z(r) 上所有点坐标都为偶数，可整体除以 2 映射到 Z(r/2)。
    这表明：
      - 方向（单位向量）不变，球面三角形的角度与球面盈余 E 不变；
      - 半径从 r/2 变为 r，面积按 r^2 缩放，故 A(r) = 4*A(r/2)。

    所以对每个奇数 r，只需求一次 A(r)，再把 2r,4r,8r,... 的贡献按 4^k 累加。
    """
    total = 0.0

    # 只枚举奇数半径 r
    for r in range(1, R + 1, 2):
        a = min_area_for_odd_r(r)

        # 把 r,2r,4r,... <=R 的 A 值累加：
        # A(2^k r) = 4^k * A(r)
        rr = r
        mul = 1.0
        while rr <= R:
            total += a * mul
            rr *= 2
            mul *= 4.0

    return total


ans = solve(N)
print(f"{ans:.6f}")