Problem 689: Binary Series

View on Project Euler

Project Euler Problem 689 Solution

EulerSolve provides an optimized solution for Project Euler Problem 689, Binary Series, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary Let $$Z=\sum_{n\ge 1}\frac{B_n}{n^2},\qquad B_n\in\{0,1\},\qquad \Pr(B_n=0)=\Pr(B_n=1)=\tfrac12.$$ The problem asks for the tail probability $$p(a)=\Pr(Z\gt a)$$ at \(a=0.5\), to high precision. Since \(0\le Z\le \zeta(2)=\pi^2/6\) and the series contains infinitely many independent bits, direct enumeration is hopeless. The implementation instead rewrites the distribution in a symmetric form, derives its characteristic function, and then recovers the probability by Fourier inversion. Mathematical Approach The whole method is a clean chain: center the random series, compute the characteristic function term by term, convert the tail probability into a one-dimensional integral, and then approximate that integral numerically. Step 1: Center the Series and Identify the Range Each random bit contributes either \(0\) or \(1/n^2\), so its mean is \(1/(2n^2)\). Therefore $$\mu=\mathbb E[Z]=\sum_{n\ge 1}\frac{1}{2n^2}=\frac{\zeta(2)}{2}=\frac{\pi^2}{12}.$$ Write \(B_n=\frac{1+\varepsilon_n}{2}\), where \(\varepsilon_n\in\{-1,+1\}\) with equal probability. Then $$Z=\mu+\sum_{n\ge 1}\frac{\varepsilon_n}{2n^2}.$$ So the centered variable is $$Y=Z-\mu=\sum_{n\ge 1}\frac{\varepsilon_n}{2n^2},$$ which is symmetric around \(0\)....

Detailed mathematical approach

Problem Summary

Let

$$Z=\sum_{n\ge 1}\frac{B_n}{n^2},\qquad B_n\in\{0,1\},\qquad \Pr(B_n=0)=\Pr(B_n=1)=\tfrac12.$$

The problem asks for the tail probability

$$p(a)=\Pr(Z\gt a)$$

at \(a=0.5\), to high precision. Since \(0\le Z\le \zeta(2)=\pi^2/6\) and the series contains infinitely many independent bits, direct enumeration is hopeless. The implementation instead rewrites the distribution in a symmetric form, derives its characteristic function, and then recovers the probability by Fourier inversion.

Mathematical Approach

The whole method is a clean chain: center the random series, compute the characteristic function term by term, convert the tail probability into a one-dimensional integral, and then approximate that integral numerically.

Step 1: Center the Series and Identify the Range

Each random bit contributes either \(0\) or \(1/n^2\), so its mean is \(1/(2n^2)\). Therefore

$$\mu=\mathbb E[Z]=\sum_{n\ge 1}\frac{1}{2n^2}=\frac{\zeta(2)}{2}=\frac{\pi^2}{12}.$$

Write \(B_n=\frac{1+\varepsilon_n}{2}\), where \(\varepsilon_n\in\{-1,+1\}\) with equal probability. Then

$$Z=\mu+\sum_{n\ge 1}\frac{\varepsilon_n}{2n^2}.$$

So the centered variable is

$$Y=Z-\mu=\sum_{n\ge 1}\frac{\varepsilon_n}{2n^2},$$

which is symmetric around \(0\). Also every term is nonnegative and the largest possible sum is \(\zeta(2)\), so

$$a\lt 0 \implies p(a)=1,\qquad a\ge \zeta(2)\implies p(a)=0.$$

Step 2: Derive the Characteristic Function

For a fixed \(n\), the centered summand has characteristic function

$$\mathbb E\left[e^{it\varepsilon_n/(2n^2)}\right]=\frac{e^{it/(2n^2)}+e^{-it/(2n^2)}}{2}=\cos\left(\frac{t}{2n^2}\right).$$

The summands are independent, so the characteristic function of the full series factorizes into an infinite product:

$$\varphi_Y(t)=\prod_{n\ge 1}\cos\left(\frac{t}{2n^2}\right).$$

This function is real and even, reflecting the symmetry \(Y\overset{d}= -Y\). That symmetry is the reason the final inversion formula contains only a sine factor and no imaginary part to evaluate separately.

Step 3: Convert the Tail Probability into an Integral

Set

$$x=\mu-a.$$

Then \(p(a)=\Pr(Y\gt -x)\). Gil-Pelaez inversion gives

