Project Euler 374

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

Project Euler 374

题目

Maximum Integer Partition Product

An integer partition of a number n is a way of writing \(n\) as a sum of positive integers.

Partitions that differ only in the order of their summands are considered the same. A partition of \(n\) into distinct parts is a partition of \(n\) in which every part occurs at most once.

The partitions of \(5\) into distinct parts are:

\(5, 4+1\) and \(3+2\).

Let \(f(n)\) be the maximum product of the parts of any such partition of \(n\) into distinct parts and let \(m(n)\) be the number of elements of any such partition of \(n\) with that product.

So \(f(5)=6\) and \(m(5)=2\).

For \(n=10\) the partition with the largest product is \(10=2+3+5\), which gives \(f(10)=30\) and \(m(10)=3\).

And their product, \(f(10)\cdot m(10) = 30\cdot3 = 90\)

It can be verified that

\(\sum f(n)\cdot m(n)\) for \(1 \le n \le 100 = 1683550844462\).

Find \(\sum f(n)\cdot m(n)\) for \(1 \le n \le 10^{14}\).

Give your answer modulo \(982451653\), the \(50\) millionth prime.

解决方案

下面用“局部替换增益”的思路刻画最优结构，这些结构在论文中有所介绍：

  1. 不含 1：若拆分里有 \(1\) 和 \(x\) 且 \(x+1\) 不在拆分中，则用 \((x+1)\) 替换 \((1,x)\)，总和不变而乘积变大，因为 \(1\cdot x < x+1\)。因此最大乘积拆分不会包含 \(1\)。
  2. 尽量均匀：若 \(x<y\) 在拆分中，且 \(x+1,y-1\) 都不在拆分中，则用 \((x+1,y-1)\) 替换 \((x,y)\)，乘积增大：\((x+1)(y-1)-xy=y-x-1>0.\)因此最优拆分的数会尽量靠近，呈现“几乎连续”的形态：最多只有一个缺口且缺口大小至多为 \(1\)。
  3. 大数倾向拆成 2 与剩余：若 \(x>4\) 在拆分中且 \(x-2\) 不在拆分中，则用 \((2,x-2)\) 替换 \((x)\)，总和不变而乘积增大，因为 \(2(x-2)=2x-4>x\)。这迫使最优拆分的起点在 \(2\) 或 \(3\) 附近。

综合结构2和3，最优拆分应当是形如 \(2,\dots,k\) 或 \(3,\dots,k\) 的连续段，允许“恰好少一个数”的缺口（且只会发生在非常受限的位置）。

论文的引理2.2指出，令三角数\(T_m=\dbinom{m+1}{2}=\dfrac{m(m+1)}{2}.\)对任意 \(n>1\)，存在唯一的正整数 \(m\) 与 \(l\) 使得\(n=T_m+l, 0\le l\le m,\)且\(m=\left\lfloor\dfrac{\sqrt{8n+1}-1}{2}\right\rfloor, l=n-T_m.\)

论文的定理3.1给出了最大乘积 \(a_n\)（即本题的 \(f(n)\)）的分段公式：若 \(n=T_m+l\)，则 \[ f(n)= \begin{cases} \dfrac{(m+1)!}{m-l}, & 0\le l\le m-2,\\ \dfrac{m+2}{2}m!, & l=m-1,\\ (m+1)!, & l=m. \end{cases} \]

而且证明过程同时给出了达到最大值的拆分形态：

  • 当 \(0\le l\le m-2\) 时，最优拆分为从 \(2\) 到 \((m+1)\) 的连续段缺一个数 \((m-l)\)：\(n=2+3+\cdots+(m-l-1)+(m-l+1)+\cdots+m+(m+1),\)其乘积正是 \(\dfrac{(m+1)!}{m-l}\)。
  • 当 \(l=m-1\) 时，最优拆分为\(n=3+4+\cdots+m+(m+2),\)乘积为 \(\dfrac{m+2}{2}m!\)。
  • 当 \(l=m\) 时，最优拆分为\(n=2+3+\cdots+m+(m+1),\)乘积为 \((m+1)!\)。

由上面三种最优拆分直接数项即可得到：

  • 若 \(l=m\)，拆分为 \(2,3,\dots,m+1\)，长度为 \(m\)；
  • 若 \(0\le l\le m-1\)，拆分长度均为 \(m-1\)（要么缺一个数，要么从 \(3\) 开始但最后加了 \((m+2)\)）。

因此对 \(n\ge 2\)，

\[ m(n)= \begin{cases} m-1, & 0\le l\le m-1,\\ m, & l=m. \end{cases} \] 另外 \(n=1\) 时只有拆分 \(1\)，所以 \(f(1)=1,m(1)=1\)。

