Problem 569: Prime Mountain Range

View on Project Euler

Project Euler Problem 569 Solution

EulerSolve provides an optimized solution for Project Euler Problem 569, Prime Mountain Range, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary Take the consecutive primes \(p_1,p_2,p_3,\dots\) and pair them as \((p_1,p_2),(p_3,p_4),\dots\). Starting from \((0,0)\), move upward with slope \(+1\) for horizontal distance \(p_{2k-1}\), then downward with slope \(-1\) for horizontal distance \(p_{2k}\). The endpoint of the upward segment in the \(k\)-th pair is peak \(k\). For each peak \(k\), define \(P(k)\) as the number of earlier peaks that are visible from \(k\): a peak \(j\lt k\) is visible if the segment joining peaks \(j\) and \(k\) stays strictly above every intermediate peak. The task is to compute $$\sum_{k=1}^{N} P(k),\qquad N=2{,}500{,}000.$$ Mathematical Approach The implemented method turns visibility into an ordered slope problem. Once that reformulation is in place, the search can be pruned with upper convex hull information for prefixes of the peak sequence. Step 1: Write the Peak Coordinates Explicitly Let \((x_k,y_k)\) be the coordinates of peak \(k\). The construction gives $$x_k=\sum_{i=1}^{2k-1} p_i,$$ because \(x\) advances by every prime used so far until the top of the \(k\)-th ascent. Likewise, the height is $$y_k=\sum_{i=1}^{2k-1} (-1)^{i+1} p_i,$$ since odd-indexed primes raise the mountain range and even-indexed primes lower it again....

Detailed mathematical approach

Problem Summary

Take the consecutive primes \(p_1,p_2,p_3,\dots\) and pair them as \((p_1,p_2),(p_3,p_4),\dots\). Starting from \((0,0)\), move upward with slope \(+1\) for horizontal distance \(p_{2k-1}\), then downward with slope \(-1\) for horizontal distance \(p_{2k}\). The endpoint of the upward segment in the \(k\)-th pair is peak \(k\).

For each peak \(k\), define \(P(k)\) as the number of earlier peaks that are visible from \(k\): a peak \(j\lt k\) is visible if the segment joining peaks \(j\) and \(k\) stays strictly above every intermediate peak. The task is to compute

$$\sum_{k=1}^{N} P(k),\qquad N=2{,}500{,}000.$$

Mathematical Approach

The implemented method turns visibility into an ordered slope problem. Once that reformulation is in place, the search can be pruned with upper convex hull information for prefixes of the peak sequence.

Step 1: Write the Peak Coordinates Explicitly

Let \((x_k,y_k)\) be the coordinates of peak \(k\). The construction gives

$$x_k=\sum_{i=1}^{2k-1} p_i,$$

because \(x\) advances by every prime used so far until the top of the \(k\)-th ascent. Likewise, the height is

$$y_k=\sum_{i=1}^{2k-1} (-1)^{i+1} p_i,$$

since odd-indexed primes raise the mountain range and even-indexed primes lower it again. In particular,

$$x_1 \lt x_2 \lt x_3 \lt \cdots,$$

so every visibility question is a comparison between slopes from a fixed right endpoint to points with smaller \(x\)-coordinate.

Step 2: Visibility Is Equivalent to Record-Low Slopes

Fix a peak \(k\) and define the slope from peak \(j\) to peak \(k\) by

$$s_k(j)=\frac{y_k-y_j}{x_k-x_j},\qquad 1\le j\lt k.$$

Now take an intermediate peak \(i\) with \(j\lt i\lt k\). The point \(i\) lies strictly below the chord from \(j\) to \(k\) exactly when

$$y_i \lt y_k-\frac{y_k-y_j}{x_k-x_j}(x_k-x_i).$$

After rearranging, this becomes

$$\frac{y_k-y_i}{x_k-x_i} \gt \frac{y_k-y_j}{x_k-x_j},$$

