Project Euler 827

发表于 2026-01-21 分类于 Project Euler 阅读次数：本文字数： 3.5k 阅读时长 ≈ 3 分钟

Project Euler 827

题目

Pythagorean Triple Occurrence

Define \(Q(n)\) to be the smallest number that occurs in exactly \(n\) Pythagorean triples \((a,b,c)\) where \(a \lt b \lt c\).

For example, \(15\) is the smallest number occurring in exactly 5 Pythagorean triples: \[(9,12,\mathbf{15})\quad(8,\mathbf{15},17)\quad(\mathbf{15},20,25)\quad(\mathbf{15},36,39)\quad(\mathbf{15},112,113)\]

and so \(Q(5) = 15\)

You are also given \(Q(10)=48\) and \(Q(10^3)=8064000\).

Find \(\displaystyle \sum_{k=1}^{18} Q(10^k)\). Give your answer modulo \(409120391\).

解决方案

令 \(m\) 为正整数，考虑它在所有满足 \(a<b<c\) 的勾股三元组中出现的次数。它可能作为：斜边 \(c=m\)，或者是直角边之一（即 \(a=m\) 或 \(b=m\)）。

为了得到一个可计算的公式，将 \(m\) 按模 \(4\) 的素因子拆开：\(\displaystyle{m = 2^t \prod_i p_i^{e_i} \prod_j q_j^{f_j},}\)其中 \(p_i\equiv 1\pmod 4,q_j\equiv 3\pmod 4,t\ge 0\)。

记\(\displaystyle{E=\prod_i (2e_i+1), F=\prod_j (2f_j+1).}\)下面分别计算两类出现次数。

作为斜边：\(a^2+b^2=m^2\)

等价于求正整数解 \((a,b)\) 使得\(a^2+b^2=m^2.\)对任意正整数 \(N\)，记 \(r_2(N)\) 为方程\(x^2+y^2=N\)的整数解 \((x,y)\) 的个数（有序，且包含符号）。

并且 \(r_2\) 是积性函数，且对素数幂有经典公式：

若 \(p\equiv 1\pmod 4\)，则\(r_2(p^k)=4(k+1).\)
若 \(q\equiv 3\pmod 4\)，则\(r_2(q^k)=\begin{cases}4,& k\equiv 0\pmod 2,\\0,& k\equiv 0\pmod 1,\end{cases}\)
对 \(2\) 的幂：\(r_2(2^k)=4(k\ge 0).\)

对 \(m^2\) 而言，有\(\displaystyle{m^2 = 2^{2t}\prod_i p_i^{2e_i}\prod_j q_j^{2f_j},}\) 因此 \[ r_2(m^2)=4\prod_i(2e_i+1)=4E. \]

\(r_2(m^2)\) 统计的是有序整数组对 \((x,y)\)（包括符号与顺序）。将其转换为正整数且 \(0<a<b\) 的三元组数量：

去掉符号：除以 \(4\)；
去掉顺序：\((a,b)\) 与 \((b,a)\) 视为同一条边对，除以 \(2\)；
排除退化表示 \((m,0)\) 与 \((0,m)\)：这对应 \(4\) 个有序解，因此在除法前先减去 \(4\)。

于是\(m\) 作为斜边出现的次数为

\[ H(m)=\frac{r_2(m^2)-4}{8}=\frac{4E-4}{8}=\frac{E-1}{2}. \]

作为直角边：\(m^2+x^2=y^2\)

考虑方程\(m^2+x^2=y^2, x>0,y>0.\)整理为\(y^2-x^2=m^2 \iff (y-x)(y+x)=m^2.\)

令\(u=y-x, v=y+x,\)则\(uv=m^2, 0<u<v, u\equiv v\pmod 2.\)

因此\(m\) 作为直角边出现的次数等于满足 \(uv=m^2\) 且 \(u<v\) 同奇偶的因子对数量。

若 \(m\) 为奇数（即 \(t=0\)），则 \(m^2\) 为奇数，\(u,v\) 必为奇数，奇偶限制自动满足。此时因子对数为\(\dfrac{d(m^2)-1}{2},\)其中 \(d(\cdot)\) 为约数个数。由于\(\displaystyle{d(m^2)=\prod_i(2e_i+1)\prod_j(2f_j+1)=EF,}\)得

\[C(m)=\frac{EF-1}{2}.\]

若 \(m\) 为偶数（即 \(t\ge 1\)），要使 \(u,v\) 同奇偶且满足 \(y=\dfrac{u+v}{2}\) 为整数，必须 \(u,v\) 同为偶数。写成\(u=2u', v=2v',\) 则\(u'v'=\dfrac{m^2}{4}.\)因而因子对数为\(\dfrac{d(m^2/4)-1}{2}.\)

注意\(\displaystyle{\frac{m^2}{4}=2^{2t-2}\prod_i p_i^{2e_i}\prod_j q_j^{2f_j},}\)所以\(d(m^2/4)=(2t-1)EF,\)从而

\[ C(m)=\frac{(2t-1)EF-1}{2}. \]

因此，总次数为\(R(m)=H(m)+C(m).\)

若 \(t=0\)（\(m\) 为奇数）：\(R(m)=\dfrac{E-1}{2}+\dfrac{EF-1}{2}=\dfrac12E(1+F)-1.\)
若 \(t\ge 1\)（\(m\) 为偶数）：\(R(m)=\dfrac{E-1}{2}+\dfrac{(2t-1)EF-1}{2}=\dfrac12E\bigl(1+(2t-1)F\bigr)-1.\)

