Project Euler 302

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

Project Euler 302

题目

Strong Achilles Numbers

A positive integer \(n\) is powerful if \(p^2\) is a divisor of \(n\) for every prime factor \(p\) in \(n\).

A positive integer \(n\) is a perfect power if \(n\) can be expressed as a power of another positive integer.

A positive integer \(n\) is an Achilles number if \(n\) is powerful but not a perfect power. For example, \(864\) and \(1800\) are Achilles numbers: \(864 = 2^5\cdot3^3\) and \(1800 = 2^3\cdot3^2\cdot5^2\).

We shall call a positive integer S a Strong Achilles number if both \(S\) and \(\varphi(S)\) are Achilles numbers.\(^1\)

For example, \(864\) is a Strong Achilles number: \(\varphi(864) = 288 = 2^5\cdot3^2\). However, \(1800\) isn’t a Strong Achilles number because: \(\varphi(1800) = 480 = 2^5\cdot3^1\cdot5^1\).

There are \(7\) Strong Achilles numbers below \(10^4\) and \(656\) below \(10^8\).

How many Strong Achilles numbers are there below \(10^{18}\)?

\(^1 \varphi\) denotes Euler’s totient function.

解决方案

令\(M=10^{18}\)。设\(n\)的质因数分解为\(\displaystyle{n=\prod_{i=1}^t p_i^{e_i},p_1\le p_2\le \dots\le p_t}\)。

那么，\(n\) 是 powerful 当且仅当对每个素因子 \(p_i\)，都有 \(p_i^2\mid n\)，也就是 \(\displaystyle{\min_{i=1}^t\{e_i\}\ge 2}\)。

\(n\) 是 perfect power 当且仅当存在整数 \(a\ge 2\) 和 \(m\ge 2\) 使得 \(n=a^m\)。

在素因子指数表示下，\(n\) 是 perfect power 当且仅当\(\gcd(e_1,e_2,\dots,e_t)\ge 2.\)这是因为：\(n=a^m\) 等价于每个素因子指数都能被同一个 \(m\ge 2\) 整除。

因此 \(n\) 是 Achilles 当且仅当：

  • powerful：\(\displaystyle{\min_{i=1}^t\{e_i\}\ge 2}\)
  • 但不是 perfect power：\(\gcd(e_1,\dots,e_t)=1\)

那么\(\varphi(n)\)可以计算得出：\(\varphi(n)=\displaystyle{\prod_{i=1}^t p_i^{e_i-1}\cdot \prod_{i=1}^t (p_i-1).}\)

接下来论证\(e_t\ge 3\)。

可见 \(p_t\) 为 \(n\) 的最大质因子。看 \(\varphi(n)\) 中质数 \(p_t\) 的指数来源：

  • 来自 \(\prod_{i=1}^t p_i^{e_i-1}\)：给出 \(p_t^{e_t-1}\)。
  • 来自 \(\prod_{i=1}^t (p_i-1)\)：可见，对于所有\(p_i\)，都有\(p_i-1<p_t\)，因此不可能含因子 \(p_t\)。

所以\(v_{p_t}(\varphi(n)) = e_{t}-1.\)

若要 \(\varphi(n)\) powerful，则必须 \(v_{p_t}(\varphi(n))\ge 2\)，即可得到\(e_{t}-1\ge 3.\)

这给出非常强的剪枝：\(n\) 的最大质因子必须以至少 3 次方出现。

因此当 \(n<M\) 时，\(p_t^3 \le n \le M \Longrightarrow\ p_t \le \sqrt[3]{M}.\)

所以只需要用到 \(\sqrt[3]{M}\) 以内的质数。

现在我们要数满足 \(S<M\) 且 \(S,\varphi(S)\) 都是 Achilles 的 \(S\)。因此按质数从大到小递归选择是否加入一个因子 \(p^e\)（其中 \(e\ge 2\)），构造\(\displaystyle{n=\prod p^e \le M}\)，并保证递归时后续只允许用更小质数，避免重复。

如果当前选到的指数集合是 \({e_1,\dots,e_t}\)，那么

  • \(n\) 是否 powerful 由构造保证（所有 \(e_i\ge 2\)）
  • \(n\) 是否不是 perfect power 等价于\(g_n=\gcd(e_1,\dots,e_t)=1.\)所以递归里只要维护 \(g_n\) 即可。

