Project Euler 548

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

Project Euler 548

题目

Gozinta Chains

A gozinta chain for \(n\) is a sequence \(\{1,a,b,\dots,n\}\) where each element properly divides the next.

There are eight gozinta chains for \(12:\)

\(\{1,12\} ,\{1,2,12\}, \{1,2,4,12\}, \{1,2,6,12\}, \{1,3,12\}, \{1,3,6,12\}, \{1,4,12\}\) and \(\{1,6,12\}\).

Let \(g(n)\) be the number of gozinta chains for \(n\), so \(g(12)=8\). \(g(48)=48\) and \(g(120)=132\).

Find the sum of the numbers \(n\) not exceeding \(10^{16}\) for which \(g(n)=n\).

解决方案

令\(N=10^{16}\)。对 \(n>1\)，任意一条以 \(n\) 结尾的链 \({1,\dots,d,n}\)，其倒数第二项 \(d\) 必为 \(n\) 的真因子（\(d\mid n,\ d<n\)）。去掉最后的 \(n\) 后，就得到一条以 \(d\) 结尾的 gozinta chain。

因此对 \(n>1\) 有递推：

\[ g(n)=\sum_{\substack{d\mid n\ d<n}} g(d), g(1)=1. \]

例如 \(n\) 为质数时只有 \({1,n}\) 一条链，所以 \(g(p)=1\)。

设\(n\)的分解为\(\displaystyle{n=\prod_{i=1}^k p_i^{e_i}.}\)此外，我们定义\(n\)的指数多重集\(\mathrm{sig}(n)\)表示\((e_1,e_2,\dots,e_k)\)按非增序排列得到的序列。

一个因子 \(d\mid n\) 等价于选择指数向量 \((a_1,\dots,a_k)\)，其中 \(0\le a_i\le e_i\)，对应 \(\displaystyle{d=\prod_{i=1}^k p_i^{a_i}}\)。

而 “严格整除” 等价于向量的逐坐标不减且不全相等：\((a_i)\prec(b_i)\) 当且仅当 \(a_i\le b_i\) 对所有 \(i\) 成立且至少一处严格小于。

所以 \(g(n)\) 只由指数集合 \({e_1,\dots,e_k}\)（按降序写成 \(e_1\ge e_2\ge\dots\ge e_k\)）决定，与具体用了哪些质数无关。记这个函数为\(g(e_1,\dots,e_k).\)

接下来，我们尝试枚举所有可能的指数签名，若指数签名为 \((e_1,\dots,e_k)\)，在所有具有该签名的数中，最小的一个一定是把大指数分配给小质数：\(\displaystyle{n_{\min}=\prod_{i=1}^k P(i)^{e_i}.}\)其中\(P(i)\)表示第\(i\)个质数。

若 \(n_{\min}>N\)，则任何具有该签名的 \(n\) 都会超界，因此该签名可以直接丢弃。

于是我们可以递归生成所有满足 \(n_{\min}\le N\) 的非增序列 \((e_1\ge e_2\ge\dots)\)。

但是，直接用\(\displaystyle{g(n)=\sum_{d\mid n,d<n} g(d)}\)会遇到因子数太多的问题。这里引入\(\displaystyle{G(n)=\sum_{d\mid n} g(d).}\)

对 \(n>1\)，由于 \(g(n)\) 恰好等于所有真因子的 \(g\) 之和，所以\(\displaystyle{G(n)=g(n)+\sum_{d\mid n,d<n} g(d)=2g(n),(n>1).}\)并且 \(G(1)=g(1)=1\)。

设 \(n\) 的不同质因子为 \({p_1,\dots,p_k}\)。一个因子 \(d\) 是真因子当且仅当它在至少一个质因子方向上少拿了指数，也就是对某个\(i\)，有\(d\mid \dfrac{n}{p_i}.\)

因此真因子集合是并集：\(\displaystyle{\{d\mid n,\ d<n\}=\bigcup_{i=1}^k \{d\mid n/p_i\}.}\)

对并集做容斥，就得到

\[ g(n)=\sum_{\emptyset\ne I\subseteq\{1,\dots,k\}}(-1)^{|I|+1}G\left(\frac{n}{\prod_{i\in I}p_i}\right). \]

