Project Euler 360

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

Project Euler 360

题目

Scary Sphere

Given two points \((x_1,y_1,z_1)\) and \((x_2,y_2,z_2)\) in three dimensional space, the Manhattan distance between those points is defined as \(|x_1-x_2|+|y_1-y_2|+|z_1-z_2|\).

Let \(C(r)\) be a sphere with radius \(r\) and center in the origin \(O(0,0,0)\).

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

Let \(S(r)\) be the sum of the Manhattan distances of all elements of \(I(r)\) to the origin \(O\).

E.g. \(S(45)=34518\).

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

解决方案

本题半径\(r=10^{10}=2^{10}\cdot 5^{10}.\)

关键观察：若 \(x^2+y^2+z^2\equiv 0\pmod 4\)，则 \(x,y,z\) 必须全为偶数。因为平方模 \(4\) 只能是 \(0\) 或 \(1\)，三个平方和要是 \(0\)，就不可能出现任何 \(1\)。

而 \(r^2\) 含有因子 \(2^{20}\)，所以\(x^2+y^2+z^2=r^2\equiv 0\pmod{4}\Rightarrow x\equiv y\equiv z\equiv 0\pmod2.\)

把 \((x,y,z)\) 同时除以 \(2\)，右侧 \(r^2\) 也除以 \(4\)，同样仍然是“和为 \(0\pmod4\)”的情形，于是可以重复 \(10\) 次，得到：\(2^{10}\mid x, 2^{10}\mid y, 2^{10}\mid z.\)

令\(x=2^{10}X, y=2^{10}Y, z=2^{10}Z,R=5^{10}.\)代回：\((2^{10})^2(X^2+Y^2+Z^2)=(2^{10}R)^2\Rightarrow X^2+Y^2+Z^2=R^2.\)

所以点集一一对应，且曼哈顿距离线性缩放：\(|x|+|y|+|z|=2^{10}(|X|+|Y|+|Z|).\)因此\(S(10^{10})=2^{10}S(R), R=5^{10}.\)

集合 \(I(R)\) 在坐标置换下不变（因为方程 \(X^2+Y^2+Z^2=R^2\) 对 \(X,Y,Z\) 完全对称），因此\(\displaystyle{\sum_{(X,Y,Z)\in I(R)}|X|=\sum_{(X,Y,Z)\in I(R)}|Y|=\sum_{(X,Y,Z)\in I(R)}|Z|.}\)于是

\[ S(R)=\sum(|X|+|Y|+|Z|) =3\sum_{(X,Y,Z)\in I(R)}|Z|. \]

接下来只需要统计每个 \(Z\) 值出现多少点。固定 \(Z=z\)，方程变为\(X^2+Y^2=R^2-z^2=(R-z)(R+z).\)设\(m(z)=R^2-z^2=(R-z)(R+z).\)令\(r_2(n)=|\{(a,b)\in\mathbb Z^2: a^2+b^2=n\}|\) 为“二平方表示数”（计入符号与顺序）。那么对固定 \(z\neq 0\)，点 \((X,Y,Z)\) 的个数是 \(r_2(m(z))\)，而 \(Z=\pm z\) 两层都存在，因此：

\[ \sum_{(X,Y,Z)\in I(R)}|Z|=2R+\sum_{z=1}^{R-1} 2z\cdot r_2(m(z)). \]

两平方和定理：若\(k\)的质因子分解满足：\(\displaystyle{k=2^e\prod_{p\equiv1\pmod4}p^{\alpha_p}\prod_{q\equiv3\pmod4}q^{2\beta_q}}\)，也就是说，对于所有\(q\equiv 3 \pmod4\)的质因子\(q\)，其对应的指数都满足\(\beta_q\)整数（对于其它质因子则不做要求），则存在平方和表示，并且无符号且有序的表示数为\(S=\{(a,b):a,b\ge 0\}\)：

\[ r_2(k)=\prod_{p\equiv1\pmod4}(\alpha_p+1). \]

不满足两平方定理的\(k\)（或者说有\(q\equiv 3\pmod 4\)的质因子的次数为奇数），那么\(r_2(k)\)=0。

注意\(m(z)=(R-z)(R+z),1\le z\le R-1,R+z\le 2R.\)又因为\(\gcd(R-z,R+z)=\gcd(R-z,2z)=\gcd(2R,R+z).\)这里 \(R\) 是奇数（纯 \(5\) 次幂），所以任何奇素数公因子都必须整除 \(R\)，从而公因子只能是 \(5\) 的幂；尤其不存在任何 \(4k+3\) 型素数会同时整除两者。

因此，判断 \(m(z)\) 是否可表示为二平方和时，所有 \(4k+3\) 型素数的指数奇偶性可以分别在 \((R-z)\) 与 \((R+z)\) 中独立检查：只要两边各自都满足“所有 \(4k+3\) 型素数指数为偶数”，乘积就满足。