由欧拉函数的定义\(\varphi(n)=\displaystyle{\prod_{i=1}^t p_i^{e_i-1}(p_i-1)}\)可知，当我们把 \(p^e\) 乘进 \(n\) 时，对 \(\varphi(n)\) 的影响是：

  • 素数 \(p\) 的指数增加 \(e-1\)
  • 分解 \(p-1\)，把其素因子指数也加到 \(\varphi(n)\) 的指数表中

因此需要维护一个哈希表即可在任意时刻判断 \(\varphi(n)\) 是否 Achilles：

  • powerful：对所有 \(i\)，有 \(v_{p_i}(\varphi(n))\ge 2\)
  • 非 perfect power：\(\gcd(v_{p_1}(\varphi(n)),v_{p_2}(\varphi(n)),\dots,v_{p_t}(\varphi(n)))=1\)

递归过程中，\(\varphi(n)\) 的指数表里可能出现某些质数 \(q\) 满足 \(v_q(\varphi(n))=1\)。如果最终要让 \(\varphi(n)\) powerful，就必须把它提升到至少 \(2\)。

但注意：我们后续递归只会加入更小质数 \(p\)。若 \(p<q\)，则 \(p-1<q\)，所以 \(q\) 不可能出现在 \(p-1\) 的分解中；而 \(p^{e-1}\) 只会影响质数 \(p\) 本身的指数，也不能增加 \(q\) 的指数。

因此：若当前 \(\varphi(n)\) 中存在最大质数 \(q\) 满足 \(v_q(\varphi(n))<2\)，那么后续可选的质数 \(p\) 必须满足 \(p\ge q\)，否则该分支必然无法补齐 \(q\) 的指数。也就是说，我们可以维护一个数据结构来维护最大的单个质因子\(P\)，一旦枚举到的\(p\)小于\(P\)就停止。

代码

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
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
import heapq
from bisect import bisect_left, bisect_right
from math import ceil, isqrt, gcd

from tools import Factorizer
N = 10 ** 18