$$p(a)=\frac12+\frac{1}{\pi}\int_0^{\infty}\frac{\sin(tx)}{t}\,\varphi_Y(t)\,dt.$$

So an infinite random sum has been reduced to a deterministic one-dimensional integral. The hard combinatorics disappear; the remaining work is numerical analysis.

Step 4: Truncate the Infinite Product and the Infinite Interval

The implementation cannot keep infinitely many cosine factors, so it replaces \(\varphi_Y(t)\) by

$$\varphi_{Y,N}(t)=\prod_{n=1}^{N}\cos\left(\frac{t}{2n^2}\right),$$

with a large cutoff \(N\). It also replaces the interval \([0,\infty)\) by a finite interval \([0,T]\). This is effective because the omitted cosine factors are extremely close to \(1\) once \(n\) is large, and because the oscillatory integral contributes less and less beyond a sufficiently large upper limit.

Step 5: Apply Simpson's Rule

If \(T=mh\) with even \(m\), the numerical quadrature uses

$$\int_0^T f(t)\,dt\approx \frac{h}{3}\left(f(0)+f(T)+4\sum_{j=1,3,\dots,m-1}f(jh)+2\sum_{j=2,4,\dots,m-2}f(jh)\right),$$

where

$$f(t)=\frac{\sin(tx)}{t}\,\varphi_{Y,N}(t).$$

At \(t=0\), the expression \(\sin(tx)/t\) is interpreted by continuity:

$$\lim_{t\to 0}\frac{\sin(tx)}{t}=x.$$

That removes the apparent singularity at the origin and makes Simpson's rule stable there.

Worked Example: Symmetry and the Main Numerical Value

If every bit \(B_n\) is replaced by \(1-B_n\), then

$$\sum_{n\ge 1}\frac{1-B_n}{n^2}=\zeta(2)-Z.$$

So the distribution of \(Z\) is symmetric around \(\mu=\zeta(2)/2\). It follows that

$$p(\mu)=\Pr(Z\gt \mu)=\frac12,$$

and more generally

$$p(\mu-u)+p(\mu+u)=1.$$

In particular, the thresholds \(1\) and \(\zeta(2)-1\) are symmetric around \(\mu\), so both have probability \(1/2\). With the finer numerical settings used by the implementation, the target quantity comes out as

$$p(0.5)\approx 0.565654540708545,$$

which is the value checked before the final rounded answer is printed.

How the Code Works

The C++, Python, and Java implementations all compute the same truncated Fourier integral. The C++ implementation performs the numerical work directly: it precomputes the coefficients \(1/(2n^2)\) up to the cutoff, evaluates the truncated cosine product at each Simpson node, and accumulates the weighted node values over the chosen interval.

The interior Simpson nodes are divided across available threads, so the wall-clock time improves on multi-core hardware while the mathematical result stays unchanged. The endpoint at \(t=0\) is handled through the limit above, and the number of subintervals is adjusted so Simpson's rule always uses an even count.

The Python and Java implementations are thin front ends that delegate to the same compiled numerical core, so all three languages share the same checkpoints, the same truncation strategy, and the same final probability once their numerical parameters agree.

Complexity Analysis

If the integration limit is \(T\), the step size is \(h\), and the cosine-product cutoff is \(N\), then the number of Simpson subintervals is about \(T/h\). Each node evaluation multiplies \(N\) cosine factors, so the total running time is

$$O\left(\frac{T}{h}\,N\right).$$

The memory cost is \(O(N)\) for the precomputed coefficients, plus a small extra buffer for per-thread partial sums. Multithreading changes the constant factor in runtime, not the asymptotic order.

Footnotes and References

  1. Project Euler problem page: https://projecteuler.net/problem=689
  2. Characteristic function: Wikipedia - Characteristic function
  3. Gil-Pelaez inversion theorem: Wikipedia - Gil-Pelaez theorem
  4. Simpson's rule: Wikipedia - Simpson's rule
  5. Basel problem and \(\zeta(2)=\pi^2/6\): Wikipedia - Basel problem

Problem 689 source code

C++

#include <algorithm>
#include <chrono>
#include <cerrno>
#include <cmath>
#include <cstdint>
#include <cstdlib>
#include <iomanip>
#include <iostream>
#include <limits>
#include <string>
#include <thread>
#include <vector>

