Problem 777: Lissajous Curves

Project Euler Problem 777 Solution

EulerSolve provides an optimized solution for Project Euler Problem 777, Lissajous Curves, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For each ordered coprime pair $(a,b)$ with $2 \le a,b \le n$, the corresponding Lissajous-curve contribution collapses to a rational expression determined by $a$, $b$, and by how the factors $2$ and $5$ of $10$ are distributed between the two frequencies. The goal is to compute the total $$s(n)=\sum_{\substack{2 \le a,b \le n\\ \gcd(a,b)=1}} w(a,b)$$ exactly for very large $n$. The implementations therefore keep the integer numerator $4s(n)$ throughout the calculation and divide by $4$ only at the end. Mathematical Approach Write $$C_n=\{(a,b):2 \le a,b \le n,\ \gcd(a,b)=1\},\qquad T(k)=\frac{k(k+1)}{2}.$$ The full answer is built from two global coprime sums over $C_n$, plus a correction on a smaller exceptional set determined only by divisibility with respect to $10$. Step 1: Rewrite Coprimality with Möbius Inversion The starting point is $$\mathbf 1_{\gcd(a,b)=1}=\sum_{d \mid a,\ d \mid b}\mu(d)=\sum_{d \mid \gcd(a,b)}\mu(d).$$ This converts coprime pair sums into divisor sums....

Detailed mathematical approach

Problem Summary

For each ordered coprime pair $(a,b)$ with $2 \le a,b \le n$, the corresponding Lissajous-curve contribution collapses to a rational expression determined by $a$, $b$, and by how the factors $2$ and $5$ of $10$ are distributed between the two frequencies. The goal is to compute the total

$$s(n)=\sum_{\substack{2 \le a,b \le n\\ \gcd(a,b)=1}} w(a,b)$$

exactly for very large $n$. The implementations therefore keep the integer numerator $4s(n)$ throughout the calculation and divide by $4$ only at the end.

Mathematical Approach

Write

$$C_n=\{(a,b):2 \le a,b \le n,\ \gcd(a,b)=1\},\qquad T(k)=\frac{k(k+1)}{2}.$$

The full answer is built from two global coprime sums over $C_n$, plus a correction on a smaller exceptional set determined only by divisibility with respect to $10$.

Step 1: Rewrite Coprimality with Möbius Inversion

The starting point is

$$\mathbf 1_{\gcd(a,b)=1}=\sum_{d \mid a,\ d \mid b}\mu(d)=\sum_{d \mid \gcd(a,b)}\mu(d).$$

This converts coprime pair sums into divisor sums. Over the full square $1 \le a,b \le n$, the two basic aggregates become

$$S_{ab}^{\ast}(n)=\sum_{\substack{1 \le a,b \le n\\ \gcd(a,b)=1}}ab=\sum_{d=1}^{n}\mu(d)d^2T\left(\left\lfloor\frac{n}{d}\right\rfloor\right)^2,$$

$$S_{a}^{\ast}(n)=\sum_{\substack{1 \le a,b \le n\\ \gcd(a,b)=1}}a=\sum_{d=1}^{n}\mu(d)d\,T\left(\left\lfloor\frac{n}{d}\right\rfloor\right)\left\lfloor\frac{n}{d}\right\rfloor.$$

These are exactly the two large Möbius sums computed first.

Step 2: Remove the Boundary Cases $a=1$ or $b=1$

The target sum starts at $2$, not $1$, so the boundary contributions must be subtracted. Since

$$\sum_{m=1}^{n} m = T(n),$$

we obtain

$$S_{ab}(n)=S_{ab}^{\ast}(n)-2T(n)+1,$$

$$S_{a}(n)=S_{a}^{\ast}(n)-n-T(n)+1.$$

Hence

$$S_{ab}(n)=\sum_{(a,b)\in C_n}ab,\qquad S_{a}(n)=\sum_{(a,b)\in C_n}a.$$

Because $C_n$ is symmetric under swapping the coordinates, the same $S_a(n)$ also equals $\sum_{(a,b)\in C_n} b$.

