Project Euler 386

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

Project Euler 386

题目

Maximum length of an antichain

Let \(n\) be an integer and \(S(n)\) be the set of factors of \(n\).

A subset \(A\) of \(S(n)\) is called an antichain of \(S(n)\) if \(A\) contains only one element or if none of the elements of \(A\) divides any of the other elements of \(A\).

For example: \(S(30) = \{1, 2, 3, 5, 6, 10, 15, 30\}\)

\(\{2, 5, 6\}\) is not an antichain of \(S(30)\).

\(\{2, 3, 5\}\) is an antichain of \(S(30)\).

Let \(N(n)\) be the maximum length of an antichain of \(S(n)\).

Find \(\sum N(n)\) for \(1 \le n \le 10^8\)

解决方案

令\(M=10^8\)。令\(n\)的质因子分解为\(\displaystyle{n=\prod_{i=1}^{r} p_i^{e_i}}\)，那么任意因子 \(d\mid n\) 都可以写成\(\displaystyle{d=\prod_{i=1}^{r} p_i^{a_i}, 0\le a_i\le e_i}\)，于是每个因子对应一个指数向量\((a_1,a_2,\dots,a_r)\in [0,e_1]\times[0,e_2]\times\cdots\times[0,e_r]\)，并且\(d_1\mid d_2 \iff a_i^{(1)}\le a_i^{(2)}(\forall i)\)。所以 \((S(n),\mid)\) 与链的直积同构，是一个典型的分级偏序。

可见，反链的一个天然来源是固定总指数和。定义因子 \(\displaystyle{d=\prod_{i=1}^{r} p_i^{a_i}}\) 的度数（也就是质因子个数带重数）：\(\displaystyle{\Omega(d)=\sum_{i=1}^r a_i.}\)如果两个因子 \(d_1\mid d_2\) 且 \(d_1\ne d_2\)，那么必有某个坐标严格变大，所以\(\Omega(d_1)<\Omega(d_2)\)。

