Project Euler 771

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

Project Euler 771

题目

Pseudo Geometric Sequences

We define a pseudo-geometric sequence to be a finite sequence \(a_0, a_1, \dotsc, a_n\) of positive integers, satisfying the following conditions:

  • \(n \geq 4\), i.e. the sequence has at least 5 terms.

  • \(0 < a_0 < a_1 < \dotsc < a_n\), i.e. the sequence is strictly increasing.

  • \(| a_i^2 - a_{i - 1}a_{i + 1} | \le 2\) for \(1 \le i \le n-1\).

Let \(G(N)\) be the number of different pseudo-geometric sequences whose terms do not exceed \(N\).

For example, \(G(6) = 4\), as the following \(4\) sequences give a complete list:

\[1, 2, 3, 4, 5 \qquad 1, 2, 3, 4, 6 \qquad 2, 3, 4, 5, 6 \qquad 1, 2, 3, 4, 5, 6\]

Also, \(G(10) = 26, G(100) = 4710\) and \(G(1000) = 496805\).

Find \(G(10^{18})\). Give your answer modulo \(1\,000\,000\,007\).

解决方案

对每个 \(i\) 定义整数\(c_i=a_i^2-a_{i-1}a_{i+1}, c_i\in\{-2,-1,0,1,2\}.\)于是\(a_{i+1}=\dfrac{a_i^2-c_i}{a_{i-1}}.\)因此一旦给定相邻两项 \((a_{i-1},a_i)\)，下一项只能从 \(5\) 个候选 \(c_i\) 里尝试得到。

固定正整数 \(x<y\)，考虑满足\(|y^2-xz|\le 2\)的整数 \(z\)。等价于\(\dfrac{y^2-2}{x}\le z\le \dfrac{y^2+2}{x},\)该区间长度为\(\dfrac{4}{x}.\)若 \(x\ge 5\)，则区间长度 \(<1\)，因此最多只有一个整数落在其中。

由此我们可以导出两条结论：

  • 结论 1（唯一延伸）：当某一步的前一项 \(a_{i-1}\ge 5\) 时，若存在合法的 \(a_{i+1}\)，它是唯一的（不会分叉）。
  • 结论 2（唯一回溯）：同理把序列反向看，也得到，当 \(a_{i+1}\ge 5\) 时，若存在合法的 \(a_{i-1}\)，它也唯一。

这意味着：一条非平凡的伪等比数列，一旦进入数值不小的区间，整体形态几乎被锁死；真正会分叉的选择只可能发生在非常开头的小数值处。

接下来我们介绍两类解。

公差为\(1\)的等差数列

若 \(a_i=i+t\)（连续整数），则\(a_i^2-a_{i-1}a_{i+1}=(t+i)^2-(t+i-1)(t+i+1)=1,\)显然满足要求。

我们有更强的引理（连续性强制传播）：若某处出现三连项 \((x,x+1,x+2)\)，则之后的项被强制为 \(x+3,x+4,\dots\)，即序列从此只能继续连续整数。

证明要点：由唯一延伸结论，当分母够大时只能得到唯一整数解，而 \(x+3\) 总在允许区间内且必须 \(>x+2\)，所以被强制。

因此长链条的主体来自连续整数序列，我们可以直接计数。长度至少 \(5\) 等价于终点与起点差至少 \(4\)。设起点 \(s\)，终点 \(t\)，则\(1\le s<t\le N, t-s\ge 4.\)计数为

\[G_{1}(N)=\sum_{s=1}^{N-4}(N-(s+4)+1)=\frac{(N-4)(N-3)}{2}.\]

严格等比数列

若处处 \(c_i=0\)，则满足\(a_{i+1}a_{i-1}=a_i^2,\)这正是相邻三项成等比的条件，因此序列是严格等比数列。

严格递增的整数等比数列允许公比是既约分数 \(\dfrac{k}{\ell}>1\)（\(\gcd(k,\ell)=1\)）。长度为 \(m+1\)（\(m\ge 4\)）时，所有项必须为整数，可写作\((a\ell^{m},a k\ell^{m-1},a k^{2}\ell^{m-2},\dots,a k^{m}),\)其中 \(k>\ell\ge 1\)，且 \(a\ge 1\)。

因为 \(\ell<k\)，最大项是最后一项 \(a k^m\)，只需\(a k^{m}\le N\Longleftrightarrow a\le \left\lfloor\dfrac{N}{k^m}\right\rfloor.\)

对固定 \(k\)，满足 \(1\le \ell<k\)、\(\gcd(k,\ell)=1\) 的 \(\ell\) 的个数是欧拉函数 \(\varphi(k)\)。于是

\[G_{2}(N)=\sum_{k\ge 2}\varphi(k)\sum_{m\ge 4}\left\lfloor\frac{N}{k^m}\right\rfloor.\]