namespace {

using u32 = std::uint32_t;

constexpr long double kDefaultA = 0.5L;
constexpr long double kDefaultTMax = 2000.0L;
constexpr long double kDefaultStep = 0.04L;
constexpr int kDefaultCutoff = 10'000;

constexpr long double kCoarseTMax = 1200.0L;
constexpr long double kCoarseStep = 0.05L;
constexpr int kCoarseCutoff = 4'000;

constexpr long double kFineTMax = 2000.0L;
constexpr long double kFineStep = 0.04L;
constexpr int kFineCutoff = 10'000;
constexpr long double kCheckpointExpectedMain = 0.565654540708545L;

struct Options {
    long double a = kDefaultA;
    long double t_max = kDefaultTMax;
    long double step = kDefaultStep;
    int cutoff = kDefaultCutoff;
    bool run_checkpoints = true;
    bool allow_multithreading = true;
    unsigned requested_threads = 0U;
};

struct Kernel {
    explicit Kernel(const int cutoff)
        : cutoff_(cutoff), inv_half_squares_(static_cast<std::size_t>(cutoff) + 1ULL, 0.0L) {
        for (int n = 1; n <= cutoff_; ++n) {
            const long double nn = static_cast<long double>(n);
            inv_half_squares_[static_cast<std::size_t>(n)] = 0.5L / (nn * nn);
        }
    }

    long double psi(const long double t) const {
        long double p = 1.0L;
        for (int n = 1; n <= cutoff_; ++n) {
            p *= cosl(t * inv_half_squares_[static_cast<std::size_t>(n)]);
        }
        return p;
    }

  private:
    int cutoff_;
    std::vector<long double> inv_half_squares_;
};

unsigned choose_thread_count(const bool allow_multithreading,
                             const unsigned requested_threads,
                             const std::size_t workload) {
    if (!allow_multithreading || workload < 2ULL) {
        return 1U;
    }

    unsigned threads = requested_threads;
    if (threads == 0U) {
        threads = std::thread::hardware_concurrency();
        if (threads == 0U) {
            threads = 1U;
        }
    }

    return std::max(1U, std::min<unsigned>(threads, static_cast<unsigned>(workload)));
}

bool parse_u32_after_prefix(const std::string& arg, const char* prefix, u32& value) {
    const std::string p(prefix);
    if (arg.rfind(p, 0) != 0U) {
        return false;
    }

    const std::string tail = arg.substr(p.size());
    if (tail.empty()) {
        return false;
    }

    std::uint64_t parsed = 0ULL;
    for (const char c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }
        parsed = parsed * 10ULL + static_cast<std::uint64_t>(c - '0');
        if (parsed > static_cast<std::uint64_t>(std::numeric_limits<u32>::max())) {
            return false;
        }
    }

    value = static_cast<u32>(parsed);
    return true;
}

bool parse_unsigned_after_prefix(const std::string& arg,
                                 const char* prefix,
                                 unsigned& value) {
    u32 parsed = 0U;
    if (!parse_u32_after_prefix(arg, prefix, parsed)) {
        return false;
    }
    value = static_cast<unsigned>(parsed);
    return true;
}

bool parse_long_double_after_prefix(const std::string& arg,
                                    const char* prefix,
                                    long double& value) {
    const std::string p(prefix);
    if (arg.rfind(p, 0) != 0U) {
        return false;
    }

    const std::string tail = arg.substr(p.size());
    if (tail.empty()) {
        return false;
    }

    char* end = nullptr;
    errno = 0;
    const long double parsed = std::strtold(tail.c_str(), &end);
    if (end == tail.c_str() || *end != '\0' || errno == ERANGE || !std::isfinite(parsed)) {
        return false;
    }

    value = parsed;
    return true;
}

bool parse_arguments(const int argc, char** argv, Options& options) {
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);

        if (arg == "--skip-checkpoints") {
            options.run_checkpoints = false;
            continue;
        }
        if (arg == "--single-thread") {
            options.allow_multithreading = false;
            continue;
        }

        long double parsed_ld = 0.0L;
        if (parse_long_double_after_prefix(arg, "--a=", parsed_ld)) {
            options.a = parsed_ld;
            continue;
        }
        if (parse_long_double_after_prefix(arg, "--t-max=", parsed_ld)) {
            options.t_max = parsed_ld;
            continue;
        }
        if (parse_long_double_after_prefix(arg, "--step=", parsed_ld)) {
            options.step = parsed_ld;
            continue;
        }

        u32 parsed_u32 = 0U;
        if (parse_u32_after_prefix(arg, "--cutoff=", parsed_u32)) {
            if (parsed_u32 > static_cast<u32>(std::numeric_limits<int>::max())) {
                std::cerr << "--cutoff is too large.\n";
                return false;
            }
            options.cutoff = static_cast<int>(parsed_u32);
            continue;
        }

        unsigned parsed_unsigned = 0U;
        if (parse_unsigned_after_prefix(arg, "--threads=", parsed_unsigned)) {
            options.requested_threads = parsed_unsigned;
            continue;
        }

        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }

    if (options.step <= 0.0L) {
        std::cerr << "--step must be positive.\n";
        return false;
    }
    if (options.t_max <= 0.0L) {
        std::cerr << "--t-max must be positive.\n";
        return false;
    }
    if (options.cutoff < 1) {
        std::cerr << "--cutoff must be at least 1.\n";
        return false;
    }

    return true;
}