or in the new notation, \(s_k(i)\gt s_k(j)\). Therefore peak \(j\) is visible from peak \(k\) if and only if

$$s_k(j)\lt s_k(i)\qquad\text{for every }i\text{ with }j\lt i\lt k.$$

So when we scan \(j=k-1,k-2,\dots,1\), the visible peaks are exactly the indices where the slope becomes a new strict minimum. The nearest peak \(k-1\) is always visible, because there is no intermediate peak between \(k-1\) and \(k\).

Step 3: Compare Slopes with Exact Integer Arithmetic

The coordinates become large, so dividing two integers and comparing floating-point values would be risky. Because all denominators are positive, slope comparisons can be rewritten as

$$s_k(j_1)\lt s_k(j_2)\iff (y_k-y_{j_1})(x_k-x_{j_2}) \lt (y_k-y_{j_2})(x_k-x_{j_1}).$$

This is mathematically exact and avoids rounding error. The compiled implementations use wide integer arithmetic for these cross-products, so the visibility test remains exact even at the full problem size.

Step 4: Only the Upper Hull Can Produce a Better Slope

Suppose the current visible peak for the scan of peak \(k\) is \(b\). Any next visible peak must lie somewhere in the prefix \(\{1,2,\dots,b-1\}\). For a fixed right endpoint \(k\), the minimum of \(s_k(j)\) over that prefix is attained on the upper convex hull of the prefix, not at an interior point below the hull.

Geometrically, when a line through peak \(k\) rotates downward, the first point of the prefix that it can touch is a vertex of the upper hull. Therefore:

$$\exists\, j\lt b\text{ with }s_k(j)\lt s_k(b)\iff \exists\, v\text{ on the upper hull of }\{1,\dots,b-1\}\text{ with }s_k(v)\lt s_k(b).$$

This does not immediately identify the next visible peak, but it gives a fast existence test. If no upper-hull vertex improves the slope, the scan can stop at once.

Step 5: Worked Example

The first nine peaks come from the primes \(2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61\). Their coordinates are

$$\begin{aligned} (x_1,y_1)&=(2,2),\\ (x_2,y_2)&=(10,4),\\ (x_3,y_3)&=(28,8),\\ (x_4,y_4)&=(58,12),\\ (x_5,y_5)&=(100,16),\\ (x_6,y_6)&=(160,18),\\ (x_7,y_7)&=(238,22),\\ (x_8,y_8)&=(328,26),\\ (x_9,y_9)&=(440,32). \end{aligned}$$

For \(k=9\), the slopes to earlier peaks are

$$\begin{aligned} s_9(8)&=\frac{32-26}{440-328}=\frac{6}{112},\\ s_9(7)&=\frac{32-22}{440-238}=\frac{10}{202},\\ s_9(6)&=\frac{14}{280},\\ s_9(5)&=\frac{16}{340},\\ s_9(4)&=\frac{20}{382},\\ s_9(3)&=\frac{24}{412},\\ s_9(2)&=\frac{28}{430},\\ s_9(1)&=\frac{30}{438}. \end{aligned}$$

Scanning from right to left, the record-low slopes occur at \(j=8\), then \(j=7\), then \(j=5\). Hence

$$P(9)=3.$$

The same method gives \(P(3)=1\), and the cumulative value

$$\sum_{k=1}^{100} P(k)=227$$

matches the checkpoints used by the optimized implementation.

How the Code Works

The C++, Python, and Java implementations all rely on the same geometry. First, they generate the first \(2N\) primes with an odd-only sieve large enough to cover the required range. Then they build the peak coordinates in a single pass by pairing consecutive primes into one ascent and one descent.

During that same left-to-right pass, the compiled implementations maintain the upper hull of every prefix with a monotone stack. When a new point arrives, the previous hull vertex is removed while the last three hull points fail to form a strict upper turn. The remaining predecessor link is stored so that the upper hull of any prefix can later be traversed backward in constant extra memory per point.

