Project Euler 418

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

Project Euler 418

题目

Factorisation triples

Let \(n\) be a positive integer. An integer triple \((a, b, c)\) is called a factorisation triple of \(n\) if: - \(1 \le a \le b \le c\)

  • \(a\cdot b\cdot c = n\).

Define \(f(n)\) to be \(a + b + c\) for the factorisation triple \((a, b, c)\) of \(n\) which minimises \(c / a\). One can show that this triple is unique.

For example, \(f(165) = 19, f(100100) = 142\) and \(f(20!) = 4034872\).

Find \(f(43!)\).

解决方案

约束 \(abc=n\) 且 \(a\le b\le c\) 下，若三者偏离很大（例如 \(a\) 很小而 \(c\) 很大），则 \(\dfrac{c}{a}\) 会被放大。反过来，如果 \(a,b,c\) 尽量接近，则 \(c/a\) 趋近于 \(1\)，显然更优。

用连续化的直觉：在实数域中，若仅要求 \(xyz=n\)，则 \(x=y=z=n^{1/3}\) 是最均衡的点。因此整数最优三元组通常也会出现在 \(n^{1/3}\) 的附近。记\(r=n^{1/3}.\)后续的枚举策略会围绕 \(r\) 设置一个很窄的相对区间\(\left[\dfrac{r}{\varepsilon},r\varepsilon\right]\) 来搜寻候选因子。这里 \(\varepsilon>1\) 是可调参数。

对任意可行三元组，定义\(q=\dfrac{n}{b}.\)则有\(ac=q.\)在 \(b\) 固定时，\(q\) 也固定。此时\(\dfrac{c}{a}=\dfrac{q/a}{a}=\dfrac{q}{a^2}.\)

这给出一个关键单调性：当 \(q\) 固定时，\(\dfrac{q}{a^2}\) 随 \(a\) 增大而严格减小。

因此在固定 \(b\) 的前提下，让 \(c/a\) 最小”等价于“让 \(a\) 尽可能大，但仍需满足排序约束 \(a\le b\le c\)。

因为 \(c=q/a\)，条件 \(b\le c\) 等价于\(b\le \dfrac{q}{a}\iff a\le \dfrac{q}{b}.\)再加上 \(a\le b\)，得到\(a\le \min\left(b,\left\lfloor\dfrac{q}{b}\right\rfloor\right).\)

因此，对固定 \(b\)，最优 \(a\) 必为所有满足\(a\mid q\)和\(a\le \min\left(b,\left\lfloor\dfrac{q}{b}\right\rfloor\right)\)的 \(a\) 中的最大者。随后令 \(c=q/a\) 自动满足 \(c\ge b\)。

于是三重优化被压缩为：

  1. 枚举候选的 \(b\)；
  2. 对每个 \(b\) 找到最大的可行 \(a\)；
  3. 比较得到全局最优。

如果 \((a,b,c)\) 为最优，按均衡性直觉三者应接近 \(r\)。在实际实现中，做法是：

  • 只收集所有满足\(\dfrac{r}{\varepsilon}\le d\le r\varepsilon\)的约数 \(d\mid n\) 作为候选集合 \(C\)；
  • 在 \(C\) 上枚举 \(b\) 并在同一集合中寻找 \(a\)。

这一步需要一个明确的假设：最优三元组 \((a,b,c)\) 的三条边（至少 \(a\) 与 \(b\)，通常也包括 \(c\)）都能在窗口内找到相应约数候选，从而能在候选集合中被枚举到。

同时，为避免浮点误差，立方根与区间端点应尽量用整数方法构造。可取整数立方根\(r_0=\left\lfloor n^{1/3}\right\rfloor\) 并再设置\(L=\left\lfloor \dfrac{r_0}{\varepsilon}\right\rfloor-s, U=\left\lfloor r_0\varepsilon\right\rfloor+s,\)

其中 \(s\) 是一个小的安全余量，用于覆盖 \(\lfloor\cdot\rfloor\) 截断与误差传播。

直接生成 \(n=N!\) 的全部约数不现实，因此采用经典的meet-in-the-middle思想解决：

对每个素数 \(p\le N\)，勒让德公式给出\(\displaystyle{v_p(N!)=\sum_{k\ge 1}\left\lfloor\frac{N}{p^k}\right\rfloor.}\)于是\(\displaystyle{N!=\prod_{p\le N} p^{v_p(N!)}.}\)

把素数集合分成两组 \(G_1,G_2\)（互不相交），对应指数也分成两组。分别枚举两组生成的所有约数集合 \(D_1,D_2\)，但只保留不超过上界 \(U\) 的值：

  • \(\displaystyle{D_1=\left\{ \prod_{p\in G_1} p^{e_p}\mid 0\le e_p\le v_p(N!), \prod_{p\in G_1} p^{e_p}\le U\right\}}\)；
  • \(\displaystyle{D_2=\left\{ \prod_{p\in G_2} p^{e_p}\mid 0\le e_p\le v_p(N!), \prod_{p\in G_1} p^{e_p}\le U\right\}}\)；

将 \(D_2\) 排序。对于每个 \(x\in D_1\)，要找 \(y\in D_2\) 使得\(L\le xy\le U.\)等价于\(\left\lceil\dfrac{L}{x}\right\rceil \le y \le \left\lfloor\dfrac{U}{x}\right\rfloor.\)

