Project Euler 373

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

Project Euler 373

题目

Circumscribed Circles

Every triangle has a circumscribed circle that goes through the three vertices. Consider all integer sided triangles for which the radius of the circumscribed circle is integral as well.

Let \(S(n)\) be the sum of the radii of the circumscribed circles of all such triangles for which the radius does not exceed \(n\).

\(S(100)=4950\) and \(S(1200)=1653605\).

Find \(S(10^7)\).

解决方案

设三边为 \(a,b,c\)，面积为 \(\Delta\)，外接圆半径为 \(R\)。那么有：\(R=\dfrac{abc}{4\Delta}.\)因为 \(a,b,c,R\) 都是整数，所以\(4\Delta=\dfrac{abc}{R}\in\mathbb Z.\)于是\(16\Delta^2=(4\Delta)^2\)必然是一个完全平方数。

另一方面，海伦公式给出\(16\Delta^2=(a+b+c)(a+b-c)(b+c-a)(c+a-b).\)记四个因子\(A=a+b+c,B=a+b-c,C=b+c-a,D=c+a-b,\)则非退化三角形等价于 \(A,B,C,D>0\)，且有恒等式\(A=B+C+D.\)并且\(ABCD\)是完全平方数。

由于 \(ABCD\) 是平方数，可以把 \(A,B,C,D\) 写成“三对共享因子”的形式（本质是把每个素因子的指数按“出现两次”的方式分配）：\(A=xzu, B=xtv, C=yzv, D=ytu,\)其中 \(x,y,z,t,u,v\) 为正整数。（当然可再附加一些互素条件保证“原始解”不重计，但在计数结论里最终会被吸收。）

此时\(ABCD=(xzu)(xtv)(yzv)(ytu)=(xyztuv)^2\)确实是平方数，并且立刻得到\(4\Delta=\sqrt{ABCD}=xyztuv\in\mathbb Z.\)

把 \(A=B+C+D\) 代入上式参数化：\(xzu=xtv+yzv+ytu.\)整理为\(u(xz-yt)=v(xt+yz).\)

定义三个高斯整数：\(X=x+iy, Z=z+it, U=u+iv.\)则\(XZ=(xz-yt)+i(xt+yz).\)上一步的等式 \(u(xz-yt)=v(xt+yz)\) 等价于 \(XZU\) 是纯虚数。这是非常典型的高斯整数结构：让一个乘积的实部为 \(0\)。

一个通用的构造是取任意高斯整数 \(P,Q,W\)，令\(X=iP\overline{Q}, Z=Q\overline{W}, U=W\overline{P}.\)那么\(XZU=iP\overline{Q}\cdot Q\overline{W}\cdot W\overline{P}=i\cdot N(P)N(Q)N(W)\) 确实是纯虚数（\(N(\cdot)\) 表示高斯整数的范数 \(N(a+ib)=a^2+b^2\)）。

也就是说：所有满足条件的三角形都能用三个高斯整数 \(P,Q,W\) 的范数与辐角组合来构造出来。

在上述构造下，可以把三角形的三边与面积都写成由 \(P,Q,W\) 的范数与叉积（实部/虚部）组成的整数表达式；最终外接圆半径会简化成一个非常干净的形式：\(8R = N(P)N(Q)N(W).\)这个结论的推导过程如下：

设\(P=p_1+ip_2, Q=q_1+iq_2, W=w_1+iw_2(p_1,p_2,q_1,q_2,w_1,w_2\in\mathbb Z),\)那么有\(X=x+iy, Z=z+it, U=u+iv,\)则逐项展开可得（全是整数）：

$ x=p_1q_2-p_2q_1, y=p_1q_1+p_2q_2,z=q_1w_1+q_2w_2, t=q_2w_1-q_1w_2,u=w_1p_1+w_2p_2, v=w_2p_1-w_1p_2. $

同时范数满足乘法性：\(N(X)=x^2+y^2=N(P)N(Q),N(Z)=z^2+t^2=N(Q)N(W),N(U)=u^2+v^2=N(W)N(P).\)

根据上文，我们可以得到反解为\(a=\dfrac{B+D}{2}, b=\dfrac{B+C}{2}, c=\dfrac{C+D}{2}.\)代入 \(B,C,D\) 得到三边：\(a=\dfrac{xtv+ytu}{2}=\dfrac{t(xv+yu)}{2},b=\dfrac{xtv+yzv}{2}=\dfrac{v(xt+yz)}{2},c=\dfrac{yzv+ytu}{2}=\dfrac{y(zv+tu)}{2}.\)

而在当前构造 \(X=iP\overline Q,Z=Q\overline W,U=W\overline P\) 下，它们还能进一步简化成非常干净的形式：

先注意\(XU=iP\overline Q\cdot W\overline P=iN(P)W\overline Q,\)因此\(xv+yu=N(P)\cdot z.\)同理\(XZ=iP\overline Q\cdot Q\overline W=iN(Q)P\overline W\Rightarrow xt+yz=N(Q)\cdot u,\) 以及\(ZU=Q\overline W\cdot W\overline P=N(W)Q\overline P\Rightarrow zv+tu=N(W)\cdot x.\)

