Project Euler 764

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

Project Euler 764

题目

Asymmetric Diophantine Equation

Consider the following Diophantine equation: \[16x^2+y^4=z^2\] where \(x\), \(y\) and \(z\) are positive integers.

Let \(S(N) = \displaystyle{\sum(x+y+z)}\) where the sum is over all solutions \((x,y,z)\) such that \(1 \leq x,y,z \leq N\) and \(\gcd(x,y,z)=1\).

For \(N=100\), there are only two such solutions: \((3,4,20)\) and \((10,3,41)\). So \(S(10^2)=81\).

You are also given that \(S(10^4)=112851\) (with \(26\) solutions), and \(S(10^7)\equiv 248876211 \pmod{10^9}\).

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

解决方案

令\(N=10^{16}\)。我们将方程改写为\((4x)^2+(y^2)^2=z^2.\)令\(A=4x, B=y^2, C=z,\)则 \((A,B,C)\) 是一组勾股三元组：\(A^2+B^2=C^2.\)

勾股三元组的标准参数化（Euclid 公式）为：存在正整数 \(k,m,n\)（\(m>n\)、\(\gcd(m,n)=1\)、\(m\not\equiv n\pmod 2\)），使得

\[(A,B,C)=\bigl(k\cdot 2mn,k\cdot (m^2-n^2),k\cdot (m^2+n^2)\bigr),\]

或交换两条直角边：

\[(A,B,C)=\bigl(k\cdot (m^2-n^2),k\cdot 2mn,k\cdot (m^2+n^2)\bigr).\]

因此，本题分成两个放置方式的分支讨论。

第一类：令 \(4x=2kmn\), \(y^2=k(m^2-n^2)\)

此时\(x=\dfrac{k mn}{2}, y^2=k(m^2-n^2), z=k(m^2+n^2).\)

若某素数 \(p\mid k\)，则 \(p\mid y^2\Rightarrow p\mid y\)，同时 \(p\mid z\)，且 \(p\mid x\)（因为 \(x\) 也含因子 \(k\)）。这会导致 \(\gcd(x,y,z)\ge p\)，与题设互素矛盾。故必须\(k=1.\)

于是得到\(4x=2mn, y^2=m^2-n^2, z=m^2+n^2.\)

接下来，我们强制 \(y^2=m^2-n^2\) 为完全平方数。我们要求存在整数 \(y\) 使\(m^2-n^2=y^2.\)改写为\((m-n)(m+n)=y^2.\)设\(m-n=a^2, m+n=b^2,\) 则 \(y=ab\)。为了使 \(m,n\) 为整数，必须 \(a^2,b^2\) 同奇偶。又因为 \(m,n\) 在原始勾股参数中异奇偶，所以这里自然落到 \(a,b\) 同为奇数 的情形。

由上式解出\(m=\dfrac{a^2+b^2}{2}, n=\dfrac{b^2-a^2}{2}.\)并且若 \(\gcd(a,b)=1\)，则\(\gcd(m,n)=\dfrac{\gcd(m+n,m-n)}{2}=\dfrac{\gcd(b^2,a^2)}{2}=1,\)满足原始参数条件。

计算 \(x,z\)得到：\(x=\dfrac{mn}{2}=\dfrac{1}{2}\cdot \dfrac{a^2+b^2}{2}\cdot \dfrac{b^2-a^2}{2}=\dfrac{b^4-a^4}{8},z=m^2+n^2=\dfrac{(a^2+b^2)^2+(b^2-a^2)^2}{4}=\dfrac{a^4+b^4}{2}.\)

于是第一类所有解可写为

\[(x,y,z)=\left(\frac{b^4-a^4}{8},ab,\frac{a^4+b^4}{2}\right),\]

其中\(a,b\) 为正奇数，\(\gcd(a,b)=1\)且 \(b>a\)（保证 \(x>0\)）。

这一类解的 \(y=ab\) 为奇数，因此它们全部对应 奇数 \(y\)。

第二类：令 \(4x=k(m^2-n^2)\), \(y^2=2kmn\)

此时\(4x=k(m^2-n^2), y^2=2kmn, z=k(m^2+n^2),\)其中 \(\gcd(m,n)=1\)、\(m\not\equiv n\pmod 2\)，故 \(m^2-n^2\) 为奇数。

因为左边 \(4x\) 必被 \(4\) 整除，而右边 \(k(m^2-n^2)\) 中 \((m^2-n^2)\) 为奇数，所以必须 \(4\mid k\)。

写 \(k=4k'\)。则\(x=k'(m^2-n^2), y^2=8k'mn, z=4k'(m^2+n^2).\)若 \(k'>1\)，则 \(k'\mid x\)、且 \(k'\mid z\)，并且从 \(y^2\) 得到 \(k'\mid y\)，因此 \(\gcd(x,y,z)\ge k'>1\)，这是不允许的。所以\(k'=1\Rightarrow k=4.\)

于是得到\(x=m^2-n^2, y^2=8mn, z=4(m^2+n^2).\)