Step 3: Isolate the Base-10 Exceptional Set

Define the bucket map

$$g(m)=\gcd(m,10)\in\{1,2,5,10\}.$$

The correction layer depends only on whether one frequency contributes the factor $2$ of $10$ and the other contributes the factor $5$. That condition is

$$E_n=\{(a,b)\in C_n:g(a)g(b)=10\}.$$

The possible bucket matches are therefore

$$1 \leftrightarrow 10,\qquad 2 \leftrightarrow 5.$$

So the entire non-generic part of the problem reduces to counting coprime pairs in four complementary residue-type buckets.

Step 4: Count the Exceptional Pairs Without Enumerating Them

For any subset $U \subseteq \{1,\dots,n\}$, define for fixed $a$

$$N_U(a)=\sum_{d \mid a}\mu(d)\,\#\{b \in U:d \mid b\},$$

$$W_U(a)=\sum_{d \mid a}\mu(d)\sum_{\substack{b \in U\\ d \mid b}} b.$$

These give the count and weighted sum of $b$-values in $U$ that are coprime to $a$. The four needed subsets are $g(b)=10$, $g(b)=5$, $g(b)=2$, and $g(b)=1$, and each one has a closed-form arithmetic progression formula.

If $g(b)=10$, then $b$ is a multiple of $10$. With $L=\operatorname{lcm}(d,10)$,

$$C_{10}(d)=\left\lfloor\frac{n}{L}\right\rfloor,\qquad Q_{10}(d)=L\,T\left(\left\lfloor\frac{n}{L}\right\rfloor\right).$$

If $g(b)=5$, then $b$ is an odd multiple of $5$. When $2 \mid d$ there are no such multiples. Otherwise, with $L=\operatorname{lcm}(d,5)$, $t=\left\lfloor n/L \right\rfloor$, and $k=\left\lceil t/2 \right\rceil$,

$$C_{5}(d)=k,\qquad Q_{5}(d)=Lk^2,$$

because $1+3+\cdots+(2k-1)=k^2$.

If $g(b)=2$, then $b$ is even but not divisible by $5$. When $5 \mid d$ the count is zero. Otherwise, with $L=\operatorname{lcm}(d,2)$, $t=\left\lfloor n/L \right\rfloor$, and $t_5=\left\lfloor n/(5L) \right\rfloor$,

$$C_{2}(d)=t-t_5,\qquad Q_{2}(d)=L\,T(t)-5L\,T(t_5).$$

If $g(b)=1$, then $b$ is coprime to $10$. When $\gcd(d,10)\ne 1$ the count is zero. Otherwise, with $x=\left\lfloor n/d \right\rfloor$, inclusion-exclusion on divisibility by $2$ and $5$ gives

$$C_{1}(d)=x-\left\lfloor\frac{x}{2}\right\rfloor-\left\lfloor\frac{x}{5}\right\rfloor+\left\lfloor\frac{x}{10}\right\rfloor,$$

$$Q_{1}(d)=d\left(T(x)-2T\left(\left\lfloor\frac{x}{2}\right\rfloor\right)-5T\left(\left\lfloor\frac{x}{5}\right\rfloor\right)+10T\left(\left\lfloor\frac{x}{10}\right\rfloor\right)\right).$$

For $g(a)=10$, the value $b=1$ must be removed afterward, because the global sum begins at $2$.

Step 5: Assemble the Exact Total

Summing the exceptional counts and weighted sums over all $a$ gives

$$N_{10}(n)=|E_n|,\qquad S_{a,10}(n)=\sum_{(a,b)\in E_n} a,\qquad S_{ab,10}(n)=\sum_{(a,b)\in E_n} ab.$$

The implementations then combine everything as

$$4s(n)=8S_{ab}(n)-12S_a(n)+6S_{a,10}(n)+4N_{10}(n)-6S_{ab,10}(n).$$

Equivalently,

$$s(n)=2S_{ab}(n)-3S_a(n)+\frac{6S_{a,10}(n)+4N_{10}(n)-6S_{ab,10}(n)}{4}.$$

This can also be read pairwise. Every $(a,b)\in C_n$ contributes