接下来只剩一个小细节：对 \(p\equiv 1\pmod4\) 的素数指数，\((\alpha+1)\) 是“乘积上指数求和后再加一”；当 \(p\neq 5\) 时，\(p\) 不会同时出现在 \((R-z)\) 与 \((R+z)\) 的 \(\gcd\) 中，所以 \((\alpha+1)\) 的乘积是完全可乘的；唯一可能同时出现的 \(p\equiv1\pmod4\) 素数是 \(5\)。

设\(A=R-z, B=R+z, e_A=v_5(A), e_B=v_5(B).\)并定义\(\displaystyle{F(n)=\prod_{p\equiv 1 \pmod 4}(\nu_p(n)+1),}\)其中 \(\nu_p(n)\) 是 \(p\) 在 \(n\) 中的指数（注意 \(p=5\) 也算在内）。

那么除 \(5\) 以外的部分完全相乘，而 \(5\) 的因子需要把 \((e_A+1)(e_B+1)\) 替换为 \((e_A+e_B+1)\)，因此： 若\(A,B\)各自满足\(4k+3\)型素数指数全偶，那么\(r_2(AB)=4 F(A)F(B)\cdot \dfrac{e_A+e_B+1}{(e_A+1)(e_B+1)}\)，否则\(r_2(AB)=0\)。

最终，我们通过预处理\(v_5\)和\(F\)函数，通过遍历\(z\)主循环计算出\(r_2\)函数的值。在最后累加计算出如下值：

\[ S(10^{10})=2^{10}\cdot 3\left(2R+\sum_{z=1}^{R-1}2z\cdot r_2((R-z)(R+z))\right). \]

代码

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
from array import array

from tools import Factorizer

W = 10
def precompute_tables(N: int):
    """
    For 1..N compute:
      good[n] : whether all primes == 3 (mod 4) appear to even exponent in n
      F[n]    : product over primes == 1 (mod 4) of (exp+1) in n  (includes p=5)
      v5[n]   : exponent of 5 in n
    """
    spf = Factorizer(N).spf

    # exp[n] = exponent of spf[n] in n
    # base[n] = n // spf[n]^exp[n]
    exp = array("B", [0]) * (N + 1)
    base = array("I", [0]) * (N + 1)

    good = array("B", [0]) * (N + 1)
    F = array("I", [0]) * (N + 1)
    v5 = array("B", [0]) * (N + 1)

    base[1] = 1
    good[1] = 1
    F[1] = 1
    v5[1] = 0

    for n in range(2, N + 1):
        p = spf[n]
        m = n // p
        if spf[m] == p:
            e = exp[m] + 1
            exp[n] = e
            b = base[m]
            base[n] = b
        else:
            e = 1
            exp[n] = 1
            b = m
            base[n] = b

        # good[n]
        g = good[b]
        if g and (p & 3) == 3 and (e & 1):
            g = 0
        good[n] = g

        # F[n]
        f = F[b]
        if (p & 3) == 1:
            f *= (e + 1)
        F[n] = f

        # v5[n]
        if p == 5:
            v5[n] = e
        else:
            v5[n] = v5[b]

    return good, F, v5


def solve(W) -> int:
    r = 10 ** W

    # strip all factors of 2: r = two_pow * R, R odd
    two_pow = 1
    R = r
    while (R & 1) == 0:
        R >>= 1
        two_pow <<= 1

    # for this problem, R should be 5^10
    # we need tables up to 2R
    N = 2 * R
    good, F, v5 = precompute_tables(N)

    # sum_z = sum_{(X,Y,Z) in I(R)} |Z|
    # z=R gives (0,0,±R) => contribution 2R
    sum_z = 2 * R

    # For 1 <= z <= R-1:
    # m(z) = (R-z)(R+z) = A*B
    # if A and B both pass good[], then:
    # r2(A*B) = 4 * F(A)*F(B) * (eA+eB+1)/((eA+1)(eB+1))
    # contribution to sum_z is 2*z*r2 = 8*z*(...)
    for z in range(1, R):
        A = R - z
        B = R + z
        if not good[A] or not good[B]:
            continue

        eA = v5[A]
        eB = v5[B]
        num = (eA + eB + 1)
        den = (eA + 1) * (eB + 1)

        cnt = F[A] * F[B]
        cnt = (cnt * num) // den  # exact integer

        sum_z += 8 * z * cnt

    # S(R) = 3 * sum_z  (symmetry among coordinates)
    # S(r) = two_pow * S(R)  (scaling)
    return two_pow * 3 * sum_z


if __name__ == "__main__":
    print(solve(W))