因此所有满足 \(\Omega(d)=k\) 的因子构成一个反链。令\(R_k(n)=\#\{d\mid n:\Omega(d)=k\}\)

那么任何“固定 \(k\)”的集合大小就是 \(R_k(n)\)，从而\(\displaystyle{N(n)\ge \max_k R_k(n)}\)。接下来要证明其实可以取等号。

这里的偏序本质是链的直积：\([0,e_1]\times\cdots\times[0,e_r]\)，这类偏序满足 Sperner 性质：其最大反链大小等于某一层（rank level）的最大大小，也就是\(\displaystyle{N(n)=\max_k R_k(n)}\)。直观理解是：在这种坐标逐个变大的结构里，反链最宽的地方就在中间层。

接下来使用生成函数计算每一层的大小。注意 \(R_k(n)\) 只与指数 \((e_1,\dots,e_r)\) 有关，与具体质数 \(p_i\) 无关。

因为每个坐标 \(a_i\) 可以取 \(0,1,\dots,e_i\)，所以计数满足\(\displaystyle{\sum_{k\ge 0} R_k(n)x^k=\prod_{i=1}^r\left(\sum_{j=0}^{e_i}x^j\right)}\)。因此 \(R_k(n)\) 就是多项式\(\displaystyle{F(x)=\sum_{k\ge 0} R_k(n)x^k=\prod_{i=1}^r\left(\sum_{j=0}^{e_i}x^j\right)}\)的 \(x^k\) 系数。

于是：\(\displaystyle{N(n)=\max_k [x^k]F(x).}\)而 \(F(x)\) 的总次数是\(\displaystyle{M=\sum_{i=1}^r e_i=\Omega(n)}\)，并且由于映射 \(d\mapsto n/d\) 会把 \(\Omega(d)=k\) 映到 \(\Omega(n/d)=M-k\)，可得系数满足对称性：\(R_k(n)=R_{M-k}(n)\)，因此最大值一定在中间附近（\(\lfloor M/2\rfloor\) 或 \(\lceil M/2\rceil\)）。

因为 \(N(n)\) 只依赖于指数的多重集\(E_n= \{e_1,e_2,\dots,e_r\},e_1\ge e_2\ge \dots\ge e_r\)，所以我们可以把 \(1\le n\le 10^8\) 按指数类型分类。记 \(C(E)\)：满足该指数类型的整数个数，记 \(W(E)\)：该类型对应的最大层宽（即最大系数）

那么答案就是

\[ \sum_{n= 1}^{M} N(n)=\sum_{n=1}^M C(E_n)\cdot W(E_n) \]

其中\(\displaystyle{W(E)=\max_k [x^k]\prod_{e\in E}(1+x+\cdots+x^e).}\)因为 \(\Omega(n)\le \log_2 N\)，所以每个 \(E\) 的指数和非常小，计算 \(W(E)\) 用简单卷积就够了。

然而直接对 \(1,\dots,M\) 分解质因数显然太慢。核心技巧是：我们只枚举那些会出现的“指数结构”，并用质数计数函数 \(\pi(x)\) 跳过大量情况。

如果某个质数 \(p>\sqrt{r}\)，那么在任意 \(\le r\) 的分解里，它的指数必然是 \(1\)（因为 \(p^2>r\)）。

这意味着：指数 \(\ge 2\) 的质因子一定来自较小的质数，指数 \(= 1\) 的末尾质因子可以用 \(\pi(r)\) 一次性计数。

构造集合\(V=\{ \lfloor N/1\rfloor,\lfloor N/2\rfloor,\dots\}\cup\{1,2,\dots,\lfloor\sqrt N\rfloor\}\)，只需要在这些点上知道 \(\pi(v)\)，就足够支撑递归中出现的所有上界 \(r\)。

可以用经典的批量 Legendre 筛法在 \(O(\sqrt N)\) 级别时间里算完 \(V\) 上的 \(\pi\)。

在递归状态 \((i, E, r)\) 中：

  • 当前已选质数都严格递增（用 \(i\) 控制）
  • 还允许乘进去的最大剩余值是 \(r\)
  • 当前指数多重集是 \(E\)

接下来可能新增一个质数 \(p\)：

  1. 指数为 \(1\) 且它是最后一个质因子：所有可选质数个数就是 \(\pi(r)-i\)，一次性加到类型 \(E\cup\{1\}\)。
  2. 指数为 \(1\) 但后面还有别的质因子：这时最小还要再乘一个 \(\ge p\) 的质数，所以必须有 \(p^2\le r\)，因此只对 \(p\le\sqrt r\) 的小质数展开递归即可。
  3. 指数为 \(k\ge 2\)： 同理必须满足 \(p^k\le r\)，并且剩余继续递归。

这样就能把所有 \(\le N\) 的数按指数类型统计出来。

代码

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
from collections import defaultdict
from tools import Factorizer, int_sqrt

N = 10 ** 8

def pi_table(n: int):
    r = int_sqrt(n)
    V = [n // i for i in range(1, n // r + 1)]
    V += list(range(r - 1, 0, -1))
    C = {v: v - 1 for v in V}
    for q in range(2, r + 1):
        if C[q] > C[q - 1]:
            pi = C[q - 1]
            q2 = q * q
            for v in V:
                if v < q2:
                    break
                C[v] -= C[v // q] - pi
    return C


def factorization_count(n: int):
    pi = pi_table(n)
    P = Factorizer(max(100, int_sqrt(n) + 5)).primes
    T = defaultdict(int)
    T[()] = 1

    def rec(idx: int, exp: tuple, r: int):
        w = pi[r] - idx
        if w <= 0:
            return
        E1 = tuple(sorted(exp + (1,)))
        T[E1] += w

        for i in range(idx, len(P)):
            p = P[i]
            if p * p > r:
                break
            rec(i + 1, E1, r // p)

            q = p * p
            k = 2
            while q <= r:
                Ek = tuple(sorted(exp + (k,)))
                T[Ek] += 1
                if q * p <= r:
                    rec(i + 1, Ek, r // q)
                k += 1
                q *= p

    rec(0, (), n)
    return T


def width_from_exponents(E: tuple):
    coeffs = [1]
    for e in E:
        new = [0] * (len(coeffs) + e)
        for i, a in enumerate(coeffs):
            if a:
                for j in range(e + 1):
                    new[i + j] += a
        coeffs = new
    return max(coeffs)


def solve(N):
    T = factorization_count(N)
    cache = {}
    ans = 0
    for E, cnt in T.items():
        w = cache.get(E)
        if w is None:
            w = width_from_exponents(E)
            cache[E] = w
        ans += w * cnt
    return ans


print(solve(N))