Project Euler 823

题目

Factor Shuffle

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

At each round, every number in the list is divided by its smallest prime factor. Then the product of these smallest prime factors is added to the list as a new number. Finally, all numbers that become \(1\) are removed from the list.

For example, below are the first three rounds for \(n = 5\):

\[[2, 3, 4, 5] \xrightarrow{(1)} [2, 60] \xrightarrow{(2)} [30, 4] \xrightarrow{(3)} [15, 2, 4].\]

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

For example, \(S(5, 3) = 15 + 2 + 4 = 21\). Also \(S(10, 100) = 257\).

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

解决方案

对任意正整数 \(x\)，把它分解为质因子并按从大到小排序：\(\displaystyle{x=\prod_{t=1}^{\Omega(x)} q_t, q_1\ge q_2\ge \cdots \ge q_{\Omega(x)}.}\)其中 \(\Omega(x)\) 表示质因子个数（带重数）。例如 \(60=2^2\cdot 3\cdot 5 \Rightarrow [5,3,2,2].\)

注意到：\(x\) 的最小质因子正好是这个序列的最后一个元素，因此除以最小质因子等价于 弹出这个序列末尾一个元素。

因此把列表中的每个数都替换成它的质因子序列看成一个栈，每一轮操作就变成：对每个栈弹出一个元素（它的最小质因子）；把弹出的所有质因子收集起来，排序后作为一个新栈加入；空栈删掉（对应数变成 \(1\) 被移除）。

于是整道题完全等价于对一堆质因子栈做操作，最终每个栈代表的数就是该栈元素乘积。

令当前有 \(r\) 个栈，其长度分别为\((\ell_1,\ell_2,\dots,\ell_r), \ell_i\ge 1.\)

一轮操作会从每个栈弹出一个质因子，所以每个长度减 \(1\)，并新加一个栈长度为 \(r\)（因为弹出了 \(r\) 个质因子）：\((\ell_1,\dots,\ell_r)\longrightarrow(\ell_1-1,\dots,\ell_r-1,r),\)然后把减到 \(0\) 的删掉。

这正是经典的 Bulgarian solitaire 的标准版过程。

关键不变量：总纸牌数（即所有栈元素总数）保持不变：\(\displaystyle{P=\sum_{x=2}^{n}\Omega(x).}\)把三角数记为\(T_k=\dfrac{k(k+1)}{2}.\)取最小的 \(k\) 使得\(T_{k-1}<P\le T_k.\)

在这个过程里，经过不算太多轮之后（经验上至多 \(k(k-1)\) 轮），栈的结构会进入一种接近三角形的稳定形态： 共有 \(k\) 个栈，第 \(i\) 行（从 \(i=0\) 开始）长度只可能是 \(i\) 或 \(i+1\)，也就是最多缺一个格子。

对于本题 \(n\)，\(k\) 的规模大约在\(O(\sqrt{n\log n})\)，因此直接模拟 \(k(k-1)\) 轮非常轻松，足以保证进入稳定形态。

稳定后可以把第 \(i\) 行补齐成长度 \(i+1\) 的形式：若该行原本长度是 \(i\)，就在最左边（栈底）补一个 \(1\)；若长度是 \(i+1\)，不补。

令补齐后的三角表元素为 \(x_0(i,j)\)，其中\(0\le j\le i < k.\)这一行代表的数就是\(\displaystyle{A_0(i)=\prod_{j=0}^{i} x_0(i,j).}\)因为补进去的只是 \(1\)，不会改变乘积。

进入稳定三角形后，每一轮操作的效果可以描述成：对固定的列 \(j\)，它存在于行 \(j,j+1,\dots,k-1\)，共有 \(k-j\) 个位置；并且每进行一轮，这一列的元素会做一次循环位移（本质上是被弹出的最小质因子重新组成新行后回到表里）。

因此第 \(j\) 列的周期就是\(k-j.\)于是进行了 \(M\) 轮之后，任意位置 \((i,j)\) 的值满足：\(x_M(i,j)=x_0(((i-j+M)\bmod (k-j))+j,j).\)

这一步非常重要：它把多个轮次的过程变成了简单取模下标。

至此，每一行 \(i\) 代表一个数\(\displaystyle{A_M(i)=\prod_{j=0}^{i} x_M(i,j).}\)如果这一行全都是 \(1\)，则它对应的数为 \(1\)，在真实过程里会被移除，因此贡献为 \(0\)。最终有：

\[S(n,m)\equiv \sum_{i=0}^{k-1} A_M(i)\]

代码

from tools import int_sqrt

N = 10 ** 4
M = 10 ** 16
MOD = 1234567891


def solve(N: int, M: int) -> int:
    spf = list(range(N + 1))
    for i in range(2, int(N ** 0.5) + 1):
        if spf[i] == i:
            for j in range(i * i, N + 1, i):
                if spf[j] == j:
                    spf[j] = i

    T = []
    P = 0
    for x in range(2, N + 1):
        fs = []
        y = x
        while y > 1:
            p = spf[y]
            fs.append(p)
            y //= p
        fs.sort(reverse=True)
        T.append(fs)
        P += len(fs)

    k = (int_sqrt(1 + 8 * P) - 1) // 2
    if k * (k + 1) // 2 < P:
        k += 1

    t = min(M, k * (k - 1))
    for _ in range(t):
        row = [pile[-1] for pile in T]
        row.sort(reverse=True)
        for pile in T:
            pile.pop()
        T = [pile for pile in T if pile]
        T.append(row)
    M -= t

    if M == 0:
        s = 0
        for pile in T:
            a = 1
            for p in pile:
                a = (a * p) % MOD
            s = (s + a) % MOD
        return s

    k = len(T)
    full = [len(T[i]) == i + 1 for i in range(k)]

    def x0(i: int, j: int) -> int:
        row = T[i]
        if full[i]:
            return row[j]
        if j == 0:
            return 1
        return row[j - 1]

    s = 0
    for i in range(k):
        prod = 1
        all1 = True
        for j in range(i + 1):
            modv = k - j
            r = ((i - j + (M % modv)) % modv) + j
            x = x0(r, j)
            if x != 1:
                all1 = False
            prod = (prod * x) % MOD
        if not all1:
            s = (s + prod) % MOD

    return s


print(solve(N, M))