Project Euler 834

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

Project Euler 834

题目

Add and Divide

A sequence is created by starting with a positive integer \(n\) and incrementing by \((n+m)\) at the \(m^{th}\) step. If \(n=10\), the resulting sequence will be \(21,33,46,60,75,91,108,126,\ldots\).

Let \(S(n)\) be the set of indices \(m\), for which the \(m^{th}\) term in the sequence is divisible by \((n+m)\).

For example, \(S(10)=\{5,8,20,35,80\}\).

Define \(T(n)\) to be the sum of the indices in \(S(n)\). For example, \(T(10) = 148\) and \(T(10^2)=21828\).

Let \(\displaystyle U(N)=\sum_{n=3}^{N}T(n)\). You are given, \(U(10^2)=612572\).

Find \(U(1234567)\).

解决方案

从 \(n\) 开始，第 \(1\) 次加 \(n+1\)，第 \(2\) 次加 \(n+2\)，……第 \(m\) 次加 \(n+m\)。因此第 \(m\) 项为\(\displaystyle{A_m = n+\sum_{t=1}^{m}(n+t)= n + mn + \frac{m(m+1)}{2}= n(m+1)+\frac{m(m+1)}{2}.}\)整理得更紧凑的形式：

\[ A_m=\frac{(m+1)(2n+m)}{2}. \]

题目要求的是：哪些 \(m\) 满足\((n+m)\mid A_m.\)

令\(d=n+m (m=d-n).\)如果 \((n+m)\mid A_m\)，则存在整数 \(k\ge 1\) 使得\(A_m = k(n+m)=kd.\)把 \(A_m\) 的表达式代入，并乘以 \(2\) 避免分数：\(2A_m=(m+1)(2n+m).\)

用 \(m=d-n\) 替换：\(2A_m=(d-n+1)(d+n).\)展开得到一个非常关键的恒等式：\((d-n+1)(d+n)=d^2+d-n(n-1).\)

于是\(2kd=d^2+d-n(n-1),\)移项并提取公因子 \(d\)：\(d(d+1-2k)=n(n-1).\)

也就是说，若 \(m\) 是可行下标，则存在某个整数 \(k\) 满足

\[ (n+m)\bigl((n+m)+1-2k\bigr)=n(n-1). \]

上式立即推出：\(d=n+m \mid n(n-1).\)

并且由于 \(m\ge 1\)，必有\(d=n+m>n.\)因此每个可行的 \(m\) 都对应某个满足\(d\mid n(n-1), d>n\)的因子 \(d\)，且 \(m=d-n\)。但这还不够：还需要保证 \(k\) 是整数。

从\(d(d+1-2k)=n(n-1)\)得到\(d+1-2k=\dfrac{n(n-1)}{d}.\)因此\(2k=d+1-\dfrac{n(n-1)}{d}.\)要让 \(k\) 为整数，必须有\(d+1 \equiv \dfrac{n(n-1)}{d}\pmod 2.\)

注意到 \(d+1\) 与 \(d\) 奇偶相反，所以这等价于：\(d \not\equiv \frac{n(n-1)}{d}\pmod 2,\)也就是 \(d\) 与 \(\dfrac{n(n-1)}{d}\) 必须一奇一偶。

再看 \(n(n-1)\) 的 \(2\)-adic 结构。由于 \(n\) 与 \(n-1\) 互素，且其中恰好一个是偶数，因此设\(n(n-1)=2^{t}\cdot ab,\)其中\(a,b\)都是奇数。其中\(t\)满足

\[ t=v_2(n(n-1))= \begin{cases} v_2(n), & n\equiv 0\pmod 2,\\ v_2(n-1), & n\equiv 1\pmod 2. \end{cases} \]

若因子 \(d\) 含有 \(2^s\)（即 \(v_2(d)=s\)），则商 \(\dfrac{n(n-1)}{d}\) 含有 \(2^{t-s}\)。要让二者一奇一偶，必须满足：\(s=0\)，则 \(d\) 为奇数，商含 \(2^t\) 为偶数；或 \(s=t\)，则 \(d\) 含满全部 \(2^t\) 为偶数，商为奇数；而 任何 \(0<s<t\) 都会让二者同时为偶数，从而不合法。

