Project Euler 350

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

Project Euler 350

题目

Constraining the least greatest and the greatest least

A list of size \(n\) is a sequence of \(n\) natural numbers.

Examples are \((2,4,6), (2,6,4), (10,6,15,6),\) and \((11)\).

The greatest common divisor, or \(\gcd\), of a list is the largest natural number that divides all entries of the list.

Examples: \(\gcd(2,6,4) = 2, \gcd(10,6,15,6) = 1\) and \(\gcd(11) = 11\).

The least common multiple, or \(\text{lcm}\), of a list is the smallest natural number divisible by each entry of the list.

Examples: \(\text{lcm}(2,6,4) = 12, \text{lcm}(10,6,15,6) = 30\) and \(\text{lcm}(11) = 11\).

Let \(f(G, L, N)\) be the number of lists of size \(N\) with \(\gcd \ge G\) and \(\text{lcm} \le L\). For example:

\(\begin{aligned} &f(10, 100, 1) = 91.\\ &f(10, 100, 2) = 327.\\ &f(10, 100, 3) = 1135.\\ &f(10, 100, 1000) \bmod 101^4 = 3286053. \end{aligned}\)

Find \(f(10^6, 10^{12}, 10^{18}) \bmod 101^4\).

解决方案

Project Euler 350 题解

可见，\(f(G,L,N)=\#\{(a_1,\dots,a_N)\in\mathbb{N}^N:\gcd(a_1,\dots,a_N)\ge G,\mathrm{lcm}(a_1,\dots,a_N)\le L\}.\)

令\(d=\gcd(a_1,\dots,a_N), d\ge G.\)写\(a_i=db_i.\)则\(\gcd(b_1,\dots,b_N)=1,\)并且 \(\mathrm{lcm}(a_1,\dots,a_N)=d\cdot \mathrm{lcm}(b_1,\dots,b_N)\le L.\)

再令\(k=\mathrm{lcm}(b_1,\dots,b_N).\)那么对固定的 \(k\)，允许的 \(d\) 恰是\(G\le d\le \left\lfloor\dfrac{L}{k}\right\rfloor,\)因此可选的 \(d\) 个数为\(\left\lfloor\dfrac{L}{k}\right\rfloor-G+1,\)前提是 \(\left\lfloor\dfrac{L}{k}\right\rfloor\ge G\)，也即 \(k\le \dfrac{L}{G}.\)

于是得到核心分解公式： \[f(G,L,N)=\sum_{k=1}^{\lfloor L/G\rfloor} F(k,N)\left(\left\lfloor\frac{L}{k}\right\rfloor-G+1\right),\]

其中\(F(k,N)=\#\{(b_1,\dots,b_N):\gcd(b_1,\dots,b_N)=1,\mathrm{lcm}(b_1,\dots,b_N)=k\}.\)

所以问题变成：对所有 \(1\le k\le M\)（其中 \(M=\lfloor L/G\rfloor\)），快速计算 \(F(k,N)\bmod P\)。在本题参数下，\(M=\dfrac{L}{G}=10^{6}\)

把\(k\)分解质因数：\(\displaystyle{k=\prod_{i=1}^s p_i^{e_i}}\)。由于每个 \(b_j\) 必须整除 \(k\)，可写为 \(\displaystyle{b_j=\prod_{i=1}^s p_i^{\alpha_{ij}}, 0\le \alpha_{ij}\le e_i.}\)

而

  • \(\gcd(b_1,\dots,b_N)=1\) 等价于对每个 \(p_i\)，\(\displaystyle{\min_j \{\alpha_{ij}\}=0;}\)
  • \(\mathrm{lcm}(b_1,\dots,b_N)=k\) 等价于对每个 \(p_i\)，\(\displaystyle{\max_j\{\alpha_{ij}\}=e_i.}\)

这些条件对不同素数 \(p_i\) 完全独立，因此\(\displaystyle{F(k,N)=\prod_{i=1}^s H(e_i,N),}\)其中 \(H(e,N)\) 是长度为 \(N\) 的序列 \((\alpha_1,\dots,\alpha_N)\)，满足\(\displaystyle{0\le \alpha_j\le e, \min_j\{\alpha_j\}=0, \max_j\{\alpha_j\}=e}\) 的方案数。

接下来求求 \(H(e,N)\) 的闭式解。可见，总共有 \((e+1)^N\) 个序列。

设事件

  • \(A\)：所有 \(\alpha_j\ge 1\)（即最小值不为 \(0\)），则 \(|A|=e^N\)；
  • \(B\)：所有 \(\alpha_j\le e-1\)（即最大值不为 \(e\)），则 \(|B|=e^N\)；
  • \(A\cap B\)：\(1\le \alpha_j\le e-1\)，则 \(|A\cap B|=(e-1)^N\)。

所以

\[H(e,N)=(e+1)^N-|A\cup B|=(e+1)^N-2e^N+(e-1)^N.\]

于是有 \[F(k,N)=\prod_{i=1}^s((e_i+1)^N-2(e_i)^N+(e_i-1)^N).\]

注意这里 \(e_i\ge 1\)；当 \(e_i=1\) 时 \((e_i-1)^N=0^N=0\)（因为 \(N>0\)）。

代码

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
from typing import List
from tools import Factorizer

MOD = 101 ** 4
G = 10 ** 6
L = 10 ** 12
N = 10 ** 18

def solve(G: int, L: int, N: int, mod: int) -> int:
    M = L // G
    spf = Factorizer(M).spf

    # exponent upper bound (safe)
    max_e = 0
    t = M
    while t:
        t //= 2
        max_e += 1
    max_e += 2  # small cushion

    powN = [pow(a, N, mod) for a in range(max_e + 3)]

    # delta[e] = (e+1)^N - 2*e^N + (e-1)^N  (e>=1)
    delta = [0] * (max_e + 3)
    for e in range(1, max_e + 2):
        delta[e] = (powN[e + 1] - 2 * powN[e] + powN[e - 1]) % mod

    ans = 0
    for k in range(1, M + 1):
        x = k
        cur = 1
        while x > 1:
            p = spf[x]
            e = 0
            while x % p == 0:
                x //= p
                e += 1
            cur = (cur * delta[e]) % mod

        weight = (L // k) - G + 1
        ans = (ans + (weight % mod) * cur) % mod

    return ans

print(solve(G, L, N, MOD))