Problem 591: Best Approximations by Quadratic Integers

View on Project Euler

Project Euler Problem 591 Solution

EulerSolve provides an optimized solution for Project Euler Problem 591, Best Approximations by Quadratic Integers, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For each non-square integer \(d \lt 100\), we approximate \(\pi\) by numbers of the form $$a+b\sqrt d,\qquad |a|,|b|\le 10^{13}.$$ Among all such pairs \((a,b)\), we choose the one that minimizes $$\left|\pi-\left(a+b\sqrt d\right)\right|.$$ The quantity called the integral part is the coefficient \(a\), not the floor of the real number \(a+b\sqrt d\). The required output is therefore $$\sum_{d\in S} |a_d|,$$ where \(S\) is the set of non-square integers from \(2\) to \(99\), and \(a_d\) is the integer coefficient belonging to the best approximation for that specific \(d\). Mathematical Approach Write \(\alpha=\sqrt d\) and \(\tau=\pi-3\). The central observation is that once \(b\) is fixed, the optimal \(a\) is forced to be the nearest integer to \(\pi-b\alpha\). So the real difficulty is to find the right \(b\) without scanning all \(2\cdot 10^{13}+1\) possibilities. Step 1: Reduce the problem to a nearest-integer search For fixed \(b\), define $$x_b=\pi-b\alpha.$$ Then minimizing \(\left|\pi-(a+b\alpha)\right|\) over all admissible \(a\) is the same as minimizing $$|x_b-a|$$ over integers \(a\) with \(|a|\le 10^{13}\)....

Detailed mathematical approach

Problem Summary

For each non-square integer \(d \lt 100\), we approximate \(\pi\) by numbers of the form

$$a+b\sqrt d,\qquad |a|,|b|\le 10^{13}.$$

Among all such pairs \((a,b)\), we choose the one that minimizes

$$\left|\pi-\left(a+b\sqrt d\right)\right|.$$

The quantity called the integral part is the coefficient \(a\), not the floor of the real number \(a+b\sqrt d\). The required output is therefore

$$\sum_{d\in S} |a_d|,$$

where \(S\) is the set of non-square integers from \(2\) to \(99\), and \(a_d\) is the integer coefficient belonging to the best approximation for that specific \(d\).

Mathematical Approach

Write \(\alpha=\sqrt d\) and \(\tau=\pi-3\). The central observation is that once \(b\) is fixed, the optimal \(a\) is forced to be the nearest integer to \(\pi-b\alpha\). So the real difficulty is to find the right \(b\) without scanning all \(2\cdot 10^{13}+1\) possibilities.

Step 1: Reduce the problem to a nearest-integer search

For fixed \(b\), define

$$x_b=\pi-b\alpha.$$

Then minimizing \(\left|\pi-(a+b\alpha)\right|\) over all admissible \(a\) is the same as minimizing

$$|x_b-a|$$

over integers \(a\) with \(|a|\le 10^{13}\). Therefore the best \(a\) is the nearest integer to \(x_b\), and the two-variable problem becomes

$$\min_{|b|\le 10^{13}} \left\| \pi-b\alpha \right\|_{\mathbb Z},$$

where \(\|y\|_{\mathbb Z}\) means the distance from \(y\) to the nearest integer.

There is also a strong bound on \(b\). The trivial candidate \(a=3\), \(b=0\) has error \(\tau=\pi-3\lt 1\), so any true optimum must have error below \(1\). If \(|b|\alpha \gt 10^{13}+\pi+1\), then for every admissible \(a\),

$$\left|a+b\alpha-\pi\right|\ge |b|\alpha-|a|-\pi \gt 1,$$

which cannot be optimal. The implementation therefore uses the safe bound

$$|b|\le B,\qquad B=\min\!\left(10^{13},\left\lfloor\frac{10^{13}+\pi+2}{\alpha}\right\rfloor\right).$$

Step 2: Use convergents of \(\sqrt d\)

Because \(\alpha=\sqrt d\) is a quadratic irrational, its continued fraction is periodic, and its convergents \(p/q\) give exceptionally good rational approximations. For each convergent we write