To compute \(P(k)\), the implementation starts with peak \(k-1\), which is automatically visible. It then repeatedly asks whether the earlier prefix contains any point with a smaller slope than the current visible peak. That question is answered by walking only along the stored upper-hull chain of the prefix. If the answer is no, the search ends. If the answer is yes, the implementation walks left until it reaches the first index whose slope is smaller than the current best; that index is the next visible peak, and the process repeats.

The C++ version also checks the first small range against a direct full scan and known checkpoint totals before launching the production run. The C++ and Java versions split the remaining \(k\)-values across several worker threads and add the partial sums. The Python implementation does not reimplement the geometry in pure Python; it delegates to the compiled optimized solver and returns the resulting numeric answer.

Complexity Analysis

Let \(M\) be the sieve limit used to obtain the first \(2N\) primes; in the implementation this is a fixed bound large enough for \(N=2{,}500{,}000\). The sieve costs \(O(M\log\log M)\) time and \(O(M)\) memory in its chosen representation, while building coordinates and prefix hull links costs \(O(N)\) time and \(O(N)\) memory. A direct visibility scan over all pairs of peaks would be \(O(N^2)\). The implemented method keeps the same \(O(N)\) storage for geometric data but prunes most searches through upper-hull existence tests, making the actual Project Euler instance tractable. The code does not present a formal worst-case bound better than quadratic for the search phase, so the honest claim is strong practical acceleration rather than a proved asymptotic improvement.

Footnotes and References

  1. Problem page: Project Euler 569
  2. Prime numbers: Wikipedia - Prime number
  3. Convex hull: Wikipedia - Convex hull
  4. Sieve of Eratosthenes: Wikipedia - Sieve of Eratosthenes
  5. Cross product and orientation tests: Wikipedia - Cross product

Problem 569 source code

C++

#include <cassert>
#include <cstdint>
#include <iostream>
#include <thread>
#include <vector>