所以三边可以写成：\(a=\dfrac{t\cdot N(P)\cdot z}{2},b=\dfrac{v\cdot N(Q)\cdot u}{2},c=\dfrac{y\cdot N(W)\cdot x}{2}.\)

可见，若你已选取 \(P,Q,W\) 使得 \(x,y,z,t,u,v>0\)，则可直接写\(\Delta=\dfrac{xyztuv}{4}.\)

由简化后的三边式，有\(abc=\dfrac{tN(P)z}{2}\cdot\dfrac{vN(Q)u}{2}\cdot\dfrac{yN(W)x}{2} =\dfrac{xyztuv\cdot N(P)N(Q)N(W)}{8}\).

因此，基于外接圆的面积公式，有：

\[ 8R=N(P)N(Q)N(W). \]

所以最终 \(R\) 的本质来自和平方数的乘积。

也就是说，原问题等价于以下计数-求和问题：对每个正整数 \(R\le n\)，计数满足 \(8R=N(P)N(Q)N(W)\) 的高斯整数三元组 \((P,Q,W)\) 所诱导出的整数三角形 \((a,b,c)\) 的个数（按三角形同构去重），并计算加权和

\[S(n)=\sum_{R\le n} R\cdot T(R),\]

其中 \(T(R)\) 表示外接圆半径恰为 \(R\) 的此类三角形的数量（通常以无序边三元组 \({a,b,c}\) 为单位计数）。

在高斯整数环 \(\mathbb Z[i]\) 中，满足 \(p\equiv3\pmod4\) 的素数不会分裂（仍为不可约）。因此若一个整数能写成两平方和（也即某个高斯整数的范数），那么它的分解中所有 \(p\equiv3\pmod4\) 的素因子都只能以偶次幂出现。由此可见：

  • 在统计“本质不同的三角形形状”（原始解的角度/组合类型）时，真正起作用的仅是 \(R\) 中满足 \(p\equiv1\pmod4\) 的素因子指数；
  • 其余因子，尤其是 \(p\equiv3\pmod4\) 的部分，只会以“平方因子”的方式进入，从效果上等同于整体缩放，不会改变同一“基本半径类型”所对应的形状数量。

于是我们把半径分成：\(R = r\cdot k,\)其中 基本半径 \(r\) 只含 \(p\equiv 1\pmod 4\) 的素因子（允许高次），而 \(k\) 代表整体缩放（把三边、面积、半径都乘同一个整数，仍保持整数性）。

因此计算 \(S(n)\) 的标准做法是：

  1. 枚举所有基本半径 \(r\le n\)；
  2. 对每个 \(r\)，算出其对应“原始三角形个数” \(w(r)\)；
  3. 把所有缩放倍数的半径和加起来：\(\displaystyle{\sum_{t=1}^{\lfloor n/r\rfloor} (tr)=r\cdot \frac{M(M+1)}2,M=\Bigl\lfloor\frac nr\Bigr\rfloor.}\)
  4. 总和：\(\displaystyle{S(n)=\sum_{r\le n} w(r)\cdot r\cdot \frac{M(M+1)}2.}\)

设基本半径\(\displaystyle{r=\prod_{i=1}^{\omega} p_i^{e_i}, p_i\equiv 1\pmod 4,}\)其中 \(\omega\) 为不同素因子个数。

在高斯整数环 \(\mathbb Z[i]\) 中，每个 \(p_i\) 都分裂为共轭的一对高斯素数\(p_i=\pi_i\overline{\pi_i}.\) 因此 \(p_i^{e_i}\) 的所有高斯因子（忽略单位元 \({\pm1,\pm i}\)）可以写成\(\pi_i^{j}\overline{\pi_i}^{e_i-j}, j=0,1,\dots,e_i.\)

把 \(\pi_i\) 的辐角记为 \(\alpha_i=\arg(\pi_i)\in(0,\pi/2)\)，则上述因子的辐角（模 \(\pi/2\)）等价于 \((2j-e_i)\alpha_i(\bmod\pi/2).\)除去落在坐标轴方向（即辐角为 \(0\) 或 \(\pi/2\) 的倍数）的“边界情形”后，这一族给出大约 \(2e_i\) 个互不等价的“非退化辐角选择”；当 \(e_i\) 为偶数时存在 \(j=e_i/2\) 对应的对称中心（\(2j-e_i=0\)），而当 \(e_i\) 为奇数时不存在这一中心，从而会在后面的整体计数里引入一个额外的奇偶差异。

将 \(\omega\) 个素因子的贡献视作相互独立并进行组合，同时再把下列显然的等价关系商掉：

  • 乘以单位元与整体共轭带来的等价（对应刚体旋转/镜像，不改变边长集合）；
  • 三角形三顶点置换带来的对称（对应边长的排列）。

最终得到一个只依赖指数模式 \((e_1,\dots,e_\omega)\) 的闭式权重：

\[ w(r)=2\left(\prod_{i=1}^{\omega} e_i\right)6^{\omega-1}-\delta\cdot 2^{\omega-1}, \]

