Project Euler 615

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

Project Euler 615

题目

The millionth number with at least one million prime factors

Consider the natural numbers having at least \(5\) prime factors, which don’t have to be distinct.

Sorting these numbers by size gives a list which starts with:

\(\begin{aligned} &32=2\cdot2\cdot2\cdot2\cdot2\\ &48=2\cdot2\cdot2\cdot2\cdot3\\ &64=2\cdot2\cdot2\cdot2\cdot2\cdot2\\ &72=2\cdot2\cdot2\cdot3\cdot3\\ &80=2\cdot2\cdot2\cdot2\cdot5\\ &96=2\cdot2\cdot2\cdot2\cdot2\cdot3\\ &\dots \end{aligned}\)

So, for example, the fifth number with at least \(5\) prime factors is \(80\).

Find the millionth number with at least one million prime factors.

Give your answer modulo \(123454321\).

解决方案

设质因子个数（含重）为\(\displaystyle{\Omega(n)=\sum_{p^e\parallel n} e.}\)令\(M=10^6\)表示质因子个数的下界，\(K=10^6\)表示题目所要求的序数。

显然满足 \(\Omega(n)\ge M\) 的最小数是\(2^M\)。因为 \(2\) 是最小的质数，想让 \(\Omega(n)\) 尽量大且 \(n\) 尽量小，就应该全用 \(2\)。

我们用“质数指数向量”表示一个数：设第 \(k\) 个质数为 \(p_k\)（\(p_0=2,p_1=3,p_2=5,\dots\)），用\(\mathbf e=(e_0,e_1,e_2,\dots)\) 表示\(\displaystyle{n(\mathbf e)=\prod_{k\ge 0} p_k^{e_k}.}\)

初始状态就是\(\mathbf e^{(0)}=(M,0,0,\dots)\)对应 \(2^{M}\)。

为了按从小到大的顺序枚举所有 \(\Omega(n)\ge 10^6\) 的数，我们考虑两种对 \(n\) 的“最小增量”变换：

第一类操作是把 \(n\) 乘上 \(2\)：\((e_0,e_1,\dots)\to (e_0+1,e_1,\dots).\)这会使\(\Omega(n)\to \Omega(n)+1\)且数值乘以 \(2\)，是最便宜的增加质因子个数的方式。

第二类操作是把一个质数替换成下一个质数：从分解里取出一个 \(p_k\)，换成 \(p_{k+1}\)：\(e_k\to e_k-1,e_{k+1}\to e_{k+1}+1.\)对应数值变化为\(n\to n\cdot \dfrac{p_{k+1}}{p_k}.\)这会保持 \(\Omega(n)\) 不变，但让 \(n\) 变大一些；并且在保持因子数不变的所有变换里，这是非常自然的一种单步增长。

可见，从 \(2^{M}\) 出发足够覆盖全部答案。原因在于，任取一个满足 \(\Omega(n)\ge M\) 的数 \(n\)，设\(\Omega(n)=10^6+t(t\ge 0).\)把其中额外的 \(t\) 个质因子尽量“拿成 2”，写成\(n=2^t\cdot m,\Omega(m)=M.\)

而 \(m\) 是恰好含 \(M\) 个质因子的数，意味着它可以写成 \(M\) 个质数的乘积。由于任意质数 \(p_k\) 都能通过反复执行“替换到下一个质数”（从 \(2=p_0\) 一路替换到 \(p_k\)）得到，所以 \(m\) 一定能从 \(2^{M}\) 出发用若干次第二类操作得到；然后再用 \(t\) 次第一类操作乘上 \(2^t\) 即可得到 \(n\)。

直接比较大整数会很慢，我们用对数对状态 \(\mathbf e\) 定义权值\(\displaystyle{\mathrm{val}(\mathbf e)=\log n(\mathbf e)=\sum_{k\ge 0} e_k\log p_k.}\)