$$w(a,b)= \begin{cases} 2ab-\dfrac{3}{2}(a+b), & g(a)g(b)\ne 10,\\[6pt] \dfrac{2ab-3a-3b+4}{4}, & g(a)g(b)=10. \end{cases}$$

So the total is exactly

$$s(n)=\sum_{(a,b)\in C_n} w(a,b).$$

The quarter-factor in the exceptional case is precisely why the programs keep $4s(n)$ as an integer.

Worked Example: $n=10$

For $n=10$, the Möbius sums give

$$S_{ab}(10)=1554,\qquad S_a(10)=261.$$

The exceptional set $E_{10}$ contains $14$ ordered pairs, including $(2,5)$, $(3,10)$, $(5,8)$, $(10,3)$, and $(10,9)$. Its totals are

$$N_{10}(10)=14,\qquad S_{a,10}(10)=89,\qquad S_{ab,10}(10)=580.$$

Therefore

$$\begin{aligned} 4s(10) &=8\cdot 1554-12\cdot 261+6\cdot 89+4\cdot 14-6\cdot 580\\ &=12432-3132+534+56-3480\\ &=6410, \end{aligned}$$

$$s(10)=\frac{6410}{4}=1602.5.$$

This matches the checkpoint built into the implementation.

How the Code Works

The C++, Python, and Java implementations all follow the same mathematical decomposition. First they build a Möbius table up to $n$. Next they evaluate the two large divisor sums for $S_{ab}^{\ast}(n)$ and $S_a^{\ast}(n)$, then apply the boundary subtraction that removes the cases with $a=1$ or $b=1$.

After that, the implementation processes the four base-10 buckets. Instead of factoring every $a$ separately and enumerating its divisors, it loops over each divisor $d$, computes the relevant count and weighted-sum formulas once, and distributes the Möbius-weighted contribution to every multiple of $d$. This produces per-$a$ tables for the exceptional set in near-harmonic time. The C++ version can run the four independent bucket passes in parallel; the Python and Java versions preserve the same arithmetic structure and final exact value.

Complexity Analysis

The Möbius sieve is linear, so it costs $O(n)$ time and $O(n)$ memory. The dominant work is the divisor-to-multiples sweep used by the four bucket tables, whose total size is controlled by

$$\sum_{d=1}^{n}\left\lfloor\frac{n}{d}\right\rfloor = O(n\log n).$$

Therefore the overall algorithm runs in $O(n\log n)$ time and uses $O(n)$ memory. Parallel execution improves wall-clock time in the C++ implementation but does not change the asymptotic bound.

Footnotes and References

Problem page: https://projecteuler.net/problem=777
Lissajous curves: Wikipedia - Lissajous curve
Mobius inversion: Wikipedia - Mobius inversion formula
Mobius function: Wikipedia - Mobius function
Inclusion-exclusion principle: Wikipedia - Inclusion-exclusion principle
Coprime integers: Wikipedia - Coprime integers

Problem 777 source code

C++

#include <algorithm>
#include <cstdint>
#include <iomanip>
#include <iostream>
#include <numeric>
#include <sstream>
#include <string>
#include <thread>
#include <vector>

using namespace std;