注意当 \(k^4>N\) 时内层全为 \(0\)，所以只需到\(k\le \lfloor \sqrt[4]{N}\rfloor.\)

剩余零散序列

将非连续、非纯等比的序列统称为（零散型）。这部分数量相对很少，原因是：始阶段可能有最多 \(5\) 个分支（来自 \(c_i\in[-2,2]\)）。一旦分母 \(\ge 5\)；唯一性保证几乎不产生分叉；此外，序列增长很快，长度最多约 \(O(\log N)\)。

由\(a_{i+1}=\dfrac{a_i^2-c}{a_{i-1}}, c\in\{-2,-1,0,1,2\},\)枚举所有可行的 \(c\)，要求：

  • 能整除：\((a_i^2-c)\bmod a_{i-1}=0\)；
  • 严格递增：\(a_{i+1}>a_i\)；
  • 不超过上界：\(a_{i+1}\le N\)。

从 \(a_0=1\) 出发，第二项 \(a_1=b\) 后，第三项必在 \(b^2\pm 2\) 附近。继续递推时整体大致按“倍乘”增长，长度达到 \(5\) 时末项规模约为 \(b^4\)。因此若 \(b>N^{1/4}+O(1)\)，长度 \(\ge 5\) 的序列无法保持全体 \(\le N\)。

因此只需枚举\(2\le a_1\le \lfloor N^{1/4}\rfloor+5\)即可覆盖所有可能（多出的常数只是为了稳妥）。

我们最终要数的是所有起点、所有长度\(\ge 5\)的序列。当我们递归构造出一个序列前缀\(A=(a_0,a_1,\dots,a_{L-1})\)且 \(L\ge 5\) 时，所有以 \(a_{L-1}\) 结尾、长度\(\ge 5\) 的连续子段恰好是\((a_s,a_{s+1},\dots,a_{L-1}), s=0,1,\dots,L-5.\)这些子段自动满足伪等比条件（因为条件只依赖局部三元组）。

但我们不希望把连续整数序列或纯等比数列重复计入（它们已计入 \(G_{1},G_{2}\) ）。因此对每个后缀子段，我们仅当它 既不是连续整数、也 不是等比数列 时才计入 \(G_3\) 贡献。

于是最终得到：

\[G(N)=G_{1}(N)+G_{2}(N)+G_{3}(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
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
from tools import *

MOD = 1_000_000_007
N = 10 ** 18

def iroot4(n: int) -> int:
    return int_sqrt(int_sqrt(n))

def phi_sieve(n: int):
    phi = list(range(n + 1))
    for i in range(2, n + 1):
        if phi[i] == i:
            for j in range(i, n + 1, i):
                phi[j] -= phi[j] // i
    return phi

def is_geo(seq: List[int]) -> bool:
    for i in range(len(seq) - 2):
        if seq[i] * seq[i + 2] != seq[i + 1] * seq[i + 1]:
            return False
    return True

def is_consecutive(seq: List[int]) -> bool:
    for i in range(1, len(seq)):
        if seq[i] != seq[i - 1] + 1:
            return False
    return True

def count_suffix_spor(seq: List[int]) -> int:
    L = len(seq)
    res = 0
    for s in range(0, L - 4):
        sub = seq[s:]
        if is_consecutive(sub):
            continue
        if is_geo(sub):
            continue
        res += 1
    return res

def G1(N: int) -> int:
    if N < 5:
        return 0
    return (N - 4) * (N - 3) * pow(2, MOD - 2, MOD) % MOD

def G2(N: int, phi: List[int], K: int) -> int:
    ans = 0
    for k in range(2, K + 1):
        pk = k * k * k * k
        if pk > N:
            break
        while pk <= N:
            ans = (ans + phi[k] * (N // pk)) % MOD
            if pk > N // k:
                break
            pk *= k
    return ans

def G3(N: int, K: int) -> int:
    ans = 0
    def dfs(seq: List[int]):
        nonlocal ans
        L = len(seq)
        if L >= 5:
            ans = (ans + count_suffix_spor(seq)) % MOD
        if L >= 5 and is_consecutive(seq):
            return
        if L >= 5 and is_geo(seq):
            return
        a0, a1 = seq[-2], seq[-1]
        for c in (-2, -1, 0, 1, 2):
            num = a1 * a1 - c
            if num % a0:
                continue
            a2 = num // a0
            if a2 <= a1 or a2 > N:
                continue
            seq.append(a2)
            dfs(seq)
            seq.pop()
    for a1 in range(2, K + 6):
        dfs([1, a1])
    return ans

def G(N: int) -> int:
    K = iroot4(N)
    phi = phi_sieve(K)
    ans = (G1(N) + G2(N, phi, K) + G3(N, K)) % MOD
    return ans

print(G(N))