namespace {

using i64 = long long;
using u64 = std::uint64_t;

// 128-bit safe multiply-add for cross products / slope comparisons.
using i128 = __int128_t;

// Upper hull maintenance for points with strictly increasing x.
static inline i128 cross(const std::vector<i64>& x, const std::vector<i64>& y, int a, int b, int c) {
    // cross( (b-a), (c-b) )
    const i128 abx = (i128)x[b] - x[a];
    const i128 aby = (i128)y[b] - y[a];
    const i128 bcx = (i128)x[c] - x[b];
    const i128 bcy = (i128)y[c] - y[b];
    return abx * bcy - aby * bcx;
}

// Compare slopes to a fixed right endpoint k:
// slope(k, j1) < slope(k, j2)  <=>  (y_k - y_{j1})/(x_k - x_{j1}) < (y_k - y_{j2})/(x_k - x_{j2})
static inline bool slope_less(const std::vector<i64>& x,
                              const std::vector<i64>& y,
                              int k,
                              int j1,
                              int j2) {
    const i128 dy1 = (i128)y[k] - y[j1];
    const i128 dx1 = (i128)x[k] - x[j1];
    const i128 dy2 = (i128)y[k] - y[j2];
    const i128 dx2 = (i128)x[k] - x[j2];
    return dy1 * dx2 < dy2 * dx1;
}

static std::vector<int> first_primes(std::size_t need) {
    // Sieve up to 1e8: pi(1e8)=5761455 > 5e6, enough for this problem.
    const int limit = 100'000'000;
    const int half = limit / 2;  // odds: n = 2*i+1 for i in [1..half]

    std::vector<bool> comp(half + 1, false);
    std::vector<int> primes;
    primes.reserve(need);
    primes.push_back(2);

    const int sqrt_limit = 10'000;  // floor(sqrt(1e8))
    for (int i = 1; (2 * i + 1) <= sqrt_limit; ++i) {
        if (comp[i]) continue;
        const int p = 2 * i + 1;
        const int start = (p * p) / 2;
        for (int j = start; j <= half; j += p) comp[j] = true;
    }

    for (int i = 1; i <= half && primes.size() < need; ++i) {
        if (!comp[i]) primes.push_back(2 * i + 1);
    }
    assert(primes.size() >= need);
    primes.resize(need);
    return primes;
}

struct Data {
    std::vector<i64> x;
    std::vector<i64> y;
    std::vector<int> hull_prev;  // previous vertex on upper hull of prefix i (valid for all i).
};

static Data build_points_and_prefix_hulls(int N) {
    // Need primes p_1..p_{2N}.
    const std::size_t need = (std::size_t)2 * (std::size_t)N;
    const std::vector<int> p = first_primes(need);

    Data d;
    d.x.assign(N + 1, 0);
    d.y.assign(N + 1, 0);
    d.hull_prev.assign(N + 1, 0);

    i64 valley_y = 0;
    i64 cur_x = 0;
    for (int k = 1; k <= N; ++k) {
        const int up = p[(std::size_t)2 * (std::size_t)k - 2];    // p_{2k-1}
        const int down = p[(std::size_t)2 * (std::size_t)k - 1];  // p_{2k}
        cur_x += up;
        d.x[k] = cur_x;
        d.y[k] = valley_y + up;
        cur_x += down;
        valley_y += (i64)up - (i64)down;
    }

    // Build upper hull for each prefix using the standard monotone chain update, but store predecessor pointers
    // so each prefix hull can be traversed by following hull_prev pointers from its rightmost vertex.
    std::vector<int> st;
    st.reserve(N);
    for (int i = 1; i <= N; ++i) {
        while (st.size() >= 2U) {
            const int a = st[st.size() - 2];
            const int b = st[st.size() - 1];
            if (cross(d.x, d.y, a, b, i) >= 0) {
                st.pop_back();
            } else {
                break;
            }
        }
        d.hull_prev[i] = st.empty() ? 0 : st.back();
        st.push_back(i);
    }

    return d;
}

static int P_full_scan(const Data& d, int k) {
    if (k <= 1) return 0;
    int cnt = 1;
    int best = k - 1;
    for (int j = k - 2; j >= 1; --j) {
        if (slope_less(d.x, d.y, k, j, best)) {
            best = j;
            ++cnt;
        }
    }
    return cnt;
}

static int P_fast(const Data& d, int k) {
    if (k <= 1) return 0;
    int cnt = 1;
    int best = k - 1;
    int t = best;

    while (t > 1) {
        // Existence test: min slope among points in [1..t-1] is achieved on the upper hull of prefix (t-1).
        bool exists = false;
        for (int v = t - 1; v != 0; v = d.hull_prev[v]) {
            if (slope_less(d.x, d.y, k, v, best)) {
                exists = true;
                break;
            }
        }
        if (!exists) break;

        int j = t - 1;
        while (j >= 1 && !slope_less(d.x, d.y, k, j, best)) --j;
        if (j <= 0) break;  // should not happen, but keeps us safe if equality ever occurs.

        ++cnt;
        best = j;
        t = j;
    }
    return cnt;
}

static u64 sum_P_range(const Data& d, int L, int R) {
    u64 acc = 0;
    for (int k = L; k <= R; ++k) acc += (u64)P_fast(d, k);
    return acc;
}

}  // namespace