long double simpson_integral_parallel(const long double x,
                                      const long double t_max,
                                      const long double step,
                                      const int cutoff,
                                      const bool allow_multithreading,
                                      const unsigned requested_threads) {
    std::size_t segments = static_cast<std::size_t>(llround(t_max / step));
    if ((segments & 1ULL) != 0ULL) {
        ++segments;
    }

    const std::size_t interior_count = (segments > 1ULL) ? (segments - 1ULL) : 0ULL;
    const unsigned threads =
        choose_thread_count(allow_multithreading, requested_threads, interior_count);

    const Kernel kernel(cutoff);

    auto eval = [&](const std::size_t i) {
        if (i == 0ULL) {
            return x;
        }
        const long double t = static_cast<long double>(i) * step;
        return (sinl(t * x) / t) * kernel.psi(t);
    };

    std::vector<long double> partials(threads, 0.0L);
    std::vector<std::thread> pool;
    pool.reserve(threads);

    for (unsigned tid = 0U; tid < threads; ++tid) {
        const std::size_t begin_offset = (interior_count * tid) / threads;
        const std::size_t end_offset = (interior_count * (tid + 1U)) / threads;

        const std::size_t begin_index = 1ULL + begin_offset;
        const std::size_t end_index = 1ULL + end_offset;

        pool.emplace_back([&, tid, begin_index, end_index]() {
            long double local = 0.0L;
            for (std::size_t i = begin_index; i < end_index; ++i) {
                const long double fi = eval(i);
                local += ((i & 1ULL) != 0ULL) ? (4.0L * fi) : (2.0L * fi);
            }
            partials[tid] = local;
        });
    }

    for (auto& th : pool) {
        th.join();
    }

    long double total = eval(0ULL) + eval(segments);
    for (const long double part : partials) {
        total += part;
    }

    return total * step / 3.0L;
}

long double probability_greater_than(const long double a,
                                     const long double t_max,
                                     const long double step,
                                     const int cutoff,
                                     const bool allow_multithreading,
                                     const unsigned requested_threads) {
    const long double pi = acosl(-1.0L);
    const long double zeta2 = (pi * pi) / 6.0L;
    const long double mu = zeta2 / 2.0L;

    if (a < 0.0L) {
        return 1.0L;
    }
    if (a >= zeta2) {
        return 0.0L;
    }

    const long double x = mu - a;
    const long double integral =
        simpson_integral_parallel(x, t_max, step, cutoff, allow_multithreading, requested_threads);

    return 0.5L + integral / pi;
}

bool close_to(const long double lhs, const long double rhs, const long double tolerance) {
    return fabsl(lhs - rhs) <= tolerance;
}