$$q\alpha=p+\delta,$$

so \(\delta=q\alpha-p\) is very small.

This matters because changing \(b\) by \(q\) changes \(b\alpha\) by almost the integer \(p\):

$$ (b+q)\alpha=b\alpha+p+\delta.$$

So within the arithmetic progression \(b=b_0+tq\), the fractional part of \(b\alpha\) drifts only by the tiny amount \(t\delta\). That turns a huge search over \(b\) into a very small local search around each convergent.

Step 3: Split the search into residue classes modulo \(q\)

Since every convergent satisfies \(\gcd(p,q)=1\), the inverse of \(p\) modulo \(q\) exists. For each residue \(j\in\{0,1,\dots,q-1\}\), choose \(b_0\) so that

$$b_0p\equiv j \pmod q.$$

Then \(b_0p=j+kq\) for some integer \(k\). If we now write

$$b=b_0+tq,$$

we get

$$\begin{aligned} b\alpha &=b_0\alpha+tq\alpha \\ &=\frac{b_0(p+\delta)}{q}+t(p+\delta) \\ &=k+tp+\frac{j}{q}+\frac{b_0\delta}{q}+t\delta. \end{aligned}$$

Since \(\pi=3+\tau\), it follows that

$$\pi-b\alpha=\left(3-k-tp\right)+\left(\tau-\frac{j}{q}-\frac{b_0\delta}{q}-t\delta\right).$$

The first parenthesis is an integer, so the quality of the approximation is controlled entirely by the small residual

$$\operatorname{diff}(t)=\frac{j}{q}+\frac{b_0\delta}{q}-\tau+t\delta.$$

We want this quantity to be as close to \(0\) as possible.

Step 4: Locate the best \(t\) from a linear estimate

For fixed \(q\), \(p\), and residue class \(j\), the expression \(\operatorname{diff}(t)\) is linear in \(t\). Therefore the best \(t\) is near the real minimizer

$$t\approx -\frac{\frac{j}{q}+\frac{b_0\delta}{q}-\tau}{\delta}.$$

Only residue classes with \(j/q\) near \(\tau\) can compete, because \(|b_0\delta/q|\le |\delta|\) is tiny for a useful convergent. The implementations therefore inspect only a small fixed window of \(j\)-values around \(\tau q\), and for each such class they inspect only a small fixed window of \(t\)-values around the estimate above.

The admissible \(t\)-range is also clipped by \(|b|\le B\). Its endpoints are checked explicitly so that a best solution on the boundary is never missed.

Step 5: Recover \(a\) and compare candidates safely

Once a candidate \(b\) is known, we compute

$$x_b=\pi-b\alpha$$

and take the nearest integer \(a_0\) to \(x_b\). Because the computation uses finite-precision arithmetic, the implementations test the three integers

$$a_0-1,\qquad a_0,\qquad a_0+1$$

instead of trusting one rounding step blindly. This guarantees that the true nearest admissible integer is captured even if \(x_b\) lies very close to a half-integer.

The best pair is the one with smallest absolute error. If two candidates are numerically tied, the implementations keep the one with smaller \(|a|\), and then smaller \(|b|\), to make the output deterministic.

Worked Example: \(d=2\), \(N=10\)

Here \(\alpha=\sqrt 2\approx 1.414213562\) and \(\tau=\pi-3\approx 0.141592654\).

Take the convergent

$$\frac{p}{q}=\frac{3}{2},$$

for which

$$\delta=2\sqrt 2-3\approx -0.171572875.$$

Since \(\tau q\approx 0.283\), the most relevant residue is \(j=0\). Because \(p\equiv 1\pmod 2\), the corresponding base residue is \(b_0\equiv 0\pmod 2\), so we take \(b_0=0\).

Then

$$\operatorname{diff}(t)=-\tau+t\delta.$$

The linear estimate gives

$$t\approx -\frac{-\tau}{\delta}\approx -0.825,$$

so \(t=-1\) is the natural nearby candidate. This yields

$$b=b_0+tq=-2.$$

Now compute