其中\(\delta=\mathbb{1}\{\gcd(e_1,\dots,e_\omega)\not\equiv 0\pmod2\}\)是示性函数。

这个式子的含义可以分成“主项 + 修正项”两部分理解：

  • 主项 \(2\left(\prod e_i\right)6^{\omega-1}\)：把每个素因子 \(p_i^{e_i}\) 的非退化辐角选择视作 \(e_i\) 个“正代表”（与其共轭配对后由外面的系数 \(2\) 统一处理），而每引入一个新的独立素因子，其辐角贡献在三角形的三条边（或三个角）之间有 \(3\) 个位置可分配，并且可伴随一次共轭方向的二元选择，合并为 \(6\) 种本质不同的组合方式，因此得到乘子 \(6^{\omega-1}\)。

  • 修正项 \(-\delta\cdot 2^{\omega-1}\)：主项的推导隐含了一个前提——三角形的三个角都严格落在 \((0,\pi)\)，从而构造不会落到边界，并且“每个构造都与其镜像（整体共轭）成对出现”。 当 \(r\) 不是完全平方数（等价于存在奇指数）时，会出现一类额外的“轴向/边界”组合：某个由高斯积生成的关键复数落在坐标轴方向，使得对应角度成为 \(0\) 或 \(\pi\)（模 \(\pi\)），几何上表现为三点共线的退化情形，或在“镜像/共轭配对”中发生粘连，导致同一个无向三角形被主项方案重复计数一次。这类边界异常在各素因子之间是相互独立的：对每个 \(p_i^{e_i}\)，要把它的辐角贡献推到轴向，只需在\(\pi_i^{e_i}\) 与 \(\overline{\pi_i}^{e_i}\)两个极端方向之间二选一，因此异常配置的符号模式数为 \(2^\omega\)。但其中“对所有素因子同时取共轭”（把每个 \(\pi_i^{e_i}\) 换成 \(\overline{\pi_i}^{,e_i}\)）对应整体镜像，在边界情形下并不会产生新的异常类型，因此需要除去这一全局二重性，异常类型数为\(2^{\omega-1}\)。

因此只要：筛出所有 \(p\le n,p\equiv 1\pmod 4\)；用 DFS（按非降素数顺序）枚举所有 \(r\)；在 DFS 过程中维护 \(w(r)\) 所需的状态即可。

代码

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
from tools import *

N = 10_000_000

def S(n: int) -> int:
    """
    Compute S(n): sum of circumradii of all integer-sided triangles
    whose circumradius is integer and <= n.
    """
    primes = [p for p in Factorizer(n).primes if p >= 5 and p % 4 == 1]
    plen = len(primes)

    # stack items:
    # (limit0, prod0, idx0, r0, e0, pow2_0, pow6_0, all_even_prefix)
    #
    # Invariant:
    #   prod0 is current basic radius r
    #   limit0 = n // prod0
    #   e0 is exponent of the last used prime (meaningful when prod0 != 1)
    #   r0 = 2 * product of exponents of all primes BEFORE the last prime
    #   pow2_0 = 2^(omega-1), pow6_0 = 6^(omega-1)
    #   all_even_prefix: whether all finalized exponents (before last prime) are even
    stack = [(n, 1, 0, 0, 0, 0, 0, True)]

    res2 = 0  # store 2*S(n) to avoid fractions

    while stack:
        limit0, prod0, idx0, r0, e0, pow2_0, pow6_0, all_even_prefix = stack.pop()

        for i in range(idx0, plen):
            p = primes[i]
            if p > limit0:
                break

            prod = prod0 * p
            limit = limit0 // p  # = n // prod

            # update the (rcoef, exponent, pow2, pow6, all_even_prefix) for this prod
            if prod0 == 1:
                # first prime starts: omega = 1, exponent = 1
                rcoef = 2
                e = 1
                pow2 = 1
                pow6 = 1
                all_even = True  # no previous exponents
            elif i == idx0:
                # repeat the same prime -> exponent increases
                rcoef = r0
                e = e0 + 1
                pow2 = pow2_0
                pow6 = pow6_0
                all_even = all_even_prefix
            else:
                # switch to a new prime -> finalize previous exponent e0
                rcoef = r0 * e0
                e = 1
                pow2 = pow2_0 * 2
                pow6 = pow6_0 * 6
                all_even = all_even_prefix and (e0 % 2 == 0)

            # weight w(prod) = 2*(prod exponents)*6^(omega-1) - delta*2^(omega-1)
            # Here rcoef * e * pow6 equals 2*(product of all exponents)*6^(omega-1).
            weight = rcoef * e * pow6
            if (e & 1) or (not all_even):
                weight -= pow2

            # contribution: weight * prod * sum_{t=1..limit} t = weight * prod * limit*(limit+1)/2
            res2 += weight * prod * limit * (limit + 1)

            # if we can multiply further by primes >= p, push state
            if limit >= p:
                stack.append((limit, prod, i, rcoef, e, pow2, pow6, all_even))

    return res2 // 2


print(S(N))