bool run_checkpoints(const Options& options) {
    const long double pi = acosl(-1.0L);
    const long double zeta2 = (pi * pi) / 6.0L;
    const long double mu = zeta2 / 2.0L;

    const long double p_mu = probability_greater_than(mu,
                                                      kFineTMax,
                                                      kFineStep,
                                                      kFineCutoff,
                                                      options.allow_multithreading,
                                                      options.requested_threads);
    if (!close_to(p_mu, 0.5L, 2.0e-12L)) {
        std::cerr << "Checkpoint failed: p(mu) should be 0.5, got " << std::setprecision(18)
                  << static_cast<double>(p_mu) << "\n";
        return false;
    }

    const long double p_zeta2_minus_1 = probability_greater_than(zeta2 - 1.0L,
                                                                  kFineTMax,
                                                                  kFineStep,
                                                                  kFineCutoff,
                                                                  options.allow_multithreading,
                                                                  options.requested_threads);
    if (!close_to(p_zeta2_minus_1, 0.5L, 2.0e-12L)) {
        std::cerr << "Checkpoint failed: p(zeta(2)-1) should be 0.5, got "
                  << std::setprecision(18) << static_cast<double>(p_zeta2_minus_1) << "\n";
        return false;
    }

    const long double p_one = probability_greater_than(1.0L,
                                                        kFineTMax,
                                                        kFineStep,
                                                        kFineCutoff,
                                                        options.allow_multithreading,
                                                        options.requested_threads);
    if (!close_to(p_one, 0.5L, 2.0e-12L)) {
        std::cerr << "Checkpoint failed: p(1) should be 0.5, got " << std::setprecision(18)
                  << static_cast<double>(p_one) << "\n";
        return false;
    }

    const long double symmetry_shift = 0.2L;
    const long double p_left = probability_greater_than(mu - symmetry_shift,
                                                         kCoarseTMax,
                                                         kCoarseStep,
                                                         kCoarseCutoff,
                                                         options.allow_multithreading,
                                                         options.requested_threads);
    const long double p_right = probability_greater_than(mu + symmetry_shift,
                                                          kCoarseTMax,
                                                          kCoarseStep,
                                                          kCoarseCutoff,
                                                          options.allow_multithreading,
                                                          options.requested_threads);
    if (!close_to(p_left + p_right, 1.0L, 2.0e-11L)) {
        std::cerr << "Checkpoint failed: symmetry mismatch, p(mu-u)+p(mu+u)="
                  << std::setprecision(18) << static_cast<double>(p_left + p_right) << "\n";
        return false;
    }

    const long double coarse_main = probability_greater_than(0.5L,
                                                              kCoarseTMax,
                                                              kCoarseStep,
                                                              kCoarseCutoff,
                                                              options.allow_multithreading,
                                                              options.requested_threads);
    const long double fine_main = probability_greater_than(0.5L,
                                                            kFineTMax,
                                                            kFineStep,
                                                            kFineCutoff,
                                                            options.allow_multithreading,
                                                            options.requested_threads);

    if (!close_to(fine_main, kCheckpointExpectedMain, 5.0e-11L)) {
        std::cerr << "Checkpoint failed: p(0.5) reference mismatch, got "
                  << std::setprecision(18) << static_cast<double>(fine_main) << "\n";
        return false;
    }

    if (!close_to(coarse_main, fine_main, 2.0e-10L)) {
        std::cerr << "Checkpoint failed: coarse/fine mismatch, coarse=" << std::setprecision(18)
                  << static_cast<double>(coarse_main) << ", fine=" << static_cast<double>(fine_main)
                  << "\n";
        return false;
    }

    const unsigned thread_probe =
        choose_thread_count(true, options.requested_threads, static_cast<std::size_t>(50'000ULL));
    if (options.allow_multithreading && thread_probe > 1U) {
        const long double single_thread =
            probability_greater_than(0.5L, kCoarseTMax, kCoarseStep, kCoarseCutoff, false, 1U);
        const long double multi_thread = probability_greater_than(0.5L,
                                                                   kCoarseTMax,
                                                                   kCoarseStep,
                                                                   kCoarseCutoff,
                                                                   true,
                                                                   thread_probe);
        if (!close_to(single_thread, multi_thread, 1.0e-13L)) {
            std::cerr << "Checkpoint failed: single/multi thread mismatch, single="
                      << std::setprecision(18) << static_cast<double>(single_thread)
                      << ", multi=" << static_cast<double>(multi_thread) << "\n";
            return false;
        }
    }

    return true;
}

} // namespace

int main(int argc, char** argv) {
    Options options;
    if (!parse_arguments(argc, argv, options)) {
        return 1;
    }

    const auto start = std::chrono::steady_clock::now();

    if (options.run_checkpoints && !run_checkpoints(options)) {
        return 1;
    }

    const long double probability = probability_greater_than(options.a,
                                                              options.t_max,
                                                              options.step,
                                                              options.cutoff,
                                                              options.allow_multithreading,
                                                              options.requested_threads);

    const auto finish = std::chrono::steady_clock::now();
    const double elapsed = std::chrono::duration<double>(finish - start).count();

    if (options.run_checkpoints) {
        std::cout << "Checkpoints passed.\n";
    }

    std::cout << std::fixed << std::setprecision(12)
              << "p(" << static_cast<double>(options.a) << ") = "
              << static_cast<double>(probability) << "\n";

    std::cout << std::fixed << std::setprecision(8)
              << "Rounded to 8 decimals: " << static_cast<double>(probability) << "\n";

    std::cout << std::fixed << std::setprecision(6)
              << "Elapsed: " << elapsed << " s\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.nio.file.*;
import java.util.*;
import java.util.regex.*;

public class Euler689 {
    private static final Pattern ANSWER_RE = Pattern.compile("answer\\s*:\\s*(.+)$", Pattern.CASE_INSENSITIVE);
    private static final Pattern EQUAL_RE = Pattern.compile("=\\s*(.+)$");

    private static String parseOutput(String stdout) {
        String[] lines = stdout.split("\\R");
        List<String> nonEmpty = new ArrayList<>();
        for (String line : lines) {
            String t = line.trim();
            if (!t.isEmpty()) {
                nonEmpty.add(t);
            }
        }
        if (nonEmpty.isEmpty()) {
            return "";
        }

        List<String> answers = new ArrayList<>();
        List<String> equals = new ArrayList<>();
        for (String line : nonEmpty) {
            Matcher m1 = ANSWER_RE.matcher(line);
            if (m1.find()) {
                answers.add(m1.group(1).trim());
            }
            Matcher m2 = EQUAL_RE.matcher(line);
            if (m2.find()) {
                equals.add(m2.group(1).trim());
            }
        }

        if (!answers.isEmpty()) {
            return answers.get(answers.size() - 1);
        }
        if (!equals.isEmpty()) {
            return equals.get(equals.size() - 1);
        }
        return nonEmpty.get(nonEmpty.size() - 1);
    }

    private static String pickCompiler() throws Exception {
        for (String compiler : List.of("clang++", "g++")) {
            Process probe = new ProcessBuilder("bash", "-lc", "command -v " + compiler)
                    .redirectErrorStream(true)
                    .start();
            String out = new String(probe.getInputStream().readAllBytes());
            int rc = probe.waitFor();
            if (rc == 0 && !out.trim().isEmpty()) {
                return compiler;
            }
        }
        throw new RuntimeException("No C++ compiler found (clang++/g++).");
    }

    private static Path cppSource(Path root) {
        return root.resolve("solutionsCpp").resolve("Euler689.cpp");
    }

    private static boolean shouldSkipCheckpoints(Path root) {
        Path src = cppSource(root);
        try {
            String text = Files.readString(src);
            return text.contains("--skip-checkpoints");
        } catch (Exception ex) {
            return false;
        }
    }

    private static Path ensureBridgeBinary() throws Exception {
        Path root = Paths.get(System.getProperty("user.dir"));
        Path src = cppSource(root);
        Path bin = root.resolve("solutionsCpp").resolve(".euler689_java_bridge");

        boolean rebuild = Files.notExists(bin)
                || Files.getLastModifiedTime(src).compareTo(Files.getLastModifiedTime(bin)) > 0;

        if (rebuild) {
            String compiler = pickCompiler();
            Process compile = new ProcessBuilder(
                    compiler,
                    "-std=c++17",
                    "-O2",
                    src.toString(),
                    "-o",
                    bin.toString())
                    .inheritIO()
                    .start();
            if (compile.waitFor() != 0) {
                throw new RuntimeException("Failed to compile Euler689 C++ bridge.");
            }
        }

        return bin;
    }

    private static String runBridge(Path bin, Path root, Path srcDir) throws Exception {
        List<String> cmd = new ArrayList<>();
        cmd.add(bin.toString());
        if (shouldSkipCheckpoints(root)) {
            cmd.add("--skip-checkpoints");
        }

        Process first = new ProcessBuilder(cmd)
                .directory(root.toFile())
                .redirectErrorStream(true)
                .start();
        String out = new String(first.getInputStream().readAllBytes());
        int rc = first.waitFor();
        if (rc == 0) {
            return out;
        }

        Process second = new ProcessBuilder(cmd)
                .directory(srcDir.toFile())
                .redirectErrorStream(true)
                .start();
        String out2 = new String(second.getInputStream().readAllBytes());
        int rc2 = second.waitFor();
        if (rc2 == 0) {
            return out2;
        }

        throw new RuntimeException("Euler689 C++ bridge failed.\n" + out + "\n" + out2);
    }

    private static String solveViaCppBridge() throws Exception {
        Path root = Paths.get(System.getProperty("user.dir"));
        Path src = cppSource(root);
        Path bin = ensureBridgeBinary();
        String out = runBridge(bin, root, src.getParent());
        String parsed = parseOutput(out);
        if (parsed.isEmpty()) {
            throw new RuntimeException("Euler689 C++ bridge produced empty output.");
        }
        return parsed;
    }

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