Project Euler 738

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

Project Euler 738

题目

Counting Ordered Factorisations

Define \(d(n,k)\) to be the number of ways to write \(n\) as a product of \(k\) ordered integers

\[n = x_1\times x_2\times x_3\times \ldots\times x_k \qquad 1\le x_1\le x_2\le\ldots\le x_k\]

Further define \(D(N,K)\) to be the sum of \(d(n,k)\) for \(1\le n\le N\) and \(1\le k\le K\).

You are given that \(D(10, 10) = 153\) and \(D(100, 100) = 35384\).

Find \(D(10^{10},10^{10})\) giving your answer modulo \(1\,000\,000\,007\).

解决方案

令\(N=10^{10},K=10^{10}\)。题目定义\(d(n,k)=\#\left\{(x_1,\dots,x_k)\mid x_1x_2\cdots x_k=n,1\le x_1\le\cdots\le x_k\right\}.\)并令\(\displaystyle{D(N,K)=\sum_{n=1}^{N}\sum_{k=1}^{K} d(n,k).}\)

对固定的 \(k\) 考虑内层求和：\(\displaystyle{E(N,k)=\sum_{n=1}^{N} d(n,k).}\)

注意到：每一个满足 \(x_1\le\cdots\le x_k\) 的 \(k\)-元组都有唯一的乘积 \(n=x_1\cdots x_k\)，而它会在 \(\sum_{n\le N} d(n,k)\) 中被恰好计数一次，当且仅当该乘积 \(\le N\)。

因此\(E(N,k)=\#\left\{(x_1,\dots,x_k)\mid 1\le x_1\le\cdots\le x_k,x_1x_2\cdots x_k\le N\right\}.\)

于是题目目标变为\(\displaystyle{D(N,K)=\sum_{k=1}^{K} E(N,k).}\)

为了处理非降序的约束，定义带下界的计数函数：\(F(N,k,L)=\#\left\{(x_1,\dots,x_k)\mid L\le x_1\le\cdots\le x_k,x_1\cdots x_k\le N\right\}.\)显然我们要的就是\(E(N,k)=F(N,k,1).\)

当 \(k=1\) 时，只剩一个数 \(x_1\)，条件为 \(L\le x_1\le N\)，所以\(F(N,1,L)=\max(N-L+1,0).\)

当 \(k\ge2\) 时，固定第一个数 \(x_1=a\)。由于非降序，剩余部分必须满足\(a\le x_2\le\cdots\le x_k, x_2\cdots x_k\le\left\lfloor \dfrac{N}{a}\right\rfloor.\)因此对应的方案数就是\(F\left(\left\lfloor\dfrac{N}{a}\right\rfloor,k-1,a\right).\)

把 \(a\) 从最小允许值 \(L\) 枚举即可：\(\displaystyle{F(N,k,L)=\sum_{a=L}^{\infty} F\left(\left\lfloor\frac{N}{a}\right\rfloor,k-1,a\right).}\)

当然并不需要枚举到无穷。因为所有因子都满足 \(x_i\ge a\)，若存在解则必须有\(a^k\le x_1x_2\cdots x_k \le N,\)因此 \(a\le \lfloor N^{1/k}\rfloor\)。所以最终递推为

\[ F(N,k,L)=\sum_{a=L}^{\lfloor N^{1/k}\rfloor}F\left(\left\lfloor\frac{N}{a}\right\rfloor,k-1,a\right) \]

并有剪枝\(L^k>N \Rightarrow F(N,k,L)=0.\)

令\(k_{\max}=\left\lfloor \log_2 N \right\rfloor.\)对任意满足 \(x_1\le\cdots\le x_k\) 且乘积 \(\le N\) 的元组，如果所有 \(x_i\ge 2\)，则乘积至少为 \(2^k\)。当 \(k\ge k_{\max}+1\) 时有 \(2^k>N\)，所以不可能所有因子都 \(\ge2\)，从而必然存在因子等于 \(1\)。

又因为序列是非降序，若出现 \(1\)，它只能出现在最前面；因此对 \(k\ge k_{\max}+1\)，每个可行元组都满足 \(x_1=1\)。于是我们得到一个非常关键的双射：

  • 从长度 \(k\) 的可行元组去掉第一个 \(1\)，得到长度 \(k-1\) 的可行元组；
  • 反过来，在任意长度 \(k-1\) 的可行元组前面加一个 \(1\)，得到长度 \(k\) 的可行元组。

因此对所有 \(k\ge k_{\max}+1\) 有\(E(N,k)=E(N,k-1),\)从而\(E(N,k)=E(N,k_{\max}),\forall k\ge k_{\max}.\)

于是当 \(K\ge k_{\max}\) 时：

\[ D(N,K)=\sum_{k=1}^{k_{\max}-1}E(N,k)+(K-k_{\max}+1)\cdot E(N,k_{\max}) \]

因此在本题中，只要单独算出 \(E(N,1),E(N,2),\dots,E(N,k_{\max})\) 就够了，之后的 \(k\) 全部相同。

值得注意的是，我们会单独优化 \(k=2\)：\(F\displaystyle{(N,2,L)=\sum_{a=L}^{\lfloor\sqrt N\rfloor}(\lfloor N/a\rfloor-a+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
from tools import int_sqrt
from functools import cache
import gmpy2

N = 10**10
K = 10 ** 10
MOD = 1_000_000_007

def iroot(n: int, k: int) -> int:
    return int(gmpy2.iroot(n, k)[0])

@cache
def F(N: int, k: int, L: int) -> int:
    if k == 1:
        if N < L:
            return 0
        return (N - L + 1) % MOD

    if k == 2:
        s = int_sqrt(N)
        if L > s:
            return 0
        qsum = 0
        for a in range(L, s + 1):
            qsum += N // a
        sz = s - L + 1
        a_sum = (L + s) * sz // 2
        return (qsum + sz - a_sum) % MOD

    r = iroot(N, k)
    if r < L:
        return 0

    s = 0
    for a in range(L, r + 1):
        s += F(N // a, k - 1, a)
        s %= MOD

    return s


def D(N: int, K: int) -> int:
    r = N.bit_length() - 1
    ans = 0
    for k in range(1, r + 1):
        ans = (ans + F(N, k, 1)) % MOD
    ans = (ans + ((K - r) % MOD) * F(N, r, 1)) % MOD
    return ans


print(D(N, K))