两种操作都会使 \(n\) 变大，所以也都会让 \(\mathrm{val}\) 变大（权值严格上升）。

于是我们可以用一个最小堆做最短路式枚举（本质是BFS）：

  • 堆中始终存放尚未输出的候选状态，键为 \(\mathrm{val}(\mathbf e)\)；
  • 每次弹出堆顶，就是目前还没输出的最小 \(n\)；
  • 对该状态执行所有可能的两类操作，生成新状态入堆并去重。

这样弹出的序列就按 \(n\) 从小到大排列。

当我们弹出第 \(K\) 个状态时，对应的 \(n\) 就是题目要的第 \(K\) 个数。

代码

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
from heapq import heappush, heappop
from bisect import bisect_left

from tools import *

MOD = 123454321
M = 10**6   # 至少 1e6 个质因子
K = 10**6   # 第 1e6 个数


# 经验上只会用到很靠前的质数，下界给足即可（> 23500 个质数就很稳）
# 这里取 400000，质数数量约 33860
PRIMES = Factorizer(4 * 10 ** 5).primes
LOGP = [math.log(p) for p in PRIMES]
LOG2 = LOGP[0]

# 为了在 visited 里省内存，用无冲突整数编码状态
# 只要 base > 最大可能出现的 (index, exponent, e0) 即可
BASE = 2_000_003


def code_state(e0: int, vec: tuple):
    """
    vec: ((idx1, exp1), (idx2, exp2), ...) 只记录 idx>=1 的非零项，按 idx 升序
    编码为一个大整数，确保无碰撞。
    """
    x = e0
    for idx, exp in vec:
        x = x * BASE + idx
        x = x * BASE + exp
    return x


def vec_inc(vec: tuple, idx: int, delta: int):
    """
    在 vec 上对某个 idx>=1 的指数加 delta（可能为负），返回新的 vec（仍有序且无 0 项）
    """
    lst = list(vec)
    i = bisect_left(lst, (idx, 0))
    if i < len(lst) and lst[i][0] == idx:
        ne = lst[i][1] + delta
        if ne == 0:
            lst.pop(i)
        else:
            lst[i] = (idx, ne)
    else:
        # 原来没有该 idx
        if delta != 0:
            lst.insert(i, (idx, delta))
    return tuple(lst)


# 初始状态：2^M
e0 = M
vec = tuple()  # 只记录 idx>=1 的非零指数
logv = e0 * LOG2

hq = [(logv, e0, vec)]
vis = set()
vis.add(code_state(e0, vec))

# 弹出前 K-1 个，下一次弹出即第 K 个
for _ in range(K - 1):
    logv, e0, vec = heappop(hq)

    # 操作 A：乘以 2
    ne0 = e0 + 1
    nvec = vec
    c = code_state(ne0, nvec)
    if c not in vis:
        vis.add(c)
        heappush(hq, (logv + LOG2, ne0, nvec))

    # 操作 B(0)：把一个 2 变成 3
    if e0 > 0:
        ne0 = e0 - 1
        nvec = vec_inc(vec, 1, +1)  # idx=1 对应质数 3
        c = code_state(ne0, nvec)
        if c not in vis:
            vis.add(c)
            heappush(hq, (logv - LOG2 + LOGP[1], ne0, nvec))

    # 操作 B(k>=1)：把一个 p_k 变成 p_{k+1}
    for idx, exp in vec:
        # idx -> idx+1
        nvec1 = vec_inc(vec, idx, -1)
        nvec2 = vec_inc(nvec1, idx + 1, +1)
        c = code_state(e0, nvec2)
        if c not in vis:
            vis.add(c)
            heappush(hq, (logv - LOGP[idx] + LOGP[idx + 1], e0, nvec2))

# 第 K 个状态
logv, e0, vec = heappop(hq)

# 计算模值
ans = pow(2, e0, MOD)
for idx, exp in vec:
    ans = (ans * pow(PRIMES[idx], exp, MOD)) % MOD
print(ans)