Project Euler 858

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

Project Euler 858

题目

LCM

Define \(G(N) = \sum_S \operatorname{lcm}(S)\) where \(S\) ranges through all subsets of \(\{1, \dots, N\}\) and \(\operatorname{lcm}\) denotes the lowest common multiple. Note that the \(\operatorname{lcm}\) of the empty set is \(1\).

You are given \(G(5) = 528\) and \(G(20) = 8463108648960\).

Find \(G(800)\). Give your answer modulo \(10^9 + 7\).

解决方案

由题目可知，因为元素 \(1\) 对最小公倍数没有影响，任意 \(T\subseteq\{2,\dots,N\}\) 对应两种子集（含 1 / 不含 1）且 \(\operatorname{lcm}\) 相同，所以

\[ G(N)=2\sum_{T\subseteq\{2,\dots,N\}} \operatorname{lcm}(T). \]

因此只需计算右侧，再乘 \(2\)。

若按精确 \(\operatorname{lcm}\) 值计算每个值产生的数量，加入新数 \(x\) 时每个旧状态 \(L\) 会产生不选 \(x\)：仍为 \(L\)；选 \(x\)：变为 \(\operatorname{lcm}(L,x)\)。这在定义上很自然，但状态数（不同 \(\operatorname{lcm}\) 的个数）增长过快，无法直接承受。

为此，我们只保留对未来有用的质因子指数。我们按某个处理顺序逐个加入整数。对某一时刻，记尚未处理的整数集合为 \(R\)。定义每个质数的未来上界指数\(\displaystyle B_R(q)=\max_{m\in R}\{v_q(m)\},\)并定义未来模版\(\displaystyle M_R=\prod_{q} q^{B_R(q)}.\)

对已选子集 \(T\)（并且其\(\operatorname{lcm}\)为\(L\)），那么只记录投影键\(\kappa_R(T)=\gcd(L,M_R).\)

引理 1（投影充分性）：对任意未来要加入的整数 \(x\in R\)，有\(\gcd(L,x)=\gcd(\kappa_R(T),x).\)

证明： 逐质数看指数。对任意质数 \(q\)：\(v_q(x)\le B_R(q),\)所以\(\min(v_q(L),v_q(x))=\min(\min(v_q(L),B_R(q)),v_q(x))=v_q(\gcd(\kappa_R(T),x)).\)各质数同时成立，即得结论。

这个结论说明：未来每一步转移里真正需要的 \(\gcd(L,x)\)，只靠投影键就能还原，不需要完整 \(L\)。

因此，我们以加权和而非计数存状态，令状态值为该类子集 \(\operatorname{lcm}\) 的总和：

\[ F_R(k)=\sum_{\kappa_R(T)=k}\operatorname{lcm}(T). \]

加入新整数 \(x\) 时，对任意键 \(k\)：

  1. 不选 \(x\)：贡献原值，键按新的 \(M\) 重新投影；
  2. 选 \(x\)：根据引理1，每个 \(\operatorname{lcm}\) 乘上\(\dfrac{x}{\gcd(\operatorname{lcm}(T),x)}=\dfrac{x}{\gcd(k,x)},\)

故整类总和整体乘同一系数。因此单步更新完全可在键与组和层面进行，不需要展开到单个子集。

因此，我们将 \(2,\dots,N\) 按最大质因子从大到小分层处理。好处是：当处理完质数 \(p\) 这一层后，剩余未处理整数都不再含 \(p\)。于是 \(p\) 的指数对未来已无作用，可以从投影模版 \(M_R\) 中彻底删去。

因为\(p\)来说，后续每个 \(x\) 都满足 \(v_p(x)=0\)。由引理 1，后续转移只依赖 \(\gcd(\cdot,x)\)；而 \(p\)-幂对该 gcd 永远不起作用，删去不影响任何未来转移。

代码

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
from collections import defaultdict
from tools import gcd

N = 800
MOD = 10**9 + 7
def update(dp, n, cap, mod):
    ndp = defaultdict(int)
    if mod is None:
        for l, v in dp.items():
            ndp[gcd(cap, l)] += v
            ng = n // gcd(l, n)
            ndp[gcd(cap, l * ng)] += v * ng
    else:
        for l, v in dp.items():
            a = gcd(cap, l)
            ndp[a] = (ndp[a] + v) % mod
            ng = n // gcd(l, n)
            b = gcd(cap, l * ng)
            ndp[b] = (ndp[b] + v * ng) % mod
    return ndp

def G(N):
    fac = list(range(N + 1))
    cap = 1
    for p in range(2, N + 1):
        if fac[p] != p:
            continue
        for d in range(2 * p, N + 1, p):
            if fac[d] == d:
                fac[d] = p
        t = N // p
        while t >= p:
            t //= p
            cap *= p

    dp = defaultdict(int)
    dp[1] = 1

    for p in range(N, 1, -1):
        if fac[p] != p:
            continue

        cap *= p

        for d in range((N // p) * p, p, -p):
            if fac[d]:
                dp = update(dp, d, cap, MOD)
                fac[d] = 0

        while cap % p == 0:
            cap //= p

        dp = update(dp, p, cap, MOD)

    s = sum(dp.values())
    return (s * 2) % MOD


print(G(N))