Problem 576: Irrational Jumps

Problem 576 source code

C++

#include <algorithm>
#include <cmath>
#include <cstdint>
#include <iomanip>
#include <iostream>
#include <map>
#include <numeric>
#include <vector>

// Project Euler 576: Irrational Jumps
//
// For fixed irrational jump length l, the positions are x_k = frac(k*l).
// A gap of length g starting at d catches the walker at the first k with x_k in [d, d+g).
// Equivalently, for fixed x_k, it catches gaps with d in (x_k - g, x_k], intersected with (0, 1-g).
//
// For each l we build the piecewise-constant function T_l(d) = first hit time (number of jumps),
// by "painting" the interval of d-values hit at time k onto the still-uncovered domain, increasing k.
// Then S(l,g,d) = l * T_l(d).
//
// For M(n,g), we need max over d of sum_{p<=n, p prime} S(sqrt(1/p), g, d).
// Each S_p(d) is piecewise-constant; their sum is too. The maximum occurs on one of the induced
// subintervals, so we sweep all segment-start events across primes.

using u64 = std::uint64_t;

static std::vector<int> primes_up_to(int n) {
    std::vector<bool> is_prime(n + 1, true);
    if (n >= 0) is_prime[0] = false;
    if (n >= 1) is_prime[1] = false;
    for (int i = 2; 1LL * i * i <= n; ++i) {
        if (!is_prime[i]) continue;
        for (int j = i * i; j <= n; j += i) is_prime[j] = false;
    }
    std::vector<int> ps;
    for (int i = 2; i <= n; ++i)
        if (is_prime[i]) ps.push_back(i);
    return ps;
}

struct Seg {
    long double l = 0;
    long double r = 0;
    int k = 0;  // first hit time
};

static std::vector<Seg> build_first_hit_segments(long double alpha, long double g) {
    // Domain for d: (0, 1-g). We'll work on [0, 1-g] and ignore measure-zero endpoint issues.
    const long double domL = 0.0L;
    const long double domR = 1.0L - g;
    const long double eps = 1e-21L;

    std::map<long double, long double> uncovered;
    uncovered.emplace(domL, domR);
    long double remaining = domR - domL;

    std::vector<Seg> segs;
    segs.reserve((size_t)(1.0L / g) + 10);

    long double x = 0.0L;
    int k = 0;
    while (!uncovered.empty() && remaining > eps) {
        ++k;
        x += alpha;
        x -= floorl(x);  // frac

        long double L = x - g;
        long double U = x;
        if (L < domL) L = domL;
        if (U > domR) U = domR;
        if (L + eps >= U) continue;

        // Iterate uncovered intervals that intersect [L, U].
        auto it = uncovered.upper_bound(L);
        if (it != uncovered.begin()) --it;
        while (it != uncovered.end()) {
            const long double a = it->first;
            const long double b = it->second;
            if (b <= L + eps) {
                ++it;
                continue;
            }
            if (a >= U - eps) break;

            const long double ol = std::max(a, L);
            const long double orr = std::min(b, U);
            if (ol + eps < orr) {
                segs.push_back(Seg{ol, orr, k});
                remaining -= (orr - ol);
            }

            // Remove current interval and reinsert leftovers.
            auto it_erase = it++;
            uncovered.erase(it_erase);
            if (a + eps < L) uncovered.emplace(a, std::min(L, b));
            if (U + eps < b) uncovered.emplace(std::max(U, a), b);
        }
    }

    // Sort and merge adjacent segments with the same k (should be rare).
    std::sort(segs.begin(), segs.end(), [](const Seg& A, const Seg& B) { return A.l < B.l; });
    std::vector<Seg> merged;
    merged.reserve(segs.size());
    for (const auto& s : segs) {
        if (merged.empty()) {
            merged.push_back(s);
            continue;
        }
        auto& last = merged.back();
        if (s.k == last.k && fabsl(s.l - last.r) < 1e-18L) {
            last.r = s.r;
        } else {
            merged.push_back(s);
        }
    }

    // Ensure coverage includes the left boundary by filling any tiny gaps (numerical).
    // The exact maximum is unaffected by eps-sized gaps.
    return merged;
}

struct Event {
    long double pos = 0;
    int idx = 0;
    long double new_val = 0;
};

static long double compute_M(int n, long double g) {
    const auto ps = primes_up_to(n);
    const long double domL = 0.0L;
    const long double domR = 1.0L - g;

    const int m = (int)ps.size();
    std::vector<long double> weight(m);
    for (int i = 0; i < m; ++i) weight[i] = sqrtl(1.0L / (long double)ps[i]);

    std::vector<long double> cur(m, 0.0L);
    std::vector<Event> events;
    events.reserve(2000000);

    for (int i = 0; i < m; ++i) {
        const long double alpha = weight[i];
        auto segs = build_first_hit_segments(alpha, g);
        if (segs.empty() || segs.front().l > domL + 1e-15L) {
            std::cerr << "Segment construction failed to cover domain for p=" << ps[i] << "\n";
            std::exit(1);
        }

        // Find the segment covering domL (should be first).
        int j0 = 0;
        while (j0 + 1 < (int)segs.size() && segs[j0].r <= domL) ++j0;
        cur[i] = weight[i] * (long double)segs[j0].k;

        // Create events at subsequent segment starts.
        for (int j = j0 + 1; j < (int)segs.size(); ++j) {
            events.push_back(Event{segs[j].l, i, weight[i] * (long double)segs[j].k});
        }
    }

    std::sort(events.begin(), events.end(), [](const Event& a, const Event& b) { return a.pos < b.pos; });

    long double sum = 0.0L;
    for (long double v : cur) sum += v;

    long double best = sum;
    long double pos_prev = domL;

    size_t idx = 0;
    while (idx < events.size()) {
        const long double pos = events[idx].pos;
        if (pos > pos_prev + 1e-18L) best = std::max(best, sum);

        // apply all events at this pos
        while (idx < events.size() && fabsl(events[idx].pos - pos) < 1e-18L) {
            const int i = events[idx].idx;
            sum += events[idx].new_val - cur[i];
            cur[i] = events[idx].new_val;
            ++idx;
        }
        pos_prev = pos;
        if (pos_prev > domR) break;
    }
    best = std::max(best, sum);
    return best;
}