因此可行的 \(d\) 必须满足：\(v_2(d)\in\{0,t\}.\)也就是说：\(d\) 要么是 \(n(n-1)\) 的 奇因子，要么是“带上全部 \(2^t\)”的 合格偶因子。

这时可以写成\(a=\dfrac{n}{2^{v_2(n)}},b=\dfrac{n-1}{2^{v_2(n-1)}}.\)由于 \(n\) 与 \(n-1\) 互素，故 \(a\) 与 \(b\) 也互素，且均为奇数。那么 \(n(n-1)\) 的所有奇因子恰好是\(\{d_1d_2: d_1\mid a,d_2\mid b\}.\)设一个奇因子为\(p=d_1d_2.\)那么所有可行的 \(d\) 只有两类：

  1. 奇因子本身：\(d=p\)
  2. 带满全部 \(2^t\) 的偶因子：\(d=2^t p\)

并且需要满足 \(d>n\)，从而对应的 \(m=d-n>0\)。因此，对奇因子部分：需要 \(p>n\)；对合格偶因子部分：需要 \(2^t p>n\)，等价于 \(p>\left\lfloor \dfrac{n}{2^t}\right\rfloor\)。

每个合法的 \(m\) 都是 \(m=d-n\)，因此\(\displaystyle{T(n)=\sum_{d\in D_n,d>n}(d-n)=\left(\sum_{d\in D_n,d>n} d\right)-n\cdot\#\{d\in D_n:d>n\}.}\)其中 \(D_n\) 表示所有“可行”的 \(d\) 的集合（也就是满足条件的因子集合）。

把 \(d\) 用 \(p=d_1d_2\) 表示后，我们需要两次统计：

  • 阈值 \(L=n\) 的奇因子对：\(\displaystyle{\sum_{d_1\mid a}\sum_{d_2\mid b,d_1d_2>L}(d_1d_2-n)}\)
  • 阈值 \(L=\left\lfloor \dfrac{n}{2^t}\right\rfloor\) 的合格偶因子对：\(\displaystyle{\sum_{d_1\mid a}\sum_{d_2\mid b,d_1d_2>L}(2^t d_1d_2-n)}\)

于是只要能高效地求出对于给定阈值 \(L\)，所有满足 \(d_1d_2>L\) 的乘积之和与个数，就能得到 \(T(n)\)，进而求出 \(U(N)\)。

将 \(a\) 的因子表记为升序数组 \(A\)，\(b\) 的因子表记为升序数组 \(B\)。对固定阈值 \(L\)，我们希望统计所有满足\(x\in A,y\in B,xy>L\)的贡献。因此，使用双指针法我们可以完成这个过程。

最终预处理每个奇数的因子后，我们可以通过上述的双指针法求出\(T(n)\)即可。

代码

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
N = 1234567

def solve(N):
    M = N // 2 + 1
    div = [[] for _ in range(M)]
    for d in range(1, N + 1, 2):
        step = d << 1
        for x in range(d, N + 1, step):
            div[x >> 1].append(d)

    def tz(x):
        return (x & -x).bit_length() - 1

    ans = 0

    for n in range(3, N + 1):
        if n & 1 == 0:
            t = tz(n)
            a = n >> t
            b = n - 1
        else:
            x = n - 1
            t = tz(x)
            a = n
            b = x >> t

        A = div[a >> 1]
        B = div[b >> 1]

        j = len(B) - 1
        cnt = 0
        s = 0
        for da in A:
            while j >= 0 and da * B[j] > n:
                s += B[j]
                cnt += 1
                j -= 1
            ans += da * s - n * cnt

        limit = n >> t
        j = len(B) - 1
        cnt = 0
        s = 0
        for da in A:
            while j >= 0 and da * B[j] > limit:
                s += B[j]
                cnt += 1
                j -= 1
            ans += (da * s) << t
            ans -= n * cnt

    return ans

print(solve(N)