我们接下来强制 \(8mn\) 为完全平方。由 \(\gcd(m,n)=1\) 且一奇一偶，不失一般性，可设 \(m\) 为偶、\(n\) 为奇。由于互素，\(mn\) 的平方因子只能分别来自两边。

若 \(m\) 为偶，则把 \(m=2u\)，则\(y^2=8mn=16un.\)要成为平方，必须 \(u\) 与 \(n\) 都是完全平方（因为 \(\gcd(u,n)=1\) 且 \(n\) 为奇数）。所以设\(u=a^2, n=b^2.\)其中\(b\)为奇数。

于是 \(m=2a^2\)，并且\(y=4ab.\)计算 \(x,z\)：\(x=m^2-n^2=(2a^2)^2-b^4=4a^4-b^4,z=4(m^2+n^2)=4\bigl((2a^2)^2+b^4\bigr)=4(4a^4+b^4).\)

为了覆盖 \(x>0\) 的所有情况，我们取绝对值即可（本题 \(x\) 定义为正整数）：\(x=|4a^4-b^4|.\)因此第二类原始解为

\[(x,y,z)=\left(|4a^4-b^4|,4ab,4(4a^4+b^4)\right),\]

其中\(b\) 必为奇数，且\(\gcd(a,b)=1\)。

这一类解的 \(y\) 必为偶数（甚至是 \(4\) 的倍数），因此与第一类完全不重叠。

两类解都满足 \(z\ge 4x\)，所以只要 \(z\le N\)，通常就自动有 \(x\le N\)；并且在 \(N\) 的量级下 \(y\) 远小于 \(N\)，最终约束可直接用 \(z\le N\) 来控制枚举范围。

但是即使严格按定义 \(x,y,z\le N\)，下面的推导与实现仍然正确，因为我们直接用 \(z\) 给出的必要条件来枚举（它是最紧的主约束）。

两类解都需要满足参数互素（第一类要求 \(\gcd(a,b)=1\)；第二类同样要求 \(\gcd(a,b)=1\)，并额外限制 \(b\) 为奇数）。为了把互素这一条件从枚举中剥离出来，可以使用 Möbius 反演。

具体地，设 \(F(u,v)\) 是定义在正整数对上的权重函数，并且它只在一个有限范围内取非零值（例如带有 \(z(u,v)\le N\) 之类的截断条件，使得总和实际是有限的）。那么有恒等式：

\[ \sum_{\gcd(u,v)=1} F(u,v)=\sum_{g\ge 1}\mu(g)\sum_{u,v}F(gu,gv). \]

在本题中，\(F(u,v)\) 就可以理解为由参数对 \((u,v)\) 生成的解对总和的贡献，即\(F(u,v)=(x(u,v)+y(u,v)+z(u,v)),\)。

关键点在于：两类构造中 \(x\) 与 \(z\) 对参数是四次齐次，而 \(y\) 是二次齐次，因此缩放 \((u,v)\mapsto(gu,gv)\) 时有\(x,z\) 乘 \(g^4\),\(y\) 乘 \(g^2\)。

于是对固定 \(g\)，只需把上界 \(N\) 替换为 \(\lfloor N/g^4\rfloor\)，就能把“互素参数对的求和”转化为“无互素约束的区间求和”，从而实现高效计算。

核心技巧是：对固定的 \(b\)，允许的 \(a\) 总是形如\(1\le a\le A_{\max}(b)\)（并带有奇偶步长）也就是说，只要固定了\(b\)，\(a\)就会有一个由\(b\)确定的上界\(A_{\max}(b)\)，因此所有表达式对 \(a\) 的求和都能用前缀和 \(\sum a\)、\(\sum a^4\) 在 \(O(1)\) 时间完成。

我们定义：\(\displaystyle{P_4(t)=\sum_{i=1}^t i^4=\frac{t(t+1)(2t+1)(3t^2+3t-1)}{30},Q_4(t)=\sum_{i=1}^t (2i-1)^4=\frac{t(2t-1)(2t+1)(12t^2-7)}{15}.}\)

第一类解对固定 \(b\) 的求和

第一类的所有解满足形式：\(x=\dfrac{b^4-a^4}{8},\ y=ab,\ z=\dfrac{a^4+b^4}{2},\)其中 \(a,b\) 为奇数，且 \(a<b\)，并满足 \(\dfrac{a^4+b^4}{2}\le \dfrac{N}{g^4}\)。

对固定 \((g,b)\)，由上界约束可得 \(a\) 的最大可能值为\(a\le \sqrt[4]{\dfrac{2N}{g^4}-b^4}.\)令满足上述约束的奇数 \(a\) 的个数为 \(t=t(g,b)\)，即 \(a=1,3,\dots,2t(g,b)-1\)，则 \(\displaystyle{\sum a=t(g,b)^2,\ \sum a^4=Q_4(t(g,b)).}\)