为了统一处理，令\(g(t)=\begin{cases}1,& t=0,\\2t-1,& t\ge 1,\end{cases}\)则总公式写为

\[ R(m)=\frac12E(1+g(t)F)-1, \]

并且 \(g(t)\) 永远为正奇数。

要求 \(R(m)=n\)，即\(\dfrac12E(1+gF)-1=n\iff E(1+gF)=2(n+1).\)

令\(S=n+1, T=2S,\)则约束为\(E(1+gF)=T.\)

其中，\(E\) 是若干奇数的乘积（每个因子形如 \(2e+1\)），因此 \(E\) 为奇数；\(g\) 为正奇数；\(F\) 为奇数。因此 \(1+gF\) 必为偶数。

由\(m\)的分解可知，当 \(E,F,t\) 固定后，最小化 \(m\) 的策略是：先确定指数多重集 \(\{e_i\}\) 与 \(\{f_j\}\)；再将最大的指数分配给最小的素数。

约束 \(E=\prod(2e_i+1)\)（\(F=\prod(2f_i+1)\)也同理）等价于将 \(E\) 做一个奇数乘法划分：\(E = r_1r_2\cdots r_k, r_i\ge 3,r_i\equiv 1\pmod2,\)并令\(e_i=\frac{r_i-1}{2}.\)

每种划分对应一个指数序列，再根据指数降序、素数升序的规则即可得到最小的 \(\prod p_i^{e_i}\)（同理得到最小的 \(\prod q_j^{f_j}\)）。

因此，对每个候选 \(E\)（或 \(F\)），可以用带记忆化的递归在其约数结构上找到最优划分。

最终，我们从\(E(1+gF)=T\)出发：

枚举所有奇数因子 \(E\mid S\)（等价于 \(E\mid T\) 且 \(E\) 为奇数）；
定义\(B=\dfrac{T}{E}, D=B-1\)，那么\(B\)为偶数，\(D\)为奇数；
则需要\(D=gF,\)因而枚举所有奇数因子 \(g\mid D\)，并令\(F=\dfrac{D}{g}.\)
由 \(g\) 反推出 \(2\) 的指数 \(t\)：若 \(g=1\)，取 \(t=0\)（这一定优于 \(t=1\)）；若 \(g>1\)，则\(g=2t-1\iff t=\dfrac{g+1}{2}.\)

对每个枚举出来的三元组 \((E,g,F)\)，构造的最小数为

\[ m_{\min}(E,g,F)=2^t\cdot \min\left\{\prod p_i^{e_i}\right\}\cdot \min\left\{\prod q_j^{f_j}\right\}, \]

也就是说，两侧的最小化分别在所有满足 \(\prod(2e_i+1)=E\) 与 \(\prod(2f_j+1)=F\) 的指数方案中进行。

最终\(Q(n)\)则是上面的枚举\(m_{\min}(E,g,F)\)结果中的最小值。

代码

from functools import lru_cache
from math import log
from tools import divisors, get_prime

MOD = 409120391
M = 18




def primes_mod4(cnt: int, m: int):
    return [p for p in get_prime(10 ** 5) if p % 4 == m][:cnt]


P1 = primes_mod4(100, 1)
P3 = primes_mod4(100, 3)
L1 = [log(p) for p in P1]
L3 = [log(p) for p in P3]
L2 = log(2)


def build_minimizer(P, L):
    @lru_cache(None)
    def odd_divs_ge3(x: int):
        ds = divisors(x)
        return tuple(d for d in ds if (d & 1) and d >= 3)

    @lru_cache(None)
    def best(prod: int, lim: int, idx: int):
        if prod == 1:
            return 0.0, ()
        best_log = float('inf')
        best_fac = ()
        ds = odd_divs_ge3(prod)
        for d in reversed(ds):
            if d > lim:
                continue
            e = (d - 1) >> 1
            sub_log, sub_fac = best(prod // d, d, idx + 1)
            val = e * L[idx] + sub_log
            if val < best_log:
                best_log = val
                best_fac = (d,) + sub_fac
        return best_log, best_fac

    @lru_cache(None)
    def solve(prod: int):
        if prod == 1:
            return 0.0, 1
        lg, facs = best(prod, prod, 0)
        v = 1
        for i, d in enumerate(facs):
            e = (d - 1) >> 1
            v = (v * pow(P[i], e, MOD)) % MOD
        return lg, v

    return solve


min_part_1 = build_minimizer(P1, L1)
min_part_3 = build_minimizer(P3, L3)


def Q_value(n: int) -> int:
    S = n + 1
    T = 2 * S

    best_log = float('inf')
    best_val = 0

    for E in divisors(S):
        if (E & 1) == 0:
            continue
        B = T // E
        D = B - 1

        logE, modE = min_part_1(E)

        for g in divisors(D):
            F = D // g
            if g == 1:
                t = 0
            else:
                t = (g + 1) >> 1

            logF, modF = min_part_3(F)

            cur_log = logE + logF + t * L2
            if cur_log < best_log:
                best_log = cur_log
                best_val = pow(2, t, MOD) * modE * modF % MOD

    return best_val


ans = 0
p10 = 1
for _ in range(M):
    p10 *= 10
    ans = (ans + Q_value(p10)) % MOD
print(ans)