对签名 \(\mathbf e=(e_1,\dots,e_k)\)，把除以若干个不同质因子理解成对若干个坐标各减 \(1\)”。记 \(\mathbf 1_I\) 为在 \(I\subseteq \{1,\dots,k\}\) 上为 \(1\)、其余为 \(0\) 的向量，对于非空签名，则有\(G(\mathbf e)=2g(\mathbf e)\)，因此可以写出递推式：

\[ G(\mathbf e)= \left \{\begin{aligned} &1 & & \text{if}\quad e=\varnothing \\ &2\sum_{\emptyset\ne I\subseteq\{1,\dots,k\}}(-1)^{|I|+1}G(\mathbf e-\mathbf 1_I) & & \text{else} \end{aligned}\right. \]

实现时把 \(\mathbf e-\mathbf 1_I\) 做完后重新排序并去掉 \(0\)，即可回到签名的规范形。

最终我们对每个签名 \(\mathbf e\) 计算 \(m=g(\mathbf e)\)。

  • 若 \(m>N\)：不可能是答案，跳过。
  • 若 \(m\le N\)：再分解 \(m\)，取其指数多重集 \(\mathrm{sig}(m)\)。如果 \(\mathrm{sig}(m)=\mathbf e\)，则 \(g(m)=g(\mathbf e)=m\)，说明 \(m\) 本身就是一个满足条件的 \(n\)。

把所有这样的 \(m\) 去重求和即可（并注意 \(n=1\) 也满足 \(g(1)=1\)）。

代码

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
import itertools
from functools import lru_cache
from math import comb

from tools import factorization, Factorizer

N = 10 ** 16

# ---------------------------
# prime utilities
# ---------------------------

PRIMES = Factorizer(100).primes

# ---------------------------
# generate all exponent signatures e1>=e2>=... with product(prod p_i^e_i) <= LIMIT
# where primes are assigned in increasing order (minimizes the number for a signature)
# ---------------------------
def gen_signatures(limit=N, primes=PRIMES):
    sigs = []

    def dfs(i: int, max_exp: int, cur_prod: int, exps: list):
        if exps:
            sigs.append(tuple(exps))
        if i >= len(primes):
            return
        p = primes[i]
        prod = cur_prod
        for e in range(1, max_exp + 1):
            prod *= p
            if prod > limit:
                break
            exps.append(e)
            dfs(i + 1, e, prod, exps)
            exps.pop()

    dfs(0, 60, 1, [])
    return sigs

SIGS = gen_signatures()

# ---------------------------
# G(signature) recursion with grouping (multiset compression)
# signature is a tuple of positive ints in non-increasing order.
# empty tuple corresponds to n=1.
# Define G(n)=sum_{d|n} g(d). For n>1, G(n)=2g(n); and G(1)=1.
# Recurrence:
#   G(e)=2 * sum_{nonempty I} (-1)^{|I|+1} G(e - 1_I)
# where e-1_I means decrement selected coordinates by 1 then re-sort and drop zeros.
# ---------------------------
def groups(sig):
    res = []
    i = 0
    while i < len(sig):
        v = sig[i]
        j = i
        while j < len(sig) and sig[j] == v:
            j += 1
        res.append((v, j - i))
        i = j
    return res

@lru_cache(maxsize=None)
def G(sig):
    if not sig:
        return 1
    gs = groups(sig)
    ranges = [range(c + 1) for v, c in gs]

    total = 0
    for ts in itertools.product(*ranges):
        s = sum(ts)
        if s == 0:
            continue

        coeff = 1
        for (v, c), t in zip(gs, ts):
            coeff *= comb(c, t)

        # sign = (-1)^{s+1}: if s odd -> +, if s even -> -
        if s % 2 == 0:
            coeff = -coeff

        child = []
        for (v, c), t in zip(gs, ts):
            if c - t:
                child.extend([v] * (c - t))
            if t and v > 1:
                child.extend([v - 1] * t)

        child.sort(reverse=True)
        total += coeff * G(tuple(child))

    return 2 * total

def g_from_sig(sig):
    return 1 if not sig else G(sig) // 2


def signature_of(n: int):
    return tuple(sorted([v[1] for v in factorization(n)],reverse=True))
    # return tuple(exps)

# ---------------------------
# main
# ---------------------------
def solve():
    sols = {1}
    for sig in SIGS:
        val = g_from_sig(sig)
        if val > N:
            continue
        if signature_of(val) == sig:
            sols.add(val)
    return sum(sols)

if __name__ == "__main__":
    print(solve())  # 12144044603581281