Project Euler 646

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

Project Euler 646

题目

Bounded Divisors

Let \(n\) be a natural number and \(p_1^{\alpha_1}\cdot p_2^{\alpha_2}\dots p_k^{\alpha_k}\) its prime factorisation.

Define the Liouville function \(\lambda(n)\) as \(\lambda(n) = (-1)^{\sum\limits_{i=1}^{k}\alpha_i}\).

(i.e. \(-1\) if the sum of the exponents \(\alpha_i\) is odd and \(1\) if the sum of the exponents is even. )

Let \(S(n,L,H)\) be the sum \(\lambda(d) \cdot d\) over all divisors \(d\) of \(n\) for which \(L \leq d \leq H\).

You are given:

\(\bullet\, S(10! , 100, 1000) = 1457\)

\(\bullet\, S(15!, 10^3, 10^5) = -107974\)

\(\bullet\, S(30!,10^8, 10^{12}) = 9766732243224\).

Find \(S(70!,10^{20}, 10^{60})\) and give your answer modulo \(1\,000\,000\,007\).

解决方案

令\(N=70,L=10^{20},H=10^{60},M=1000000007\)。

定义\(\gamma(d)=\lambda(d)d.\)若\(d\)的质因子分解为 \(\displaystyle{d=\prod_{i=1}^{k} p_i^{e_i},}\)则\(\displaystyle{\gamma(d)=\left(\prod_{i=1}^{k}p_i^{e_i}\right)\left(\prod_{i=1}^{k}(-1)^{e_i}\right)=\prod_{i=1}^{k}\left((-1)^{e_i}p_i^{e_i}\right)=\prod_{i=1}^{k}(-p_i)^{e_i}.}\)因此 \(\gamma\) 是完全乘法性的（更准确：对互素的数乘法性足够）：若 \(\gcd(a,b)=1\)，则\(\gamma(ab)=\gamma(a)\gamma(b).\)

将区间求和写成前缀差：令\(\displaystyle{F(T)=\sum_{\substack{d\mid N!\\d\le T}}\gamma(d),}\)则\(S(N!,L,H)=F(H)-F(L-1).\) 在本题中端点 \(L,H\) 不是 \(N!\) 的因子时，\(F(L-1)=F(L)\)，

于是可用\(S(N!,L,H)=F(H)-F(L).\)实际实现里直接用对数阈值 \(\log H\) 与 \(\log L\) 即可。

对每个素数 \(p\le N\)，令 \(\beta_p=v_p(N!)\)，由 Legendre 公式得 \(\displaystyle{\beta_p=\sum_{t\ge 1}\left\lfloor\frac{N}{p^t}\right\rfloor.}\)因此\(\displaystyle{N! = \prod_{p\le N} p^{\beta_p}.}\) 每个因子 \(d\mid N!\) 唯一对应一组指数选择 \(0\le e_p\le \beta_p\)，并且\(\displaystyle{d=\prod_{p\le N} p^{e_p}, \gamma(d)=\prod_{p\le N}(-p)^{e_p}.}\)

接下来考虑 Meet-in-the-middle思想解决。选取一个划分，将素数集合拆成不交两部分 \(A,B\)，定义\(\displaystyle{P=\prod_{p\in A}p^{\beta_p}, Q=\prod_{p\in B}p^{\beta_p},}\)则 \(N!=PQ\) 且 \(\gcd(P,Q)=1\)。

任意 \(d\mid N!\) 有唯一表示\(d=xy, x\mid P,y\mid Q,\)并且由乘法性\(\gamma(d)=\gamma(xy)=\gamma(x)\gamma(y).\)

于是前缀和\(F(T)\)可写为 \[ F(T)=\sum_{\substack{x\mid P}}\gamma(x)\sum_{\substack{y\mid Q\\xy\le T}}\gamma(y). \]

直接比较 \(xy\le T\) 会遇到 \(T\) 高达 \(H\) 的大整数。改用对数：\(xy\le T\Longleftrightarrow \log x+\log y\le \log T.\)

做法：

  1. 枚举所有 \(y\mid Q\)，保存二元组 \((\log y ,\gamma(y))\)。
  2. 按 \(\log y\) 升序排序，做前缀和数组\(\displaystyle{Y_{s}[k]=\sum_{i=0}^k \gamma(y_i).}\)
  3. 对每个 \(x\mid P\)，计算阈值\(\log y\le \log T-\log x,\)在排序后的 \(\log y\) 上二分出最大下标 \(k\)，贡献为\(\gamma(x)\cdot Y_{s}[k].\)

当然，我们并不需要生成超过 \(T\) 的数：因为 \(x\ge 1,y\ge 1\)，若 \(x>T\) 则任何组合都不可能满足 \(xy\le T\)。因此在生成每一侧因子时，可以用对数剪枝：当把已有项乘上一次 \(p\)，对数增加 \(\log p\)；若超过 \(\log T\) 就停止继续扩展该条链。此外，\(\gamma\) 的更新是把当前模值乘上 \((-p)\bmod M\)。

代码

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
from math import log
from bisect import bisect_right
from tools import Factorizer

MOD = 1000000007
N = 70
L = 10 ** 20
R = 10 ** 60


def v_p(n, p):
    s = 0
    while n:
        n //= p
        s += n
    return s


def prefix_sum_bound(n, max_log):
    ps = Factorizer(n).primes
    es = [v_p(n, p) for p in ps]
    lps = [log(p) for p in ps]

    L = [(1, 0.0)]
    R = [(1, 0.0)]
    l = 0
    r = len(ps) - 1

    while l <= r:
        if len(L) < len(R):
            divs = L
            idx = l
            l += 1
        else:
            divs = R
            idx = r
            r -= 1

        p = ps[idx]
        e = es[idx]
        lp = lps[idx]
        negp = MOD - (p % MOD)

        old = len(divs)
        for t in range(old):
            v0, lg0 = divs[t]
            vv = v0
            llg = lg0
            for _ in range(e):
                llg += lp
                if llg > max_log:
                    break
                vv = (vv * negp) % MOD
                divs.append((vv, llg))

    L.sort(key=lambda x: x[1])
    logs = [lg for _, lg in L]
    pref = [0] * len(L)
    s = 0
    for i, (v, _) in enumerate(L):
        s += v
        s %= MOD
        pref[i] = s

    ans = 0
    for v, lg in R:
        target = max_log - lg
        k = bisect_right(logs, target) - 1
        if k >= 0:
            ans = (ans + v * pref[k]) % MOD
    return ans


def solve(N, L, R):
    r = log(R)
    l = log(L - 1)
    a = prefix_sum_bound(N, r)
    b = prefix_sum_bound(N, l)
    return (a - b) % MOD


if __name__ == "__main__":
    print(solve(N,L,R))