$$\pi-b\sqrt 2=\pi+2\sqrt 2\approx 5.969019778,$$

whose nearest integer is \(a=6\). Hence the approximation is

$$6-2\sqrt 2\approx 3.171572875,$$

with error about

$$|\,\pi-(6-2\sqrt 2)\,|\approx 0.029980221.$$

This matches the hard-coded sanity check used by the implementations.

How the Code Works

The C++, Python, and Java implementations all follow the same mathematical search. For each non-square \(d\lt 100\), they compute \(\sqrt d\) with high precision, derive the safe bound \(B\), and generate continued-fraction convergents until the denominator exceeds that bound.

For every convergent, the implementation forms the small error term \(\delta=q\sqrt d-p\), looks only at residue classes \(j\) near \((\pi-3)q\), converts each class into a base value \(b_0\) by modular inversion, and derives the permitted interval of \(t\)-values from the constraint \(|b|\le B\).

Instead of scanning every possible \(t\), the implementation evaluates only a short interval around the linear estimate for the best \(t\), plus the interval endpoints. Each resulting \(b\) is converted into a nearby integer coefficient \(a\), and the error is checked explicitly.

The C++ version distributes the independent \(d\)-values across worker threads. The Python version performs the same search directly with decimal arithmetic. The Java version delegates execution to the same compiled numeric core, so its output agrees with the C++ behavior.

Complexity Analysis

For fixed \(d\), the denominators of the convergents of a quadratic irrational grow exponentially, so only \(O(\log B)\) convergents are needed before the denominator passes the search bound. The inspected windows in \(j\) and \(t\) are fixed-size constants, so each convergent contributes only \(O(1)\) candidate evaluations.

Therefore the practical running time per \(d\) is \(O(\log B)\), which here is effectively \(O(\log N)\), and the memory usage is \(O(\log B)\) if the convergents are stored explicitly. Since there are only finitely many non-square values \(d\lt 100\), the total cost is this per-\(d\) work multiplied by a small constant, with straightforward parallelization in C++.

Footnotes and References

  1. Problem page: Project Euler 591
  2. Continued fractions: Wikipedia - Continued fraction
  3. Convergents of continued fractions: Wikipedia - Convergent
  4. Quadratic irrationals: Wikipedia - Quadratic irrational
  5. Diophantine approximation: Wikipedia - Diophantine approximation

Problem 591 source code

C++

#include <algorithm>
#include <atomic>
#include <cmath>
#include <cstdint>
#include <cstdlib>
#include <iostream>
#include <string>
#include <thread>
#include <vector>

#include <boost/math/constants/constants.hpp>
#include <boost/multiprecision/cpp_dec_float.hpp>

using std::cerr;
using std::cout;
using std::string;
using boost::multiprecision::cpp_dec_float_100;

using Big = cpp_dec_float_100;

struct Approx {
    long long a = 0;
    long long b = 0;
    Big err = 0;
    bool valid = false;
};

static inline long long llabs_ll(long long x) { return x < 0 ? -x : x; }

static long long floor_div(long long a, long long b) {
    __int128 A = static_cast<__int128>(a);
    __int128 B = static_cast<__int128>(b);
    __int128 q = A / B;
    __int128 r = A % B;
    if (r != 0 && A < 0) q -= 1;
    return static_cast<long long>(q);
}

static long long ceil_div(long long a, long long b) {
    __int128 A = static_cast<__int128>(a);
    __int128 B = static_cast<__int128>(b);
    __int128 q = A / B;
    __int128 r = A % B;
    if (r != 0 && A > 0) q += 1;
    return static_cast<long long>(q);
}

static long long floor_ll(const Big& x) {
    long long t = x.convert_to<long long>();
    if (x < 0 && Big(t) != x) return t - 1;
    return t;
}

static long long ceil_ll(const Big& x) {
    long long t = x.convert_to<long long>();
    if (x > 0 && Big(t) != x) return t + 1;
    return t;
}

static long long round_ll(const Big& x) {
    static const Big half("0.5");
    if (x >= 0) return floor_ll(x + half);
    return ceil_ll(x - half);
}

