Project Euler 569

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

Project Euler 569

题目

Prime Mountain Range

A mountain range consists of a line of mountains with slopes of exactly \(45°\), and heights governed by the prime numbers, \(p_n\). The up-slope of the \(k^{\text{th}}\) mountain is of height \(p_{2k-1}\), and the downslope is \(p_{2k}\). The first few foot-hills of this range are illustrated below.

Tenzing sets out to climb each one in turn, starting from the lowest. At the top of each peak, he looks back and counts how many of the previous peaks he can see. In the example above, the eye-line from the third mountain is drawn in red, showing that he can only see the peak of the second mountain from this viewpoint. Similarly, from the \(9^{\text{th}}\) mountain, he can see three peaks, those of the \(5^{\text{th}},7^{\text{th}}\) and \(8^{\text{th}}\) mountain.

Let \(P(k)\) be the number of peaks that are visible looking back from the \(k^{\text{th}}\) mountain. Hence \(P(3)=1\) and \(P(9)=3\).

Also \(\displaystyle \sum_{k=1}^{100} P(k) = 227\).

Find \(\displaystyle \sum_{k=1}^{2500000} P(k)\).

解决方案

令\(K=2500000\)。设素数序列为 \(p_1=2,p_2=3,p_3=5,\dots\)。第 \(k\) 座山上坡高度是 \(p_{2k-1}\)，下坡高度是 \(p_{2k}\)，且坡度恒为 \(45^\circ\)，因此每段斜坡的水平位移等于竖直位移。

把每座山的峰顶记为点 \(P_k=(x_k,y_k)\)。从第 \(k-1\) 个峰到第 \(k\) 个峰，会先下坡 \(p_{2k-2}\)，再上坡 \(p_{2k-1}\)，于是：\(\Delta x_k = p_{2k-2}+p_{2k-1},\Delta y_k = -p_{2k-2}+p_{2k-1}.\)

并且第一个峰是从起点上坡 \(p_1\) 得到：\(x_1=p_1, y_1=p_1.\)因此递推为：

\[ x_k = x_{k-1}+p_{2k-2}+p_{2k-1},y_k = y_{k-1}-p_{2k-2}+p_{2k-1},(k\ge 2). \]

特别地，因为 \(p_{2k-1}>p_{2k-2}\)，所以 \(\Delta y_k>0\)，从而\(y_1<y_2<\cdots<y_K, x_1<x_2<\cdots<x_K,\)峰顶点列在平面上严格向右上增长。

从第 \(k\) 个峰往左看，定义它到某个更早峰 \(i<k\) 的视线斜率：\(s_k(i)=\dfrac{y_k-y_i}{x_k-x_i}>0.\)

可见，中间峰 \(j\) 是否遮挡 \(i\) 完全由斜率比较决定。考虑三点 \(i<j<k\)。点 \(P_j\) 在直线 \(\overline{P_iP_k}\) 上方（或在线上）当且仅当\(y_j \ge y_i + \dfrac{y_k-y_i}{x_k-x_i}(x_j-x_i).\)把它改写成以 \(k\) 为右端点的形式可得：\(\dfrac{y_k-y_j}{x_k-x_j} \le \dfrac{y_k-y_i}{x_k-x_i}\Longleftrightarrow s_k(j)\le s_k(i).\)

所以，若存在 \(j\in(i,k)\) 使 \(s_k(j)\le s_k(i)\)，则 \(j\) 在线段 \(\overline{P_iP_k}\) 上方/线上，遮挡 \(i\)，于是 \(i\) 不可见；反之，若对所有 \(j\in(i,k)\) 都有 \(s_k(j)>s_k(i)\)，则所有中间峰都严格在线段下方，\(i\) 可见。

于是得到等价刻画：\(i\) 从 \(k\) 可见，当且仅当\(s_k(i)\)严格小于所有 \(i<j<k\) 的 \(s_k(j)\)。

这意味着：从右往左扫描 \(i=k-1,k-2,\dots,1\)，每当 \(s_k(i)\) 创下目前最小斜率的新纪录时，峰 \(i\) 就可见。因此 \(P(k)\) 就是这些严格更小的新纪录的个数。

不过，如果直接枚举，需要花费\(O(K^2)\)的时间，并不可行。

令第 \(k\) 个峰回望时可见的峰编号（从左到右）为：\(v_1 < v_2 < \cdots < v_t = k-1.\)它们对应的斜率满足：\(s_k(v_1) < s_k(v_2) < \cdots < s_k(v_t)\)

由于可见峰序列的斜率严格单调，上述序列本质上形成了以 \(P_k\) 为观察点的一条单调链：从右向左枚举时，只有当 \(s_k(i)\) 创下更小的新纪录才会入链。因此可以把可见峰列表看成一个单调栈：栈内元素按编号递增，而对应斜率 \(s_k(\cdot)\) 严格递增（等价地，从右向左扫描时新纪录严格递减）。后续算法的核心就是维护并复用这条单调链。

引理：对任意相邻的两个可见峰 \(v_r < v_{r+1}\)，峰 \(v_r\) 从峰 \(v_{r+1}\) 处也是可见的。我们接下来证明这个引理。