int main() {
    std::ios::sync_with_stdio(false);
    std::cin.tie(nullptr);

    constexpr int N = 2'500'000;
    const Data d = build_points_and_prefix_hulls(N);

    // Validation by direct scan for the first few peaks.
    u64 sumP = 0;
    for (int k = 1; k <= 200; ++k) {
        const int pk_fast = P_fast(d, k);
        const int pk_full = P_full_scan(d, k);
        assert(pk_fast == pk_full);
        sumP += (u64)pk_fast;
        if (k == 3) assert(pk_fast == 1);
        if (k == 9) assert(pk_fast == 3);
        if (k == 100) assert(sumP == 227);
    }

    const int start = 201;
    if (start <= N) {
        const unsigned hw = std::max(1u, std::thread::hardware_concurrency());
        const unsigned threads = std::min<unsigned>(hw, 8u);
        std::vector<std::thread> pool;
        std::vector<u64> partial(threads, 0);

        const int total = N - start + 1;
        const int chunk = (total + (int)threads - 1) / (int)threads;
        for (unsigned ti = 0; ti < threads; ++ti) {
            const int L = start + (int)ti * chunk;
            const int R = std::min(N, L + chunk - 1);
            if (L > R) continue;
            pool.emplace_back([&d, &partial, ti, L, R]() { partial[ti] = sum_P_range(d, L, R); });
        }
        for (auto& th : pool) th.join();
        for (u64 v : partial) sumP += v;
    }

    std::cout << sumP << '\n';
    return 0;
}

Python

from __future__ import annotations

import re
import shutil
import subprocess
from pathlib import Path

ANSWER_RE = re.compile(r"answer\s*:\s*(.+)$", re.IGNORECASE)
EQUAL_RE = re.compile(r"=\s*(.+)$")


def parse_output(stdout: str) -> str:
    lines = [line.strip() for line in stdout.splitlines() if line.strip()]
    if not lines:
        return ""
    answers = []
    equals = []
    for line in lines:
        m1 = ANSWER_RE.search(line)
        if m1:
            answers.append(m1.group(1).strip())
        m2 = EQUAL_RE.search(line)
        if m2:
            equals.append(m2.group(1).strip())
    if answers:
        return answers[-1]
    if equals:
        return equals[-1]
    return lines[-1]


def should_skip_cpp_checkpoints(src: Path) -> bool:
    try:
        text = src.read_text(encoding="utf-8", errors="ignore")
    except OSError:
        return False
    return "--skip-checkpoints" in text


def run_cpp(binary: Path, src: Path, root: Path) -> str:
    cmd = [str(binary)]
    if should_skip_cpp_checkpoints(src):
        cmd.append("--skip-checkpoints")

    try:
        return subprocess.check_output(cmd, text=True, cwd=root)
    except subprocess.CalledProcessError:
        return subprocess.check_output(cmd, text=True, cwd=src.parent)


def solve() -> str:
    problem_id = __file__.split("Euler")[-1].split(".")[0]
    root = Path(__file__).resolve().parent.parent
    src = root / "solutionsCpp" / f"Euler{problem_id}.cpp"
    binary = root / "solutionsCpp" / f".euler{problem_id}_py_bridge"

    if not binary.exists() or src.stat().st_mtime > binary.stat().st_mtime:
        compiler = shutil.which("clang++") or shutil.which("g++")
        if not compiler:
            raise RuntimeError("No C++ compiler found (clang++/g++).")
        subprocess.check_call([compiler, "-std=c++17", "-O2", str(src), "-o", str(binary)])

    output = run_cpp(binary=binary, src=src, root=root)
    parsed = parse_output(output)
    if not parsed:
        raise RuntimeError(f"Euler{problem_id} bridge produced empty output.")
    return parsed


if __name__ == "__main__":
    print(solve())

Java

import java.util.*;
import java.util.concurrent.*;

public class Euler569 {

    static int mulCmp(long A, long B, long C, long D) {
        long hi1 = Math.multiplyHigh(A, B);
        long lo1 = A * B;
        long hi2 = Math.multiplyHigh(C, D);
        long lo2 = C * D;
        if (hi1 != hi2)
            return hi1 > hi2 ? 1 : -1;
        if (lo1 == lo2)
            return 0;
        return (lo1 + Long.MIN_VALUE > lo2 + Long.MIN_VALUE) ? 1 : -1;
    }

    static boolean slopeLess(long[] x, long[] y, int k, int j1, int j2) {
        long dy1 = y[k] - y[j1];
        long dx1 = x[k] - x[j1];
        long dy2 = y[k] - y[j2];
        long dx2 = x[k] - x[j2];
        return mulCmp(dy1, dx2, dy2, dx1) < 0;
    }