namespace {

using i128 = __int128_t;

inline long long sum_upto(long long n) {
    return static_cast<long long>((static_cast<i128>(n) * (n + 1)) / 2);
}

vector<int> mobius_sieve(int n) {
    vector<int> mu(n + 1, 1);
    vector<int> primes;
    vector<bool> composite(n + 1, false);
    mu[0] = 0;
    for (int i = 2; i <= n; ++i) {
        if (!composite[i]) {
            primes.push_back(i);
            mu[i] = -1;
        }
        for (int p : primes) {
            long long v = static_cast<long long>(i) * p;
            if (v > n) break;
            composite[static_cast<int>(v)] = true;
            if (i % p == 0) {
                mu[static_cast<int>(v)] = 0;
                break;
            }
            mu[static_cast<int>(v)] = -mu[i];
        }
    }
    return mu;
}

enum class SubsetType {
    S10,
    S5Odd,
    S2Not5,
    SCoprime10,
};

inline void count_sum_S10(int n, int d, long long& cnt, long long& sum) {
    int g = 1;
    if (d % 10 == 0) g = 10;
    else if (d % 5 == 0) g = 5;
    else if (d % 2 == 0) g = 2;
    long long L = static_cast<long long>(d / g) * 10LL;
    long long t = n / L;
    cnt = t;
    sum = static_cast<long long>((static_cast<i128>(L) * t * (t + 1)) / 2);
}

inline void count_sum_S5Odd(int n, int d, long long& cnt, long long& sum) {
    if (d % 2 == 0) {
        cnt = 0;
        sum = 0;
        return;
    }
    int g = (d % 5 == 0) ? 5 : 1;
    long long L = static_cast<long long>(d / g) * 5LL;
    long long t = n / L;
    long long k = (t + 1) / 2;
    cnt = k;
    sum = L * k * k;
}

inline void count_sum_S2Not5(int n, int d, long long& cnt, long long& sum) {
    if (d % 5 == 0) {
        cnt = 0;
        sum = 0;
        return;
    }
    long long L = (d % 2 == 0) ? d : 2LL * d;
    long long t = n / L;
    long long t5 = n / (5LL * L);
    cnt = t - t5;
    long long sumL = static_cast<long long>((static_cast<i128>(L) * t * (t + 1)) / 2);
    long long sum5L = static_cast<long long>((static_cast<i128>(5LL * L) * t5 * (t5 + 1)) / 2);
    sum = sumL - sum5L;
}

inline void count_sum_SCoprime10(int n, int d, long long& cnt, long long& sum) {
    if (d % 2 == 0 || d % 5 == 0) {
        cnt = 0;
        sum = 0;
        return;
    }
    long long x = n / d;
    long long t2 = x / 2;
    long long t5 = x / 5;
    long long t10 = x / 10;
    cnt = x - t2 - t5 + t10;
    long long sum_all = sum_upto(x);
    long long sum2 = static_cast<long long>(2LL * sum_upto(t2));
    long long sum5 = static_cast<long long>(5LL * sum_upto(t5));
    long long sum10 = static_cast<long long>(10LL * sum_upto(t10));
    long long sum_t = sum_all - sum2 - sum5 + sum10;
    sum = static_cast<long long>(static_cast<i128>(d) * sum_t);
}

void compute_subset(int n, const vector<int>& mu, SubsetType type,
                    vector<long long>& cnt, vector<long long>& sum) {
    for (int d = 1; d <= n; ++d) {
        int mud = mu[d];
        if (mud == 0) continue;
        long long c = 0;
        long long s = 0;
        switch (type) {
            case SubsetType::S10:
                count_sum_S10(n, d, c, s);
                break;
            case SubsetType::S5Odd:
                count_sum_S5Odd(n, d, c, s);
                break;
            case SubsetType::S2Not5:
                count_sum_S2Not5(n, d, c, s);
                break;
            case SubsetType::SCoprime10:
                count_sum_SCoprime10(n, d, c, s);
                break;
        }
        if (c == 0) continue;
        for (int a = d; a <= n; a += d) {
            cnt[a] += static_cast<long long>(mud) * c;
            sum[a] += static_cast<long long>(mud) * s;
        }
    }
}

inline int gcd10_bucket(int a) {
    if (a % 10 == 0) return 10;
    if (a % 5 == 0) return 5;
    if (a % 2 == 0) return 2;
    return 1;
}

i128 compute_s_num(int n) {
    vector<int> mu = mobius_sieve(n);

    i128 sum_ab = 0;
    i128 sum_a = 0;
    for (int d = 1; d <= n; ++d) {
        int mud = mu[d];
        if (mud == 0) continue;
        long long k = n / d;
        long long s = sum_upto(k);
        i128 dd = d;
        i128 ss = s;
        sum_ab += static_cast<i128>(mud) * dd * dd * ss * ss;
        sum_a += static_cast<i128>(mud) * dd * ss * k;
    }

    long long s_all = sum_upto(n);
    i128 sum_ab_2n = sum_ab - 2 * static_cast<i128>(s_all) + 1;
    i128 sum_a_2n = sum_a - n - static_cast<i128>(s_all) + 1;

    vector<long long> cnt10(n + 1, 0), sum10(n + 1, 0);
    vector<long long> cnt5o(n + 1, 0), sum5o(n + 1, 0);
    vector<long long> cnt2n(n + 1, 0), sum2n(n + 1, 0);
    vector<long long> cntg1(n + 1, 0), sumg1(n + 1, 0);

    unsigned hw = thread::hardware_concurrency();
    bool use_threads = hw >= 2 && n >= 2000;
    if (use_threads) {
        thread t1(compute_subset, n, cref(mu), SubsetType::S10, ref(cnt10), ref(sum10));
        thread t2(compute_subset, n, cref(mu), SubsetType::S5Odd, ref(cnt5o), ref(sum5o));
        thread t3(compute_subset, n, cref(mu), SubsetType::S2Not5, ref(cnt2n), ref(sum2n));
        compute_subset(n, mu, SubsetType::SCoprime10, cntg1, sumg1);
        t1.join();
        t2.join();
        t3.join();
    } else {
        compute_subset(n, mu, SubsetType::S10, cnt10, sum10);
        compute_subset(n, mu, SubsetType::S5Odd, cnt5o, sum5o);
        compute_subset(n, mu, SubsetType::S2Not5, cnt2n, sum2n);
        compute_subset(n, mu, SubsetType::SCoprime10, cntg1, sumg1);
    }

    i128 count10 = 0;
    i128 sum_a10 = 0;
    i128 sum_ab10 = 0;
    for (int a = 2; a <= n; ++a) {
        long long cnt = 0;
        long long s = 0;
        int g = gcd10_bucket(a);
        if (g == 1) {
            cnt = cnt10[a];
            s = sum10[a];
        } else if (g == 2) {
            cnt = cnt5o[a];
            s = sum5o[a];
        } else if (g == 5) {
            cnt = cnt2n[a];
            s = sum2n[a];
        } else {
            cnt = cntg1[a];
            s = sumg1[a];
            cnt -= 1;  // exclude b=1
            s -= 1;
        }
        count10 += cnt;
        sum_a10 += static_cast<i128>(a) * cnt;
        sum_ab10 += static_cast<i128>(a) * s;
    }

    i128 base = 2 * sum_ab_2n - 3 * sum_a_2n;
    i128 corr = 6 * sum_a10 + 4 * count10 - 6 * sum_ab10;
    i128 total_num = 4 * base + corr;
    return total_num;
}

string format_scientific(long double value) {
    ostringstream oss;
    oss.setf(ios::scientific);
    oss << setprecision(9) << value;
    string s = oss.str();
    size_t pos = s.find('e');
    if (pos == string::npos) return s;
    char sign = s[pos + 1];
    size_t i = pos + 2;
    while (i < s.size() && s[i] == '0') ++i;
    string exp = (i < s.size()) ? s.substr(i) : "0";
    string out = s.substr(0, pos + 1);
    if (sign == '-') out.push_back('-');
    out += exp;
    return out;
}

}  // namespace