static long long modinv(long long a, long long mod) {
    long long b = mod;
    __int128 x0 = 1, x1 = 0;
    __int128 y0 = 0, y1 = 1;
    long long aa = a;
    while (b != 0) {
        long long q = aa / b;
        long long r = aa % b;
        aa = b;
        b = r;
        __int128 nx = x0 - static_cast<__int128>(q) * x1;
        x0 = x1;
        x1 = nx;
        __int128 ny = y0 - static_cast<__int128>(q) * y1;
        y0 = y1;
        y1 = ny;
    }
    if (aa != 1 && aa != -1) return -1;
    __int128 inv = x0;
    long long res = static_cast<long long>(inv % mod);
    if (res < 0) res += mod;
    return res;
}

static bool is_square(int n) {
    int r = static_cast<int>(std::sqrt(static_cast<double>(n)));
    return r * r == n;
}

static std::vector<int> cf_period_sqrt(int d) {
    int a0 = static_cast<int>(std::sqrt(static_cast<double>(d)));
    if (a0 * a0 == d) return {};
    long long m = 0;
    long long denom = 1;
    long long a = a0;
    std::vector<int> period;
    while (true) {
        m = denom * a - m;
        denom = (d - m * m) / denom;
        a = (a0 + m) / denom;
        period.push_back(static_cast<int>(a));
        if (denom == 1 && m == a0) break;
    }
    return period;
}

static std::vector<std::pair<long long, long long>> convergents_sqrt(int d, long long q_limit) {
    int a0 = static_cast<int>(std::sqrt(static_cast<double>(d)));
    std::vector<int> period = cf_period_sqrt(d);
    std::vector<std::pair<long long, long long>> conv;

    long long p_m2 = 1, q_m2 = 0;
    long long p_m1 = a0, q_m1 = 1;
    conv.push_back({p_m1, q_m1});

    if (period.empty()) return conv;

    size_t i = 0;
    while (true) {
        long long a = period[i % period.size()];
        ++i;
        __int128 p = static_cast<__int128>(a) * p_m1 + p_m2;
        __int128 q = static_cast<__int128>(a) * q_m1 + q_m2;
        long long pn = static_cast<long long>(p);
        long long qn = static_cast<long long>(q);
        conv.push_back({pn, qn});
        if (qn > q_limit) break;
        p_m2 = p_m1;
        q_m2 = q_m1;
        p_m1 = pn;
        q_m1 = qn;
    }
    return conv;
}

static Approx best_BQA_for_d(int d, long long N, int window_j = 25, int window_t = 25) {
    const Big PI = boost::math::constants::pi<Big>();
    const Big TAU = PI - Big(3);

    Big alpha = sqrt(Big(d));

    long long B = floor_ll((Big(N) + PI + Big(2)) / alpha);
    if (B > N) B = N;

    auto conv = convergents_sqrt(d, B);

    Approx best;
    best.err = Big("1e1000");

    auto eval_b = [&](long long b) {
        if (llabs_ll(b) > B || llabs_ll(b) > N) return;
        Big val = alpha * Big(b);
        Big a_real = PI - val;
        long long a0 = round_ll(a_real);

        for (long long da = -1; da <= 1; ++da) {
            long long a = a0 + da;
            if (llabs_ll(a) > N) continue;
            Big err = abs(a_real - Big(a));
            if (!best.valid || err < best.err ||
                (err == best.err && llabs_ll(a) < llabs_ll(best.a)) ||
                (err == best.err && llabs_ll(a) == llabs_ll(best.a) && llabs_ll(b) < llabs_ll(best.b))) {
                best.valid = true;
                best.a = a;
                best.b = b;
                best.err = err;
            }
        }
    };

    for (auto [p_raw, q_raw] : conv) {
        long long q = q_raw;
        long long p = p_raw;
        if (q <= 0) continue;

        long long inv = modinv((p % q + q) % q, q);
        if (inv < 0) continue;

        Big qB = Big(q);
        long long r = (TAU * qB).convert_to<long long>();

        long long j_lo = std::max(0LL, r - static_cast<long long>(window_j));
        long long j_hi = std::min(q - 1, r + static_cast<long long>(window_j));

        Big delta = qB * alpha - Big(p);
        if (delta == 0) continue;

        for (long long j = j_lo; j <= j_hi; ++j) {
            long long b0 = static_cast<long long>((static_cast<__int128>(j) * inv) % q);

            long long t_min = ceil_div(-B - b0, q);
            long long t_max = floor_div(B - b0, q);
            if (t_min > t_max) continue;

            Big diff = Big(j) / qB + (Big(b0) * delta) / qB - TAU;

            Big ratio = (-diff) / delta;
            long long t0 = round_ll(ratio);

            long long tt_lo = std::max(t_min, t0 - static_cast<long long>(window_t));
            long long tt_hi = std::min(t_max, t0 + static_cast<long long>(window_t));
            for (long long t = tt_lo; t <= tt_hi; ++t) {
                long long b = b0 + t * q;
                eval_b(b);
            }

            eval_b(b0 + t_min * q);
            eval_b(b0 + t_max * q);
        }
    }

    return best;
}