因为 \(D_2\) 已排序，可以对这个区间用二分搜索取得下标范围并批量加入候选：\(C=\{xy\mid x\in D_1,y\in D_2,L\le xy\le U\}.\)

由于 \(G_1,G_2\) 的素因子集合不交，\(xy\) 的分解在这两组之间是唯一的，因此候选不会因为不同乘法路径而重复。

最终，为在候选集合上求最优三元组，我们将 \(C\) 排序。接下来对每个 \(b\in C\)：

  1. 计算 \(q=n/b\)；
  2. 计算 \(a\) 的上界\(A_{\max}=\min\left(b,\left\lfloor\dfrac{q}{b}\right\rfloor\right);\)
  3. 在排序后的 \(C\) 中找到不超过 \(A_{\max}\) 的最大位置，从该位置向下扫描，找到第一个满足 \(a\mid q\) 的 \(a\)；
  4. 令 \(c=q/a\)，此时自动有 \(a\le b\le c\)；
  5. 用交叉乘法比较 \(\dfrac{c}{a}\) 的大小：\(\dfrac{c}{a} < \dfrac{c^\ast}{a^\ast} \iff c\cdot a^\ast < c^\ast\cdot a.\) 若更优则更新答案，并记录最小的 \(a+b+c\)（最优比值唯一时这也唯一）。

代码

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
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
import bisect

from tools import *


def vp_factorial(N: int, p: int) -> int:
    e = 0
    while N:
        N //= p
        e += N
    return e



def icbrt(n: int) -> int:
    """floor(cuberoot(n)) for n>=0 using integer Newton iteration."""
    if n < 0:
        raise ValueError("n must be non-negative")
    if n == 0:
        return 0
    # initial guess: 2^(ceil(bitlen/3))
    x = 1 << ((n.bit_length() + 2) // 3)
    while True:
        y = (2 * x + n // (x * x)) // 3
        if y >= x:
            break
        x = y
    # adjust to be exact floor
    while (x + 1) * (x + 1) * (x + 1) <= n:
        x += 1
    while x * x * x > n:
        x -= 1
    return x

def gen_divisors_limited(primes, exps, limit: int):
    """Generate all divisors formed by primes^exps but keep only <= limit."""
    divs = [1]
    for p, e in zip(primes, exps):
        new = []
        # for each existing divisor d, multiply by p^k while <= limit
        for d in divs:
            base = 1
            for _ in range(e + 1):
                v = d * base
                if v > limit:
                    break
                new.append(v)
                base *= p
        divs = new
    return divs

def choose_split(primes, exps):
    """
    Split prime list into two disjoint groups for meet-in-the-middle.
    Use log of divisor-count product to get a roughly balanced split.
    """
    logs = [math.log(e + 1) for e in exps]
    total = sum(logs)
    target = total / 2.0
    s = 0.0
    split = 0
    for i, lg in enumerate(logs):
        if s + lg > target and i > 0:
            split = i
            break
        s += lg
    if split == 0:
        split = max(1, len(primes) // 2)
    return split

def f_factorisation_triple_of_factorial(N: int, eps_num=10001, eps_den=10000, slack=2000):
    """
    Compute f(N!) by searching factorisation triples (a<=b<=c, abc=N!)
    that minimize c/a, returning (a+b+c, (a,b,c)).

    eps = eps_num/eps_den controls the cube-root window.
    """
    if N < 1:
        raise ValueError("N must be >= 1")

    n = math.factorial(N)

    # integer cube root and window [L,U]
    r = icbrt(n)
    # L = floor(r/eps), U = floor(r*eps), then add slack for safety
    L = max(1, (r * eps_den) // eps_num - slack)
    U = (r * eps_num) // eps_den + slack

    primes = Factorizer(N).primes
    exps = [vp_factorial(N, p) for p in primes]

    # meet-in-the-middle split
    cut = choose_split(primes, exps)
    p1, e1 = primes[:cut], exps[:cut]
    p2, e2 = primes[cut:], exps[cut:]

    # generate divisors limited to U
    D1 = gen_divisors_limited(p1, e1, U)
    D2 = gen_divisors_limited(p2, e2, U)
    D2.sort()

    # build candidates in [L,U]
    cand = []
    for x in D1:
        lo = (L + x - 1) // x
        hi = U // x
        if lo > hi:
            continue
        l = bisect.bisect_left(D2, lo)
        rpos = bisect.bisect_right(D2, hi)
        for y in D2[l:rpos]:
            cand.append(x * y)
    cand.sort()

    best_a = best_b = best_c = None

    # enumerate b in candidates
    for b in cand:
        # b is a divisor of n by construction, but keep the check harmless
        if n % b != 0:
            continue
        q = n // b
        limit = min(b, q // b)  # ensures a<=b and c>=b

        idx = bisect.bisect_right(cand, limit) - 1
        while idx >= 0:
            a = cand[idx]
            if q % a == 0:
                c = q // a
                if c >= b:
                    if best_a is None or c * best_a < best_c * a:
                        best_a, best_b, best_c = a, b, c
                break
            idx -= 1

    return best_a + best_b + best_c, (best_a, best_b, best_c)

if __name__ == "__main__":
    # 示例：43!（应输出 1177163565297340320）
    ans, (a, b, c) = f_factorisation_triple_of_factorial(43)
    print(ans)
    # print(a, b, c)