int main(int argc, char** argv) {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);

    int n = 1'000'000;
    if (argc >= 2) n = stoi(argv[1]);

    const i128 expect10 = 6410;      // 1602.5 * 4
    const i128 expect100 = 97026020; // 24256505 * 4

    i128 val10 = compute_s_num(10);
    if (val10 != expect10) {
        cerr << "Validation failed: s(10) expected 1602.5, got "
             << static_cast<long double>(val10) / 4.0L << "\n";
        return 1;
    }
    i128 val100 = compute_s_num(100);
    if (val100 != expect100) {
        cerr << "Validation failed: s(100) expected 24256505, got "
             << static_cast<long double>(val100) / 4.0L << "\n";
        return 1;
    }

    i128 total_num = compute_s_num(n);
    long double value = static_cast<long double>(total_num) / 4.0L;
    cout << format_scientific(value) << "\n";
    return 0;
}

Python

def solve():
    n = 1000000

    def sum_upto(n): return n * (n+1) // 2

    # Mobius sieve
    mu = [0]*(n+1); mu[1] = 1; comp = bytearray(n+1); primes = []
    for i in range(2, n+1):
        if not comp[i]: primes.append(i); mu[i] = -1
        for p in primes:
            if i*p > n: break
            comp[i*p] = 1
            if i % p == 0: mu[i*p] = 0; break
            mu[i*p] = -mu[i]

    # Compute sum_ab and sum_a
    sum_ab = 0; sum_a_val = 0
    for d in range(1, n+1):
        md = mu[d]
        if md == 0: continue
        k = n // d; s = sum_upto(k)
        sum_ab += md * d * d * s * s
        sum_a_val += md * d * s * k
    s_all = sum_upto(n)
    sum_ab_2n = sum_ab - 2*s_all + 1
    sum_a_2n = sum_a_val - n - s_all + 1

    # Subset counting
    cnt10 = [0]*(n+1); sum10 = [0]*(n+1)
    cnt5o = [0]*(n+1); sum5o = [0]*(n+1)
    cnt2n = [0]*(n+1); sum2n = [0]*(n+1)
    cntg1 = [0]*(n+1); sumg1 = [0]*(n+1)

    for d in range(1, n+1):
        md = mu[d]
        if md == 0: continue
        # S10
        g = 1
        if d%10==0: g=10
        elif d%5==0: g=5
        elif d%2==0: g=2
        L = d//g * 10; t = n//L
        if t > 0:
            c10 = t; s10_v = L*t*(t+1)//2
            for a in range(d, n+1, d): cnt10[a] += md*c10; sum10[a] += md*s10_v
        # S5Odd
        if d%2 != 0:
            g5 = 5 if d%5==0 else 1
            L5 = d//g5 * 5; t5 = n//L5; k5 = (t5+1)//2
            if k5 > 0:
                c5o = k5; s5o_v = L5*k5*k5
                for a in range(d, n+1, d): cnt5o[a] += md*c5o; sum5o[a] += md*s5o_v
        # S2Not5
        if d%5 != 0:
            L2 = d if d%2==0 else 2*d; t2 = n//L2; t25 = n//(5*L2)
            c2n = t2-t25
            if c2n > 0:
                s2n_v = L2*t2*(t2+1)//2 - 5*L2*t25*(t25+1)//2
                for a in range(d, n+1, d): cnt2n[a] += md*c2n; sum2n[a] += md*s2n_v
        # SCoprime10
        if d%2 != 0 and d%5 != 0:
            x = n//d; t2c = x//2; t5c = x//5; t10c = x//10
            cg = x - t2c - t5c + t10c
            if cg > 0:
                sg = (sum_upto(x) - 2*sum_upto(t2c) - 5*sum_upto(t5c) + 10*sum_upto(t10c))
                sg_v = d * sg
                for a in range(d, n+1, d): cntg1[a] += md*cg; sumg1[a] += md*sg_v

    count10_total = 0; sum_a10 = 0; sum_ab10 = 0
    for a in range(2, n+1):
        g = 1
        if a%10==0: g=10
        elif a%5==0: g=5
        elif a%2==0: g=2
        if g==1: cnt=cnt10[a]; s=sum10[a]
        elif g==2: cnt=cnt5o[a]; s=sum5o[a]
        elif g==5: cnt=cnt2n[a]; s=sum2n[a]
        else: cnt=cntg1[a]-1; s=sumg1[a]-1
        count10_total += cnt; sum_a10 += a*cnt; sum_ab10 += a*s

    base = 2*sum_ab_2n - 3*sum_a_2n
    corr = 6*sum_a10 + 4*count10_total - 6*sum_ab10
    total_num = 4*base + corr
    value = total_num / 4.0

    # Format in scientific notation
    import math
    if value == 0: return "0.000000000e0"
    exp = int(math.floor(math.log10(abs(value))))
    mantissa = value / 10**exp
    return f"{mantissa:.9f}e{exp}"

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

Java

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

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

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