因此在缩放 \((a,b)\mapsto(ga,gb)\) 后，有 \[\displaystyle{ \sum x' = g^4\cdot \frac{t(g,b)b^4-Q_4(t(g,b))}{8}, \sum y' = g^2\cdot bt(g,b)^2, \sum z' = g^4\cdot \frac{Q_4(t(g,b))+t(g,b),b^4}{2}, }\]

随后对 \(\displaystyle{\sum x'+y'+z'}\) 乘上 \(\mu(g)\) 并求和叠加即可。

第二类解对固定 \(b\) 的求和

第二类解满足形式：\(x=|4a^4-b^4|,\ y=4ab,\ z=4(4a^4+b^4),\)其中 \(b\) 必为奇数，且满足 \(4(4a^4+b^4)\le \dfrac{N}{g^4}\)。

对固定 \((g,b)\)（\(b\) 为奇数），由上界约束可得 \(a\) 的最大可能值为\(a\le A_{\max}(g,b)=\left\lfloor \sqrt[4]{\dfrac{N}{16g^4}-\dfrac{b^4}{4}} \right\rfloor.\)

因此允许的 \(a\) 构成连续区间 \(1\le a\le A_{\max}(g,b)\)。为处理绝对值项，将求和按分界点\(4a^4\le b^4 \iff a\le T(b)=\left\lfloor\sqrt[4]{\dfrac{b^4}{4}}\right\rfloor\)拆分为两段，其中实际分界取 \(T=\min\bigl(T(b),A_{\max}(g,b)\bigr)\)。于是\(\displaystyle{\sum_{a=1}^{A_{\max}} |4a^4-b^4|=\sum_{a=1}^{T}(b^4-4a^4)+\sum_{a=T+1}^{A_{\max}}(4a^4-b^4),}\)从而可用四次方前缀和 \(P_4\) 在 \(O(1)\) 时间内计算该项。

与此同时，\(y,z\) 也可直接用前缀和表示：\(\displaystyle{\sum_{a=1}^{A_{\max}} y=4b\sum_{a=1}^{A_{\max}} a,\sum_{a=1}^{A_{\max}} z=4\sum_{a=1}^{A_{\max}}(4a^4+b^4).}\)

因此在缩放 \((a,b)\mapsto(ga,gb)\) 后，有

\[ \sum x' = g^4\sum_{a=1}^{A_{\max}(g,b)} |4a^4-b^4|, \sum y' = g^2\cdot 4b\sum_{a=1}^{A_{\max}(g,b)} a, \sum z' = g^4\cdot 4\sum_{a=1}^{A_{\max}(g,b)}(4a^4+b^4), \]

随后对所有 \(g\) 乘上 \(\mu(g)\) 并求和叠加即可。

代码

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
from tools import *
N = 10 ** 16

def iroot4(x: int) -> int:
    return int_sqrt(int_sqrt(x))


def mobius_sieve(n: int):
    mu = [0] * (n + 1)
    mu[1] = 1
    primes = []
    is_comp = [0] * (n + 1)
    for i in range(2, n + 1):
        if not is_comp[i]:
            primes.append(i)
            mu[i] = -1
        for p in primes:
            v = i * p
            if v > n:
                break
            is_comp[v] = 1
            if i % p == 0:
                mu[v] = 0
                break
            mu[v] = -mu[i]
    return mu


def solve(N: int) -> int:
    MOD = 10 ** 9
    L = iroot4(2 * N) + 2
    mu = mobius_sieve(L)

    sum4 = [0] * (L + 1)
    for i in range(1, L + 1):
        sum4[i] = sum4[i - 1] + i ** 4

    sumOdd4 = [0] * (L + 1)
    for i in range(1, L + 1):
        sumOdd4[i] = sumOdd4[i - 1] + (2 * i - 1) ** 4

    ans = 0

    for g in range(1, L + 1, 2):
        mg = mu[g]
        if mg == 0:
            continue

        g2 = g * g
        g4 = g2 * g2

        bound2 = (2 * N) // g4
        if bound2 >= 2:
            bmax = iroot4(bound2)
            b = 1
            while b <= bmax:
                b4 = b ** 4
                rem = bound2 - b4
                if rem <= 0:
                    break
                amax = iroot4(rem)
                if amax >= b:
                    amax = b - 2
                if amax >= 1:
                    t = (amax + 1) // 2
                    if t > 0:
                        sa4 = sumOdd4[t]
                        sx = (t * b4 - sa4) // 8
                        sz = (sa4 + t * b4) // 2
                        sy = b * t * t
                        term = g4 * (sx + sz) + g2 * sy
                        ans += mg * term
                b += 2

        bound = N // (4 * g4)
        if bound >= 1:
            bmax = iroot4(bound)
            b = 1
            while b <= bmax:
                b4 = b ** 4
                rem = bound - b4
                if rem <= 0:
                    break
                U = iroot4((rem // 4))
                if U > 0:
                    t1 = iroot4(b4 // 4)
                    if t1 > U:
                        t1 = U
                    S4U = sum4[U]
                    S4t = sum4[t1]

                    sx = (2 * t1 - U) * b4 + 4 * S4U - 8 * S4t
                    sy = 4 * b * (U * (U + 1) // 2)
                    sz = 4 * (4 * S4U + U * b4)

                    term = g4 * (sx + sz) + g2 * sy
                    ans += mg * term
                b += 2

    return ans % MOD

print(solve(N))