int main() {
    {
        const long double got = compute_M(3, 0.06L);
        const long double expected = 29.5425L;
        if (fabsl(got - expected) > 5e-4L) {
            std::cerr << "Validation failed: M(3, 0.06)\n";
            std::cerr << std::fixed << std::setprecision(6) << (double)got << " vs " << (double)expected << "\n";
            return 1;
        }
    }
    {
        const long double got = compute_M(10, 0.01L);
        const long double expected = 266.9010L;
        if (fabsl(got - expected) > 5e-4L) {
            std::cerr << "Validation failed: M(10, 0.01)\n";
            std::cerr << std::fixed << std::setprecision(6) << (double)got << " vs " << (double)expected << "\n";
            return 1;
        }
    }

    const long double ans = compute_M(100, 0.00002L);
    std::cout << std::fixed << std::setprecision(4) << (double)ans << "\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.*;

public class Euler576 {
    static int[] primesUpTo(int n) {
        boolean[] np = new boolean[n + 1];
        np[0] = np[1] = true;
        for (int i = 2; (long) i * i <= n; i++)
            if (!np[i])
                for (int j = i * i; j <= n; j += i)
                    np[j] = true;
        int c = 0;
        for (int i = 2; i <= n; i++)
            if (!np[i])
                c++;
        int[] p = new int[c];
        c = 0;
        for (int i = 2; i <= n; i++)
            if (!np[i])
                p[c++] = i;
        return p;
    }

    static double[][] buildFirstHitSegments(double alpha, double g) {
        double domR = 1.0 - g;
        TreeMap<Double, Double> uncovered = new TreeMap<>();
        uncovered.put(0.0, domR);
        double remaining = domR;
        List<double[]> segs = new ArrayList<>();
        double x = 0;
        int k = 0;
        while (!uncovered.isEmpty() && remaining > 1e-15) {
            k++;
            x += alpha;
            x -= Math.floor(x);
            double L = Math.max(x - g, 0.0), U = Math.min(x, domR);
            if (L + 1e-18 >= U)
                continue;
            Map.Entry<Double, Double> entry = uncovered.floorEntry(L);
            if (entry == null)
                entry = uncovered.firstEntry();
            while (entry != null) {
                double a = entry.getKey(), b = entry.getValue();
                if (b <= L + 1e-18) {
                    entry = uncovered.higherEntry(a);
                    continue;
                }
                if (a >= U - 1e-18)
                    break;
                double ol = Math.max(a, L), orr = Math.min(b, U);
                if (ol + 1e-18 < orr) {
                    segs.add(new double[] { ol, orr, k });
                    remaining -= (orr - ol);
                }
                uncovered.remove(a);
                if (a + 1e-18 < L)
                    uncovered.put(a, Math.min(L, b));
                if (U + 1e-18 < b)
                    uncovered.put(Math.max(U, a), b);
                entry = uncovered.higherEntry(a);
            }
        }
        segs.sort(Comparator.comparingDouble(s -> s[0]));
        List<double[]> merged = new ArrayList<>();
        for (double[] s : segs) {
            if (merged.isEmpty()) {
                merged.add(s);
                continue;
            }
            double[] last = merged.get(merged.size() - 1);
            if ((int) s[2] == (int) last[2] && Math.abs(s[0] - last[1]) < 1e-15)
                last[1] = s[1];
            else
                merged.add(s);
        }
        return merged.toArray(new double[0][]);
    }

    public static void main(String[] args) {
        int n = 100;
        double g = 0.00002;
        int[] ps = primesUpTo(n);
        int m = ps.length;
        double[] weight = new double[m];
        for (int i = 0; i < m; i++)
            weight[i] = Math.sqrt(1.0 / ps[i]);
        double[] cur = new double[m];
        List<double[]> events = new ArrayList<>(); // [pos, idx, newVal]
        for (int i = 0; i < m; i++) {
            double[][] segs = buildFirstHitSegments(weight[i], g);
            cur[i] = weight[i] * segs[0][2];
            for (int j = 1; j < segs.length; j++)
                events.add(new double[] { segs[j][0], i, weight[i] * segs[j][2] });
        }
        events.sort(Comparator.comparingDouble(e -> e[0]));
        double sum = 0;
        for (double v : cur)
            sum += v;
        double best = sum;
        for (int idx = 0; idx < events.size();) {
            double pos = events.get(idx)[0];
            while (idx < events.size() && Math.abs(events.get(idx)[0] - pos) < 1e-18) {
                int i = (int) events.get(idx)[1];
                sum += events.get(idx)[2] - cur[i];
                cur[i] = events.get(idx)[2];
                idx++;
            }
            if (sum > best)
                best = sum;
        }
        System.out.printf("%.4f%n", best);
    }
}