static void validation_checkpoints() {
    const long long N13 = 10000000000000LL;

    {
        auto r = best_BQA_for_d(2, 10);
        if (!(r.valid && r.a == 6 && r.b == -2)) {
            cerr << "Validation failed: BQA_2(pi,10) expected 6-2*sqrt(2), got a=" << r.a
                 << " b=" << r.b << "\n";
            std::exit(1);
        }
    }
    {
        auto r = best_BQA_for_d(5, 100);
        if (!(r.valid && r.a == -55 && r.b == 26)) {
            cerr << "Validation failed: BQA_5(pi,100) expected 26*sqrt(5)-55, got a=" << r.a
                 << " b=" << r.b << "\n";
            std::exit(1);
        }
    }
    {
        auto r = best_BQA_for_d(7, 1000000);
        if (!(r.valid && r.a == 560323 && r.b == -211781)) {
            cerr << "Validation failed: BQA_7(pi,1e6) expected 560323-211781*sqrt(7), got a="
                 << r.a << " b=" << r.b << "\n";
            std::exit(1);
        }
    }
    {
        auto r = best_BQA_for_d(2, N13);
        if (!(r.valid && r.a == -6188084046055LL)) {
            cerr << "Validation failed: I_2(BQA_2(pi,1e13)) expected -6188084046055, got a="
                 << r.a << "\n";
            std::exit(1);
        }
    }

    cerr << "All validation checkpoints passed.\n";
}

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

    validation_checkpoints();

    const long long N = 10000000000000LL;

    std::vector<int> ds;
    for (int d = 2; d < 100; ++d) {
        if (is_square(d)) continue;
        ds.push_back(d);
    }

    std::vector<long long> absAs(ds.size(), 0);

    unsigned numThreads = std::thread::hardware_concurrency();
    if (numThreads == 0) numThreads = 4;

    std::atomic<size_t> idx{0};
    std::vector<std::thread> threads;
    threads.reserve(numThreads);

    for (unsigned t = 0; t < numThreads; ++t) {
        threads.emplace_back([&]() {
            while (true) {
                size_t i = idx.fetch_add(1, std::memory_order_relaxed);
                if (i >= ds.size()) break;
                int d = ds[i];
                auto r = best_BQA_for_d(d, N);
                if (!r.valid) {
                    cerr << "No valid approximation found for d=" << d << "\n";
                    std::exit(1);
                }
                absAs[i] = llabs_ll(r.a);
            }
        });
    }

    for (auto& th : threads) th.join();

    __int128 sum = 0;
    for (auto v : absAs) sum += v;

    auto print_int128 = [](__int128 x) {
        if (x == 0) {
            cout << '0';
            return;
        }
        if (x < 0) {
            cout << '-';
            x = -x;
        }
        string s;
        while (x > 0) {
            int digit = static_cast<int>(x % 10);
            s.push_back(static_cast<char>('0' + digit));
            x /= 10;
        }
        std::reverse(s.begin(), s.end());
        cout << s;
    };

    print_int128(sum);
    cout << "\n";
    return 0;
}