    static boolean crossGe0(long[] x, long[] y, int a, int b, int c) {
        long abx = x[b] - x[a];
        long aby = y[b] - y[a];
        long bcx = x[c] - x[b];
        long bcy = y[c] - y[b];
        return mulCmp(abx, bcy, aby, bcx) >= 0;
    }

    static int[] firstPrimes(int need) {
        int limit = 100000000;
        int half = limit / 2;
        byte[] comp = new byte[half + 1];
        int[] primes = new int[need];
        primes[0] = 2;
        int pCount = 1;

        int sqrtLimit = 10000;
        for (int i = 1; (2 * i + 1) <= sqrtLimit; i++) {
            if (comp[i] == 1)
                continue;
            int p = 2 * i + 1;
            int start = (p * p) / 2;
            for (int j = start; j <= half; j += p) {
                comp[j] = 1;
            }
        }

        for (int i = 1; i <= half && pCount < need; i++) {
            if (comp[i] == 0) {
                primes[pCount++] = 2 * i + 1;
            }
        }
        return primes;
    }

    static class Data {
        long[] x;
        long[] y;
        int[] hullPrev;
    }

    static Data buildData(int N) {
        int need = 2 * N;
        int[] p = firstPrimes(need);

        Data d = new Data();
        d.x = new long[N + 1];
        d.y = new long[N + 1];
        d.hullPrev = new int[N + 1];

        long valleyY = 0;
        long curX = 0;
        for (int k = 1; k <= N; k++) {
            int up = p[2 * k - 2];
            int down = p[2 * k - 1];
            curX += up;
            d.x[k] = curX;
            d.y[k] = valleyY + up;
            curX += down;
            valleyY += (long) up - down;
        }

        int[] st = new int[N];
        int stSz = 0;
        for (int i = 1; i <= N; i++) {
            while (stSz >= 2) {
                int a = st[stSz - 2];
                int b = st[stSz - 1];
                if (crossGe0(d.x, d.y, a, b, i)) {
                    stSz--;
                } else {
                    break;
                }
            }
            d.hullPrev[i] = stSz == 0 ? 0 : st[stSz - 1];
            st[stSz++] = i;
        }
        return d;
    }

    static int P_fast(Data d, int k) {
        if (k <= 1)
            return 0;
        int cnt = 1;
        int best = k - 1;
        int t = best;

        while (t > 1) {
            boolean exists = false;
            for (int v = t - 1; v != 0; v = d.hullPrev[v]) {
                if (slopeLess(d.x, d.y, k, v, best)) {
                    exists = true;
                    break;
                }
            }
            if (!exists)
                break;

            int j = t - 1;
            while (j >= 1 && !slopeLess(d.x, d.y, k, j, best))
                j--;
            if (j <= 0)
                break;

            cnt++;
            best = j;
            t = j;
        }
        return cnt;
    }

    public static String solve() {
        int N = 2500000;
        Data d = buildData(N);

        int threads = Runtime.getRuntime().availableProcessors();
        if (threads <= 0)
            threads = 1;
        ExecutorService executor = Executors.newFixedThreadPool(threads);
        List<Callable<Long>> tasks = new ArrayList<>();

        int start = 1;
        int total = N - start + 1;
        int chunk = (total + threads - 1) / threads;

        for (int ti = 0; ti < threads; ti++) {
            final int L = start + ti * chunk;
            final int R = Math.min(N, L + chunk - 1);
            if (L > R)
                continue;
            tasks.add(() -> {
                long sum = 0;
                for (int k = L; k <= R; k++) {
                    sum += P_fast(d, k);
                }
                return sum;
            });
        }

        long sumP = 0;
        try {
            for (Future<Long> res : executor.invokeAll(tasks)) {
                sumP += res.get();
            }
        } catch (Exception e) {
        }
        executor.shutdown();

        return String.valueOf(sumP);
    }

    public static void main(String[] args) {
        System.out.println(solve());
    }
}