def solve_strong_achilles(N: int = 10**18) -> int:
    # ------------------------------------------------------------
    # 目标：统计 Strong Achilles 数 S < N
    # 条件：S 是 Achilles（powerful 且不是 perfect power）
    #       phi(S) 也是 Achilles
    #
    # powerful：对每个质因子 p，指数 >= 2
    # not perfect power：所有指数的 gcd == 1
    #
    # 构造思路（常见做法）：
    #   DFS 枚举 n 的质因子分解 n = ∏ p^a，且 a>=2（保证 n powerful）
    #   在 DFS 同时维护 phi(n) 的质因子指数表 phi_exp
    #   到每个 DFS 节点判断：
    #       - n 是否 Achilles：g_n == 1（指数 gcd == 1）
    #       - phi(n) 是否 Achilles：phi_exp powerful 且 gcd(exps)==1 且至少两种质因子
    # ------------------------------------------------------------

    # 由于 n powerful，最小指数为 2，因此若 n < N，则最大质数 p 至少满足 p^2 <= N
    # 更强的剪枝：很多实现从“最大质数指数至少 3”推导出只需 primes <= N^(1/3)，
    # 这里直接采用经典上界：PMAX = ceil(N^(1/3))（对 N=1e18 正好 1e6）
    PMAX = ceil(N ** (1 / 3))
    fact = Factorizer(PMAX)  # 预筛到 PMAX，用于快速分解 <= PMAX 的整数
    primes = fact.primes
    Pn = len(primes)

    # 预处理每个质数 p 的 (p-1) 的质因子分解
    # 因为 phi(∏ p^a) = ∏ p^(a-1) * (p-1)，每次选 p 都要把 (p-1) 的因子加入 phi_exp
    # 预分解能减少 DFS 内重复 factor 的开销
    pm1_fac = [([] if p == 2 else fact.factor_small(p - 1)) for p in primes]

    # ------------------------------------------------------------
    # 维护 phi(n) 的指数表：phi_exp[q] = v_q(phi(n))
    # 注意：phi_exp 只存指数 > 0 的质因子（指数为 0 的不保存）
    # ------------------------------------------------------------
    phi_exp = {}

    # ------------------------------------------------------------
    # bad_cnt：phi(n) 里指数为 1 的质因子个数
    # powerful 要求所有指数 >= 2，因此 powerful <=> bad_cnt == 0
    #
    # 这里之所以能用 “指数==1” 来代表 “指数<2”，是因为：
    # - phi_exp 中指数不会是负数
    # - 我们构造过程中只会增加指数；指数从 0 变 1 的时候就成为 bad
    # - 只要把所有 bad 修补到 >=2，就实现 powerful
    # ------------------------------------------------------------
    bad_cnt = 0

    # ------------------------------------------------------------
    # 剪枝需要“当前 phi(n) 中最大 bad 质因子 q（即指数==1 的最大质数）”
    # 因为 DFS 后续只会加入更小的质数 p：
    #   - (p-1) 的质因子都 < p
    #   - p^(a-1) 只会增加质数 p 自己的指数
    # 因此如果某个大质数 q 已经出现且指数为 1，那么未来只靠更小的 p 无法再把 q 的指数抬高
    # => 后续允许选择的 p 必须 >= q，否则该分支必死
    #
    # 实现方式：用最大堆（Python heapq 是最小堆，所以存 -q），并 lazy 删除失效元素
    # ------------------------------------------------------------
    bad_heap = []

    def assign(p: int, new: int, log_old: dict, touched: list):
        """
        将 phi_exp[p] 设置为 new（new=0 表示删除该质因子）。
        这是 DFS 过程中“修改状态”的基础操作。

        参数：
        - p: 质数
        - new: 新指数（可为 0）
        - log_old: 当前 DFS 栈帧的“旧值记录字典”（只记录一次，便于回溯）
        - touched: 当前 DFS 栈帧中被修改过的质数列表（回溯时按逆序恢复）

        维护：
        - bad_cnt：指数==1 的质因子个数
        - bad_heap：当 new==1 时，把 p 推入堆（lazy 即可）
        """
        nonlocal bad_cnt

        old = phi_exp.get(p, 0)
        if old == new:
            return

        # 只在本 DFS 栈帧第一次触碰该质数 p 时，记录它的 old 值
        # 这样回溯时只要恢复到 old 即可，不需要记录每一次中间变化
        if p not in log_old:
            log_old[p] = old
            touched.append(p)

        # 去掉 old 对 bad_cnt 的贡献
        if old == 1:
            bad_cnt -= 1

        # 加上 new 对 bad_cnt 的贡献
        if new == 1:
            bad_cnt += 1
            # p 现在成为 bad（指数==1），候选最大 bad 质数可能更新
            heapq.heappush(bad_heap, -p)

        # 写回 phi_exp
        if new == 0:
            phi_exp.pop(p, None)
        else:
            phi_exp[p] = new

    def add_factors(fac_list, log_old, touched):
        """
        将 fac_list = [(q, e), ...] 加到 phi_exp 上
        即：phi_exp[q] += e

        fac_list 来自对 (p-1) 的分解，或其它需要累加的部分。
        """
        for q, e in fac_list:
            assign(q, phi_exp.get(q, 0) + e, log_old, touched)

    def restore(log_old: dict, touched: list):
        """
        回溯：将当前 DFS 栈帧所有 touched 的质数恢复到进入栈帧时的旧值 old

        注意：
        - 我们只保证恢复正确性，bad_heap 采用 lazy 策略：
          恢复时如果 old==1 也会 push 一次堆，可能造成重复；
          但查询最大 bad 质数时会验证 phi_exp[q]==1，不符合就弹出。
        """
        nonlocal bad_cnt

        for p in reversed(touched):
            old = log_old[p]
            cur = phi_exp.get(p, 0)
            if cur == old:
                continue

            # 去掉 cur 对 bad_cnt 的贡献
            if cur == 1:
                bad_cnt -= 1

            # 加上 old 对 bad_cnt 的贡献
            if old == 1:
                bad_cnt += 1
                heapq.heappush(bad_heap, -p)

            # 恢复 phi_exp
            if old == 0:
                phi_exp.pop(p, None)
            else:
                phi_exp[p] = old

    def current_max_bad_prime() -> int:
        """
        返回当前 phi(n) 中指数==1 的最大质数 q。
        若不存在 bad（bad_cnt==0），返回 2 作为默认下界。

        实现：
        - bad_heap 里存的是候选 bad 质数（可能已失效）
        - 反复检查堆顶 q 是否仍满足 phi_exp[q]==1
          - 若是，q 就是最大 bad 质数（因为堆按 -q 最小即 q 最大）
          - 若否，弹出并继续（lazy 清理）
        """
        while bad_heap:
            q = -bad_heap[0]
            if phi_exp.get(q, 0) == 1:
                return q
            heapq.heappop(bad_heap)
        return 2

    def phi_is_achilles() -> bool:
        """
        判断当前维护的 phi(n) 是否为 Achilles：

        1) powerful：所有指数 >= 2
           在本实现中等价于 bad_cnt == 0（不存在指数==1）
        2) not perfect power：所有指数 gcd == 1
        3) 不为单质数幂：质因子数 >= 2
           因为 p^e 一定是 perfect power，不可能是 Achilles

        注：len(phi_exp) < 2 在“正确性上”可以省略（因为 gcd 也会判掉），
            但保留它可作为 cheap early-exit，减少 gcd 遍历开销。
        """
        if bad_cnt != 0 or len(phi_exp) < 2:
            return False
        g = 0
        for e in phi_exp.values():
            g = gcd(g, e)
            if g == 1:
                return True
        return False

    ans = 0

    def dfs(now: int, last_idx: int, g_n: int):
        """
        深搜枚举 n：

        参数含义：
        - now：当前构造出的 n 的数值
        - last_idx：允许选择的质数下标范围是 primes[0:last_idx]
                    （我们从大到小选质数，保证不重不漏）
        - g_n：当前 n 的所有指数的 gcd（0 表示还没选任何质因子）
               由于所有指数 a>=2，n 是 Achilles <=> g_n == 1

        DFS 的一次扩展：
        - 选择一个质数 p（从大到小）
        - 先把 (p-1) 的质因子指数加入 phi_exp（对应 phi 里的乘子 (p-1)）
        - 再枚举 n 中 p 的指数 a = 2,3,4,...：
              n *= p^a
              phi(n) 额外乘 p^(a-1) -> phi_exp[p] 指数 + (a-1)
          这里我们用“每次 a 增加 1，就把 phi_exp[p] +1”的方式增量维护
        """
        nonlocal ans

        # 统计答案：
        # - now>1：排除空积 n=1
        # - g_n==1：n 不是 perfect power（因为指数 gcd==1），且 powerful 由构造保证
        # - phi(n) 也要 Achilles
        if now > 1 and g_n == 1 and phi_is_achilles():
            ans += 1

        if last_idx <= 0:
            return

        # 下一步至少要乘 p^2 才能保持 n powerful
        # 因此需要 now * p^2 <= N  -> p <= sqrt(N/now)
        limit_p = isqrt(N // now)
        if limit_p < 2:
            return

        # 我们只允许使用 <= primes[last_idx-1] 的质数（保持降序）
        max_p_val = min(limit_p, primes[last_idx - 1])
        max_idx = min(last_idx, bisect_right(primes, max_p_val))
        if max_idx <= 0:
            return

        # 剪枝：
        # 若当前 phi(n) 中存在最大 bad 质数 q（指数==1），
        # 则后续选取的 p 必须满足 p >= q，否则无法再增加 q 的指数（该分支必死）
        q = current_max_bad_prime()
        min_idx = bisect_left(primes, q)
        if min_idx >= max_idx:
            return

        # 从大到小选质数 p
        # 注意：这里为了效率，真正应该从 max_idx-1 下降到 min_idx
        # （否则会做很多无意义枚举）
        for idx in range(max_idx - 1, -1, -1):
            p = primes[idx]

            # -------------------------
            # 这一帧对 phi_exp 的改动日志：
            # log_old: 记录每个 touched 质数进入本帧时的 old 值
            # touched: 记录被修改过的质数，用于回溯
            # -------------------------
            log_old = {}
            touched = []

            # phi *= (p-1)：把 (p-1) 的分解因子加到 phi_exp
            add_factors(pm1_fac[idx], log_old, touched)

            # 枚举指数 a >= 2：
            # 维护 now * p^a（用 have 增量乘 p）
            # 同时 phi_exp[p] 每次 +1（因为 phi(n) 的 p 指数是 a-1）
            cur_p_exp = phi_exp.get(p, 0)

            have = now * p  # 当前 have = now * p^1
            a = 1
            while have <= N // p:
                have *= p
                a += 1  # 此时 have = now * p^a，且 a>=2

                # phi 里 p 的指数增加 1（对应 a-1 增加 1）
                cur_p_exp += 1
                assign(p, cur_p_exp, log_old, touched)

                # 更新 n 的指数 gcd：把新指数 a 纳入 gcd
                # g_n==0 表示第一次选质数：gcd(0, a) 视为 a
                dfs(have, idx, gcd(g_n, a))

            # 回溯：撤销本帧所有对 phi_exp / bad_cnt 的修改
            restore(log_old, touched)

    # 从空积开始：n=1，允许使用所有 primes，指数 gcd 初始化为 0
    dfs(1, Pn, 0)
    return ans


if __name__ == "__main__":
    print(solve_strong_achilles(N))