取任意 \(j\) 满足 \(v_r<j<v_{r+1}\)。由于 \(v_r\) 与 \(v_{r+1}\) 在 \(k\) 的可见列表中相邻，说明这些 \(j\) 在从 \(k\) 向左扫描时都不是新纪录，因此\(s_k(j) \ge s_k(v_{r+1}).\)如前文推导，\(s_k(j)\ge s_k(v_{r+1})\) 等价于点 \(P_j\) 位于直线 \(\overline{P_{v_{r+1}}P_k}\) 上或下方。

另一方面，因为 \(v_r\) 是更靠左的可见峰，有\(s_k(v_r) < s_k(v_{r+1}),\)把它写成在 \(x=x_{v_r}\) 处的高度比较，等价于点 \(P_{v_r}\) 严格位于直线 \(\overline{P_{v_{r+1}}P_k}\) 的上方。

现在比较两条过点 \(P_{v_{r+1}}\) 的直线：

  • \(\ell_1 = \overline{P_{v_{r+1}}P_k}\)，其斜率为 \(\dfrac{y_k-y_{v_{r+1}}}{x_k-x_{v_{r+1}}}=s_k(v_{r+1})\)；
  • \(\ell_2 = \overline{P_{v_r}P_{v_{r+1}}}\)，它同样过 \(P_{v_{r+1}}\)，并且在更左侧的 \(x=x_{v_r}\) 处经过一个更高的点 \(P_{v_r}\)。

对同一点 \(P_{v_{r+1}}\) 来说，一条直线若在左侧取到更高的值，必然意味着它的斜率更小；因此在区间 \([x_{v_r},x_{v_{r+1}}]\) 上，\(\ell_2\)严格在\(\ell_1\)上方。

而我们已经知道对所有 \(v_r<j<v_{r+1}\)，点 \(P_j\) 都在 \(\ell_1\) 上或下方，于是也必然在 \(\ell_2\) 下方。故所有中间峰都严格在 \(\overline{P_{v_r}P_{v_{r+1}}}\) 下方，按可见性定义，\(v_r\) 从 \(v_{r+1}\) 可见。

引理说明：对固定的 \(k\)，可见峰序列\(v_1 < v_2 < \cdots < v_t=k-1\)满足\(v_r \in V(v_{r+1})(1\le r<t),\)其中\(V\)表示在峰在\(V(\cdot)\)中也能看见，也就是说，每个可见峰都出现在其右侧相邻可见峰的可见列表里。

因此，如果我们已经为每个峰 \(m<k\) 计算并保存了它回望可见的峰列表 \(V(m)\)，那么在计算 \(V(k)\) 时，一定从 \(k-1\) 开始；每当确认某个峰 \(m\) 对 \(k\) 可见，就把 \(V(m)\) 作为更左侧候选加入待检查集合；并且由于相邻可见峰之间的嵌套关系，不会漏掉任何 \(V(k)\) 中的元素。斜率严格变小仍作为唯一筛选标准。

代码

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
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
#include <bits/stdc++.h>
#include "tools.h"
using namespace std;

typedef long long ll;

const int K = 2500000;
vector<int> gen_primes_by_count(int need)
{
    double n = (double)need;
    int limit = (int)(n * (log(n) + log(log(n))) + 20.0);
    vector<int> primes = Factorizer(limit).primes;
    primes.resize(need);
    return primes;
}

bool frac_less(ll a, ll b, ll c, ll d)
{
    return (__int128)a * d < (__int128)c * b;
}

ll solve(int K)
{
    int need_pr = 2 * K;
    vector<int> primes = gen_primes_by_count(need_pr);
    vector<ll> x(K), y(K);
    x[0] = primes[0];
    y[0] = primes[0];
    for (int i = 1; i < K; i++)
    {
        ll b = primes[2 * i - 1];
        ll a = primes[2 * i];
        x[i] = x[i - 1] + a + b;
        y[i] = y[i - 1] + a - b;
    }
    vector<int> start(K), len(K);
    vector<int> all, todo, cur;
    ll ans = 0;
    start[0] = 0;
    len[0] = 0;
    for (int n = 1; n < K; n++)
    {
        todo.clear();
        cur.clear();
        todo.push_back(n - 1);
        ll minNum = 0;
        ll minDen = 0;
        while (!todo.empty())
        {
            int m = todo.back();
            todo.pop_back();

            ll dy = y[n] - y[m];
            ll dx = x[n] - x[m];

            if (minDen == 0 || frac_less(dy, dx, minNum, minDen))
            {
                minNum = dy;
                minDen = dx;

                ans++;
                cur.push_back(m);

                int lm = len[m];
                if (lm)
                {
                    int last = todo.empty() ? INT32_MAX : todo.back();
                    int st = start[m];

                    int pos = 0;
                    if (last != INT32_MAX)
                    {
                        int l = 0, r = lm;
                        while (l < r)
                        {
                            int mid = (l + r) >> 1;
                            int v = all[st + mid];
                            if (v < last)
                                l = mid + 1;
                            else
                                r = mid;
                        }
                        if (l < lm && all[st + l] == last)
                            pos = l + 1;
                        else
                            pos = 0;
                    }

                    for (int t = pos; t < lm; t++)
                        todo.push_back(all[st + t]);
                }
            }
        }

        reverse(cur.begin(), cur.end());

        start[n] = all.size();
        len[n] = cur.size();
        all.insert(all.end(), cur.begin(), cur.end());
    }

    return ans;
}

int main()
{
    cout << solve(K) << "\n";
    return 0;
}