接下里正式开始进行计算。固定 \(m\ge 2\)，考虑这一组\(n=T_m+l, l=0,1,\dots,m,\)即从 \(T_m\) 到 \(T_{m+1}-1\) 共 \(m+1\) 个数。

  • 对 \(l=0,\dots,m-2\)：\(f(n)m(n)=(m-1)\cdot \dfrac{(m+1)!}{m-l}.\)令 \(d=m-l\)，则 \(d\) 从 \(m\) 递减到 \(2\)，于是这部分和为\(\displaystyle{(m-1)(m+1)!\sum_{d=2}^{m}\frac{1}{d}.}\)
  • 对 \(l=m-1\)：\(f(n)m(n)=(m-1)\cdot \dfrac{m+2}{2}m!.\)
  • 对 \(l=m\)：\(f(n)m(n)=m\cdot (m+1)!.\)

需要注意的是，最后一组可能不完整。令\(M=\max\{m:T_m\le N\}=\left\lfloor\dfrac{\sqrt{8N+1}-1}{2}\right\rfloor, L=N-T_M.\)

则对 \(m<M\) 的组都是完整的；对 \(m=M\) 只取 \(l=0,\dots,L\)。

当 \(L\le m-2\) 时，只需要计算\(\displaystyle{(m-1)(m+1)!\sum_{d=m-L}^{m}\frac1d,}\) 可见，\(\displaystyle{\sum_{d=m-L}^{m}\frac1d}\)可用前缀和差分实现。

6. Python 参考实现

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
# Project Euler 374
# Sum_{1<=n<=N} f(n)*m(n) mod 982451653
# f(n): maximum product of distinct-part partition
# m(n): number of parts in a partition achieving f(n)

import math
from array import array

MOD = 982451653
INV2 = (MOD + 1) // 2


def max_m_with_tri_leq(N: int) -> int:
    # largest m such that T_m = m(m+1)/2 <= N
    return (math.isqrt(1 + 8 * N) - 1) // 2


def solve(N: int = 10**14) -> int:
    if N <= 0:
        return 0

    # Handle n=1 and n=2 directly
    ans = 0
    if N >= 1:
        ans += 1  # f(1)=1, m(1)=1
    if N >= 2:
        ans += 2  # f(2)=2, m(2)=1
    ans %= MOD
    if N <= 2:
        return ans

    M = max_m_with_tri_leq(N)  # last (possibly partial) triangular-index group

    # We will build inverses inv[k] and prefix sums pref[k]=sum_{i=2..k} inv[i]
    # We need inv up to (M+1) because we will divide by (m+1) to get m! from (m+1)!.
    K = M + 2
    inv = array("I", [0]) * (K + 1)
    pref = array("I", [0]) * (K + 1)

    inv[1] = 1
    fac = 1          # fac = k! in the loop
    pref_sum = 0     # running sum of inv[2..k]

    for k in range(2, M + 2):  # k goes to M+1 inclusive
        # inv[k] via linear recurrence
        invk = (MOD - (MOD // k) * inv[MOD % k] % MOD) % MOD
        inv[k] = invk

        # update prefix sum of inverses
        pref_sum += invk
        if pref_sum >= MOD:
            pref_sum -= MOD
        pref[k] = pref_sum

        # update factorial: fac = k!
        fac = (fac * k) % MOD

        m = k - 1
        if m < 2:
            continue

        if m < M:
            # Full group for this m:
            # G_m = (m-1)(m+1)! * sum_{d=2..m} 1/d
            #     + (m-1)*(m+2)/2 * m!
            #     + m*(m+1)!
            sum_inv_2_to_m = pref[m]  # = sum_{i=2..m} inv[i]
            mfact = fac * invk % MOD  # m! = (m+1)! / (m+1) = k!/k

            term = (m - 1) * fac % MOD
            term = term * sum_inv_2_to_m % MOD
            term = (term + (m - 1) * mfact % MOD * (m + 2) % MOD * INV2) % MOD
            term = (term + m * fac) % MOD

            ans += term
            ans %= MOD
        else:
            # Last group (m == M), possibly partial
            Tm = m * (m + 1) // 2
            L = N - Tm
            if L < 0:
                break
            if L > m:
                L = m

            mfact = fac * invk % MOD  # m!

            if L <= m - 2:
                # only l=0..L terms of type (m+1)!/(m-l)
                lower = m - L  # denominator d ranges lower..m
                sum_inv = (pref[m] - pref[lower - 1]) % MOD
                term = (m - 1) * fac % MOD * sum_inv % MOD
            elif L == m - 1:
                sum_inv_2_to_m = pref[m]
                term = (m - 1) * fac % MOD * sum_inv_2_to_m % MOD
                term = (term + (m - 1) * mfact % MOD * (m + 2) % MOD * INV2) % MOD
            else:
                # L == m: full group
                sum_inv_2_to_m = pref[m]
                term = (m - 1) * fac % MOD * sum_inv_2_to_m % MOD
                term = (term + (m - 1) * mfact % MOD * (m + 2) % MOD * INV2) % MOD
                term = (term + m * fac) % MOD

            ans = (ans + term) % MOD
            break

    return ans


if __name__ == "__main__":
    print(solve(10**14))

7. 本题答案

\[ \sum_{n=1}^{10^{14}} f(n),m(n)\equiv \boxed{334420941}\pmod{982451653}. \]

代码