Project Euler 822

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

Project Euler 822

题目

Square the Smallest

A list initially contains the numbers \(2, 3, \dots, n\).

At each round, the smallest number in the list is replaced by its square. If there is more than one such number, then only one of them is replaced.

For example, below are the first three rounds for \(n = 5\): \[[2, 3, 4, 5] \xrightarrow{(1)} [4, 3, 4, 5] \xrightarrow{(2)} [4, 9, 4, 5] \xrightarrow{(3)} [16, 9, 4, 5].\]

Let \(S(n, m)\) be the sum of all numbers in the list after \(m\) rounds.

For example, \(S(5, 3) = 16 + 9 + 4 + 5 = 34\). Also \(S(10, 100) \equiv 845339386 \pmod{1234567891}\).

Find \(S(10^4, 10^{16})\). Give your answer modulo \(1234567891\).

解决方案

令\(P=1234567891\)。某个初始数 \(i\in\{2,\dots,n\}\)，如果总共被平方了 \(t\) 次，则它最终变成：

\[ i \xrightarrow{\square} i^2 \xrightarrow{\square} i^{2^2} \xrightarrow{\square}\cdots \xrightarrow{\square} i^{2^t}. \]

因此最终序列完全由每个 \(i\) 的平方次数 \(t_i\) 决定，顺序不重要。问题变为：求每个 \(i\) 在 \(m\) 轮内被选为“当前最小值”的次数 \(t_i\)。

直接比较巨大整数不现实，但比较大小只需单调变换即可。令\(f(x)=\log_2(\log_2 x),\)所以有\(f(x^2)=\log_2(\log_2(x^2))=\log_2(2\log_2 x)=1+\log_2(\log_2 x)=f(x)+1.\)

结论： 对某个元素做一次“平方”，等价于把它的 \(f\)函数值 加 1。

而每一轮选择当前最小的 \(x\)，由于 \(f\) 单调递增，这等价于：每一轮选择当前最小的 \(f\)-值，并把它加 1。于是整个过程变成了对一堆实数反复执行 “取最小 +1”。

当某一时刻满足\((\min\{(a_{n,k})^2\} \ge \max\{ a_{n,k}\},\)则把最小值平方后，它会变成不小于原最大值，因此它立刻成为新的最大值；与此同时新的最小值就是原来的第二小。

如果此时把序列按从小到大排序为\(x_1\le x_2\le \cdots \le x_{n-1},\)且满足 \(x_1^2\ge x_{n-1}\)，那么下一轮一定是平方 \(x_2\)，再下一轮平方 \(x_3\)，……，直到平方 \(x_{n-1}\)。这整段过程里每个元素恰好被平方一次，并且又回到同样结构。

最终，循环阶段条件 \(\min^2\ge \max\) 等价于\(f(\max)\le f(\min)+1\Longleftrightarrow\min f \ge \max f - 1.\)

初始最大值是 \(n\)，初始最大 \(f\)-值为\(f(n)=\log_2(\log_2 n).\)为了让 \(\min f \ge f(n)-1\)，只需要把所有小于阈值\(B=f(n)-1\)的元素的 \(f\)值补到至少 \(B\)。

对初始数 \(i\)，其初始 \(f(i)=\log_2(\log_2 i)\)。若它被平方 \(s_i\) 次，则 \(f\)值变为 \(f(i)+s_i\)。我们只需最小的整数 \(s_i\) 使得 \(f(i)+s_i\ge B.\)因此\(s_i=\max\left(0,\left\lceil B-f(i)\right\rceil\right)=\max\left(0,\left\lceil f(n)-1-f(i)\right\rceil\right).\)

令\(\displaystyle{k=\sum_{i=2}^n s_i,}\)则经过 \(k\) 轮后一定进入循环阶段。并且此时每个元素对应的真实值为\(i^{2^{s_i}}.\)

循环阶段周期长度是 \(N=n-1\)。剩余轮数为\(R=m-k.\)把它除以周期：\(R=qN+r, 0\le r<N.\)那么，每个元素都会额外被平方 \(q\) 次；按循环顺序的前 \(r\) 个元素再多平方 1 次。

于是总平方次数为\(t_i=s_i+q+\mathbf{1}[\text{pos}(i)<r],\)其中 \(\text{pos}(i)\) 是元素 \(i^{2^{s_i}}\) 在循环阶段的排序位置（从 0 开始）。

排序比较 \(i^{2^{s_i}}\) 的大小很难，但取对数即可：\(\log_2(i^{2^{s_i}})=2^{s_i}\log_2 i.\)因此只需要按\(2^{s_i}\log_2 i\) 升序排序即可得到循环顺序。

因此，对每个\(i\)，假设最终平方了\(t_i\)次，根据费马小定理，最终等价于计算\(i^{2^{t_i}}\equiv i^{2^{t_i}\bmod \varphi(P)}\pmod 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
from heapq import heappush, heappop, heapify
from math import log2, ceil
from tools import phi


N = 10 ** 4
M = 10 ** 16
MOD = 1234567891
PHI = phi(MOD)

def simulate_small(n, m):
    N = n - 1
    cnt = [0] * N
    h = [(log2(i), i - 2) for i in range(2, n + 1)]
    heapify(h)
    for _ in range(m):
        v, idx = heappop(h)
        v *= 2.0
        cnt[idx] += 1
        heappush(h, (v, idx))
    return cnt

def S(n, m):
    N = n - 1
    if m == 0:
        return sum(range(2, n + 1)) % MOD

    lg2 = [0.0] * N
    for i in range(2, n + 1):
        lg2[i - 2] = log2(i)

    vmax = log2(log2(n))
    B = vmax - 1.0

    s = [0] * N
    k = 0
    eps = 1e-12
    for i in range(2, n + 1):
        v = log2(log2(i))
        d = B - v
        if d > eps:
            si = int(ceil(d - eps))
            s[i - 2] = si
            k += si

    if m <= k:
        t = simulate_small(n, m)
        ans = 0
        for idx, ti in enumerate(t):
            base = idx + 2
            e = pow(2, ti, PHI)
            ans = (ans + pow(base, e, MOD)) % MOD
        return ans

    R = m - k
    q, r = divmod(R, N)

    order = list(range(N))
    order.sort(key=lambda idx: lg2[idx] * (1 << s[idx]))

    pos = [0] * N
    for p, idx in enumerate(order):
        pos[idx] = p

    ans = 0
    for idx in range(N):
        extra = 1 if pos[idx] < r else 0
        ti = s[idx] + q + extra
        base = idx + 2
        e = pow(2, ti, PHI)
        ans = (ans + pow(base, e, MOD)) % MOD
    return ans

print(S(N,M))
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
from heapq import heappush, heappop, heapify
from math import log2
from tools import phi

N = 10 ** 4
M = 10 ** 16
MOD = 1234567891
P = phi(MOD)
def S(n, m):
    a = list(range(2,n+1))
    o = len(a)
    q = [(log2(log2(x)),i) for i,x in enumerate(a)]
    cnt = [0 for _ in range(o)]
    heapify(q)
    mx = max(q)[0]
    while m > 0 and mx - q[0][0] > 1:
        u, idx = heappop(q)
        u += 1
        cnt[idx] += 1
        mx = max(mx, u)
        heappush(q, (u, idx))
        m -= 1
    q.sort()
    d,r = divmod(m,o)
    for i in range(o):
        idx = q[i][1]
        cnt[idx] += d + (i < r)
    ans = 0
    for i in range(o):
        ans += pow(a[i],pow(2, cnt[i], P),MOD)
    return ans % MOD
print(S(N, M))