Python

from decimal import Decimal, getcontext
import math

def solve():
    getcontext().prec = 120
    PI = Decimal('3.14159265358979323846264338327950288419716939937510'
                 '58209749445923078164062862089986280348253421170679')
    TAU = PI - 3
    N = 10**13

    def is_sq(n): r=int(math.isqrt(n)); return r*r==n

    def cf_period(d):
        a0=int(math.isqrt(d))
        if a0*a0==d: return []
        m=0; den=1; a=a0; period=[]
        while True:
            m=den*a-m; den=(d-m*m)//den; a=(a0+m)//den
            period.append(a)
            if den==1 and m==a0: break
        return period

    def convergents(d, qlim):
        a0=int(math.isqrt(d)); period=cf_period(d)
        conv=[(a0,1)]
        if not period: return conv
        pm2,qm2=1,0; pm1,qm1=a0,1; i=0
        while True:
            a=period[i%len(period)]; i+=1
            p=a*pm1+pm2; q=a*qm1+qm2
            conv.append((p,q))
            if q>qlim: break
            pm2,qm2=pm1,qm1; pm1,qm1=p,q
        return conv

    def modinv(a,m):
        b=m; x0,x1=1,0; aa=a
        while b:
            q=aa//b; aa,b=b,aa-q*b; x0,x1=x1,x0-q*x1
        if aa!=1 and aa!=-1: return -1
        return x0%m

    def best_bqa(d, NN):
        alpha=Decimal(d).sqrt()
        B=int((Decimal(NN)+PI+2)/alpha)
        if B>NN: B=NN
        conv=convergents(d,B)
        best_err=Decimal('1e1000'); best_a=0; best_b=0; found=False

        def eval_b(b):
            nonlocal best_err,best_a,best_b,found
            if abs(b)>B or abs(b)>NN: return
            val=alpha*Decimal(b); ar=PI-val; a0=int(ar+Decimal('0.5')) if ar>=0 else -int(-ar+Decimal('0.5'))
            for da in (-1,0,1):
                a=a0+da
                if abs(a)>NN: continue
                err=abs(ar-Decimal(a))
                if not found or err<best_err or (err==best_err and abs(a)<abs(best_a)):
                    found=True; best_a=a; best_b=b; best_err=err

        for p,q in conv:
            if q<=0: continue
            inv=modinv(p%q if p%q>=0 else p%q+q, q)
            if inv<0: continue
            qD=Decimal(q); r=int(TAU*qD)
            delta=qD*alpha-Decimal(p)
            if delta==0: continue
            jlo=max(0,r-25); jhi=min(q-1,r+25)
            for j in range(jlo,jhi+1):
                b0=(j*inv)%q
                tmin=(-B-b0+q-1)//q if -B-b0>=0 else -(-(- B-b0))//q
                # Proper ceil/floor div
                num=-B-b0
                tmin = num//q if num>=0 else -((-num)//q) if (-num)%q==0 else -((-num)//q+1) if num<0 else num//q
                tmin = math.ceil((-B-b0)/q) if q>0 else 0
                tmax = math.floor((B-b0)/q) if q>0 else 0
                if tmin>tmax: continue
                diff=Decimal(j)/qD+Decimal(b0)*delta/qD-TAU
                ratio=-diff/delta; t0=int(ratio+Decimal('0.5')) if ratio>=0 else -int(-ratio+Decimal('0.5'))
                ttlo=max(tmin,t0-25); tthi=min(tmax,t0+25)
                for t in range(ttlo,tthi+1):
                    eval_b(b0+t*q)
                eval_b(b0+tmin*q); eval_b(b0+tmax*q)
        return abs(best_a)

    ds = [d for d in range(2,100) if not is_sq(d)]
    total = sum(best_bqa(d,N) for d in ds)
    return str(total)

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

Java

import java.nio.file.*;
import java.util.*;
import java.util.regex.*;

public class Euler591 {
    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("Euler591.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(".euler591_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 Euler591 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("Euler591 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("Euler591 C++ bridge produced empty output.");
        }
        return parsed;
    }

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