Project Euler 799

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

Project Euler 799

题目

Pentagonal Puzzle

Pentagonal numbers are generated by the formula: \(P_n = \dfrac 12n(3n-1)\) giving the sequence:

\[1,5,12,22,35, 51,70,92,\ldots \]

Some pentagonal numbers can be expressed as the sum of two other pentagonal numbers.

For example:

\[P_8 = 92 = 22 + 70 = P_4 + P_7\]

\(3577\) is the smallest pentagonal number that can be expressed as the sum of two pentagonal numbers in two different ways

\[\begin{aligned} P_{49} = 3577 & = 3432 + 145 = P_{48} + P_{10} \\ & = 3290 + 287 = P_{47}+P_{14} \end{aligned}\]

\(107602\) is the smallest pentagonal number that can be expressed as the sum of two pentagonal numbers in three different ways.

Find the smallest pentagonal number that can be expressed as the sum of two pentagonal numbers in over \(100\) different ways.

解决方案

从\(P_n=\dfrac{1}{2}n(3n-1)\)出发，注意到\(24P_n=12n(3n-1)=36n^2-12n,\)于是\(24P_n+1=36n^2-12n+1=(6n-1)^2.\)因此 \(P_x+P_y=P_z\)等价于\(24P_x+24P_y=24P_z,\)再分别加 \(1\) 并整理得\((6x-1)^2+(6y-1)^2=(6z-1)^2+1.\)

令\(A=6x-1, B=6y-1, C=6z-1,\)则\(A^2+B^2=C^2+1.\)同时因为 \(x,y,z\) 为正整数，所以\(A\equiv B\equiv C\equiv 5\pmod 6,\)并且 \(A,B,C\) 均为正奇数。因此，对每个 \(z\)，我们只需研究\(N=C^2+1=(6z-1)^2+1,\)并统计满足\(A^2+B^2=N, A\equiv B\equiv 5\pmod 6, A,B>0\) 的不同无序对 \({A,B}\) 的个数，这个个数就等于 \(P_z\) 作为两个五边形数之和的不同表示方式数。

若奇素数 \(p\equiv 3\pmod 4\) 且 \(p\mid (C^2+1)\)，则\(C^2\equiv -1\pmod p.\)这意味着 \(-1\) 在模 \(p\) 下是二次剩余，但当 \(p\equiv 3\pmod 4\) 时，\(-1\) 不是二次剩余，矛盾。

因此\(C^2+1\)的所有奇素因子都满足\(p\equiv 1\pmod 4,\)而且由于 \(C\) 为奇数，\(C^2+1\) 一定是偶数，且\(C^2\equiv 1\pmod 8\Rightarrow C^2+1\equiv 2\pmod 8,\)所以\(C^2+1\)的有且仅有一个质因子\(2\)。

在高斯整数环 \(\mathbb Z[i]\) 中：若素数 \(p\equiv 1\pmod 4\)，则存在整数 \(a,b\) 使\(p=a^2+b^2,\)并且\(p=(a+bi)(a-bi)\).另外\(2\) 满足\(2=(1+i)(1-i).\)

令\(\displaystyle{N=\prod_{j=1}^{t} p_j^{e_j}}\)为素因子分解（其中每个奇素数 \(p_j\equiv 1\pmod 4\)，且 \(2\) 的指数为 \(1\)）。对每个分裂素数 \(p=a^2+b^2\)，记高斯素数\(\pi=a+bi, \overline\pi=a-bi.\)那么所有满足范数为 \(p^e\) 的高斯整数（忽略单位元差别）可以由\(\pi^k\overline{\pi}^{,e-k} (k=0,1,\ldots,e)\)生成，一共有 \(e+1\) 种。

由于范数满足乘法性，\(N\) 的所有平方和表示都可以通过逐素因子做笛卡尔积乘法得到。最终得到所有\(u+vi\in\mathbb Z[i], u^2+v^2=N.\)

将每个结果规范化为\((|u|,|v|)\)并令\(|u|\le |v|,\)即可对应到一个唯一无序对 \({u,v}\)。最后加上模 \(6\) 约束：\(u\equiv v\equiv 5\pmod 6,\)统计满足者的个数，得到该 \(z\) 的表示方式数。

若要表示方式数超过 \(M=100\)，则 \(N=C^2+1\) 必须有相当丰富的平方和分解，这通常意味着它含有很多 \(4k+1\) 型素因子与指数。

对每个素数 \(p\equiv 1\pmod 4\)，条件\(p\mid(C^2+1)\)等价于存在 \(i\) 满足\(i^2\equiv -1\pmod p, C\equiv \pm i\pmod p.\) 而 \(C=6z-1\)，因此\(6z\equiv 1\pm i\pmod p,\)从而得到 \(z\) 在模 \(p\) 下的两个同余类。用这些同余类做分块筛选，可以快速定位“被很多小素数整除”的候选 \(z\)，再对候选做完整分解与计数。

代码

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
from functools import cache

from tools import *
M = 100

def P(n: int) -> int:
    return n * (3 * n - 1) // 2

@cache
def sqrt_minus1_modp(p: int) -> int:
    return int(sympy.sqrt_mod(-1, p))


@cache
def decompose_prime(p: int) -> Tuple[int, int]:
    if p == 2:
        return (1, 1)
    t = sqrt_minus1_modp(p)
    r0, r1 = p, t
    while r1 * r1 > p:
        r0, r1 = r1, r0 % r1
    a = r1
    b2 = p - a * a
    b = int_sqrt(b2)
    return a, b


def gmul(x: Tuple[int, int], y: Tuple[int, int]) -> Tuple[int, int]:
    return x[0] * y[0] - x[1] * y[1], x[0] * y[1] + x[1] * y[0]


def gconj(x: Tuple[int, int]) -> Tuple[int, int]:
    return x[0], -x[1]


def count_ways(z: int) -> int:
    c = 6 * z - 1
    N = c * c + 1
    fac = factorization(N)

    reps = {(1, 0)}
    for p, e in fac:
        a, b = decompose_prime(p)
        pi = (a, b)
        pows = [(1, 0)]
        for _ in range(e):
            pows.append(gmul(pows[-1], pi))

        fs = []
        for k in range(e + 1):
            fs.append(gmul(pows[k], gconj(pows[e - k])))

        nxt = set()
        for r in reps:
            for f in fs:
                nxt.add(gmul(r, f))
        reps = nxt

    ok = set()
    for u, v in reps:
        u = abs(u)
        v = abs(v)
        if u > v:
            u, v = v, u
        if u % 6 == 5 and v % 6 == 5:
            ok.add((u, v))
    return len(ok)


def solve(M):
    prime_limit = 400
    threshold = 7
    block = 1 << 20

    ps = Factorizer(prime_limit).primes
    ps = [p for p in ps if p != 2 and p % 4 == 1]

    roots = []
    for p in ps:
        inv6 = pow(6, p - 2, p)
        i = sqrt_minus1_modp(p)
        r1 = ((1 + i) * inv6) % p
        r2 = ((1 - i) * inv6) % p
        roots.append((p, r1, r2))

    base = 1
    while True:
        cnt = bytearray(block)
        for p, r1, r2 in roots:
            s = (r1 - base) % p
            for i in range(s, block, p):
                cnt[i] += 1
            s = (r2 - base) % p
            for i in range(s, block, p):
                cnt[i] += 1

        for i, v in enumerate(cnt):
            if v < threshold:
                continue
            z = base + i
            w = count_ways(z)
            if w > M:
                return P(z)
        base += block


print(solve(M))