Problem 644: Squares on the Line

Project Euler Problem 644 Solution

EulerSolve provides an optimized solution for Project Euler Problem 644, Squares on the Line, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary A position is a line segment of length $L$. A legal move chooses a square with side $s\in\{1,\sqrt{2}\}$ and places its base on the segment, so the occupied interval has length $s$. If the left edge is at $x$, the move removes $[x,x+s]$ and leaves two independent subgames of lengths $x$ and $L-s-x$. The C++, Python, and Java implementations analyze this as a Sprague-Grundy game. They then maximize the scaled average proportion of zero-xor moves $$e(L)=\frac{L}{2}\left(\frac{m(L-1)}{L-1}+\frac{m(L-\sqrt{2})}{L-\sqrt{2}}\right),$$ and return $$f(A,B)=\max_{L\in[A,B]} e(L), \qquad (A,B)=(200,500).$$ Here $m(u)$ is the total length of placements that split the remaining segment $u$ into two subpositions with equal Grundy value. The task is therefore a continuous impartial-game calculation followed by a one-dimensional maximization problem. Mathematical Approach The derivation has two layers. First we compute the Grundy function $g(L)$ for a real segment length. Then we build a second function $m(u)$ that counts how often a move lands on xor $0$, and finally maximize the resulting explicit formula for $e(L)$. Step 1: Write the continuous Grundy recurrence For $0\le L<1$ no square fits, so the position is terminal and $$g(L)=0.$$ For larger $L$, a move of side $s\in\{1,\sqrt{2}\}$ can start at any $x$ with $0\le x\le L-s$....

Detailed mathematical approach

Problem Summary

A position is a line segment of length $L$. A legal move chooses a square with side $s\in\{1,\sqrt{2}\}$ and places its base on the segment, so the occupied interval has length $s$. If the left edge is at $x$, the move removes $[x,x+s]$ and leaves two independent subgames of lengths $x$ and $L-s-x$.

The C++, Python, and Java implementations analyze this as a Sprague-Grundy game. They then maximize the scaled average proportion of zero-xor moves

$$e(L)=\frac{L}{2}\left(\frac{m(L-1)}{L-1}+\frac{m(L-\sqrt{2})}{L-\sqrt{2}}\right),$$

and return

$$f(A,B)=\max_{L\in[A,B]} e(L), \qquad (A,B)=(200,500).$$

Here $m(u)$ is the total length of placements that split the remaining segment $u$ into two subpositions with equal Grundy value. The task is therefore a continuous impartial-game calculation followed by a one-dimensional maximization problem.

Mathematical Approach

The derivation has two layers. First we compute the Grundy function $g(L)$ for a real segment length. Then we build a second function $m(u)$ that counts how often a move lands on xor $0$, and finally maximize the resulting explicit formula for $e(L)$.

Step 1: Write the continuous Grundy recurrence

For $0\le L<1$ no square fits, so the position is terminal and

$$g(L)=0.$$

For larger $L$, a move of side $s\in\{1,\sqrt{2}\}$ can start at any $x$ with $0\le x\le L-s$. The two remaining parts are independent, so the nimber of that move is

$$g(x)\oplus g(L-s-x).$$

Therefore

$$g(L)=\operatorname{mex}\Bigl(\{g(x)\oplus g(L-1-x):0\le x\le L-1\}\cup\{g(x)\oplus g(L-\sqrt{2}-x):0\le x\le L-\sqrt{2}\}\Bigr),$$

where the second set is absent when $L<\sqrt{2}$. This is the usual Sprague-Grundy recurrence, except that the game parameter is a real length instead of an integer size.

Step 2: Explain why all breakpoints are of the form $a+b\sqrt{2}$

Every move removes either $1$ or $\sqrt{2}$ from the total occupied length before the remainder is split. Starting from $0$, the only lengths at which the move structure can change are sums built from these two generators:

$$a+b\sqrt{2},\qquad a,b\in\mathbb{Z}_{\ge 0}.$$

Because $\sqrt{2}$ is irrational, different pairs $(a,b)$ produce different real numbers. After sorting these values, any open interval between consecutive breakpoints has the same interaction pattern with all earlier intervals shifted by $1$ and by $\sqrt{2}$. Hence the reachable xor-set stays constant throughout that interval, so $g(L)$ is piecewise constant.

This is the structural simplification used by the implementations: instead of sampling arbitrary real lengths, they store merged intervals on which one Grundy value applies everywhere.

Step 3: Convert zero-xor moves into a measure function

Fix a side length $s$ and a segment length $L$. A placement at $x$ is strategically important when it leaves xor $0$, because the next player then receives a $P$-position. The condition is

$$g(x)\oplus g(L-s-x)=0,$$

which is equivalent to

$$g(x)=g(L-s-x).$$

Let $u=L-s$. The set of such placements has total length

$$m(u)=\operatorname{meas}\left(\left\{x\in[0,u]:g(x)=g(u-x)\right\}\right).$$

Since the legal interval of left endpoints has length $u$, the fraction of zero-xor placements for side $s$ is

$$\frac{m(L-s)}{L-s}.$$

The objective function is exactly the average of those two fractions, scaled by $L$:

$$e(L)=\frac{L}{2}\left(\frac{m(L-1)}{L-1}+\frac{m(L-\sqrt{2})}{L-\sqrt{2}}\right).$$

Step 4: Build $m(u)$ by convolving equal-Grundy intervals

Suppose two intervals $I=[A,B]$ and $J=[C,D]$ carry the same Grundy value. For a fixed $u$, a point $x$ contributes to $m(u)$ when $x\in I$ and $u-x\in J$. Equivalently,

$$x\in I\cap [u-D,u-C].$$

So the contribution from this pair is the overlap length

$$\lambda_{I,J}(u)=\operatorname{meas}\left(I\cap [u-D,u-C]\right).$$

As $u$ varies, $\lambda_{I,J}(u)$ is piecewise linear: it increases linearly, may stay flat, then decreases linearly. Its slope can only change when one endpoint meets another, namely at

$$u=A+C,\quad A+D,\quad B+C,\quad B+D.$$

Summing these trapezoidal contributions over all interval pairs with equal Grundy value yields a global piecewise-linear spline

$$m(u)=s_i u+c_i \qquad (u\in I_i).$$

Distinct intervals matter in both orders $(I,J)$ and $(J,I)$, so they receive double weight; a diagonal pair contributes once.

Step 5: Reduce the maximization to endpoints and derivative roots

Once $m(u)$ is piecewise linear, the function $e(L)$ becomes elementary on every interval where both shifted arguments stay inside fixed spline pieces. If

$$m(L-1)=s_1(L-1)+c_1,\qquad m(L-\sqrt{2})=s_2(L-\sqrt{2})+c_2,$$

then

$$e(L)=\frac{L}{2}\left(s_1+\frac{c_1}{L-1}+s_2+\frac{c_2}{L-\sqrt{2}}\right).$$

A derivative with the same zeros as $e'(L)$ is

$$\left(s_1+s_2\right)-\frac{c_1}{(L-1)^2}-\frac{c_2\sqrt{2}}{(L-\sqrt{2})^2}.$$

Therefore every local maximum must occur at an interval endpoint or at an interior root of this expression. The continuous search on $[A,B]$ is reduced to a finite collection of one-dimensional root checks.

Worked Example: Why $e(2)=2$

For every $u<1$, no square fits on a segment of length $u$, so $g(u)=0$. Consequently, if $0\le u\le 1$, then $g(x)=g(u-x)=0$ for almost every $x\in[0,u]$, and therefore

$$m(u)=u.$$

Now set $L=2$. Then

$$L-1=1,\qquad L-\sqrt{2}=2-\sqrt{2}<1.$$

For side $1$, every placement leaves two subsegments whose total length is $1$, so both subsegments stay below length $1$ except at measure-zero endpoints; the xor is always $0$. For side $\sqrt{2}$, the remaining total length is already below $1$, so again every placement gives xor $0$. Hence

$$m(1)=1,\qquad m(2-\sqrt{2})=2-\sqrt{2}.$$

Substituting into the formula gives

$$e(2)=\frac{2}{2}\left(\frac{1}{1}+\frac{2-\sqrt{2}}{2-\sqrt{2}}\right)=2,$$

which matches the first validation checkpoint used by the implementations.

How the Code Works

The C++ and Java implementations follow the derivation directly. They first enumerate all breakpoints $a+b\sqrt{2}$ up to the needed limit, sort them, and scan the resulting intervals from left to right. On each interval they evaluate the reachable xor-set at a midpoint, take the mex, and merge adjacent intervals whenever the Grundy value remains unchanged.

Next they regroup those intervals by Grundy value and form the piecewise-linear spline for $m(u)$. Each same-Grundy interval pair contributes four event positions where the slope can change, so after sorting the events a single sweep reconstructs every linear piece $s_i u+c_i$.

Finally they build the $L$-breakpoints induced by shifting spline boundaries by $1$ and $\sqrt{2}$. On each such interval, the implementation checks the endpoints, looks for sign changes in the derivative expression, and applies bisection when an interior critical point exists. The Python implementation is only a thin wrapper around the C++ solver, while the C++ version also verifies checkpoints before printing the final value, including $e(2)=2$, $e(4)=1.11974851$, $f(2,10)=2.61969775$, and $f(10,20)=5.99374121$.

Complexity Analysis

Let $B$ be the number of raw breakpoints $a+b\sqrt{2}$ up to the search limit, and let $S$ be the number of merged Grundy intervals after equal neighboring values are coalesced. Generating and sorting the raw breakpoints costs $O(B\log B)$.

Building the piecewise-constant Grundy table is worst-case $O(S^2)$, because each new interval may scan a substantial portion of the previously constructed intervals to collect reachable xor values. Constructing $m(u)$ from equal-Grundy interval pairs is also quadratic in the interval count in the worst case: if $n_k$ intervals carry Grundy value $k$, then the pair generation cost is $O\left(\sum_k n_k^2\right)$, followed by an event sort of the same order up to a logarithmic factor.

The final maximization pass is linear in the number of induced $L$-intervals, with only constant-time evaluations plus a fixed number of bisection iterations per interval. Memory usage is $O(S+E)$, where $E$ is the number of event records in the spline construction.

Footnotes and References

Problem page: Project Euler 644
Sprague-Grundy theorem: Wikipedia - Sprague-Grundy theorem
Minimal excluded value: Wikipedia - mex
Convolution: Wikipedia - Convolution
Bisection method: Wikipedia - Bisection method

Problem 644 source code

C++

#include <algorithm>
#include <cmath>
#include <cstdlib>
#include <iomanip>
#include <iostream>
#include <limits>
#include <thread>
#include <vector>

// The squares game reduces to placing intervals of length 1 or sqrt(2) without overlap.
struct Segment {
    long double l;
    long double r;
    int grundy;
};

struct MSpline {
    long double l;
    long double r;
    long double slope;
    long double c;
};

struct Event {
    double pos;
    double delta;
    bool operator<(const Event& other) const {
        return pos < other.pos;
    }
};

static constexpr long double kEps = 1e-12L;

static std::vector<long double> build_breakpoints(long double limit, long double sqrt2) {
    std::vector<long double> vals;
    int max_b = static_cast<int>(limit / sqrt2) + 2;
    vals.reserve(static_cast<size_t>(limit * limit / 2));
    for (int b = 0; b <= max_b; ++b) {
        long double base = sqrt2 * static_cast<long double>(b);
        if (base > limit + kEps) break;
        int a_max = static_cast<int>(std::floor(limit - base + kEps));
        for (int a = 0; a <= a_max; ++a) {
            vals.push_back(base + static_cast<long double>(a));
        }
    }
    std::sort(vals.begin(), vals.end());
    std::vector<long double> uniq;
    uniq.reserve(vals.size());
    for (long double v : vals) {
        if (uniq.empty() || v - uniq.back() > kEps) {
            uniq.push_back(v);
        }
    }
    if (uniq.empty() || limit - uniq.back() > kEps) {
        uniq.push_back(limit);
    }
    return uniq;
}

static int find_segment(const std::vector<Segment>& segs, long double x) {
    int lo = 0;
    int hi = static_cast<int>(segs.size());
    while (lo < hi) {
        int mid = (lo + hi) >> 1;
        if (segs[mid].r <= x + kEps) {
            lo = mid + 1;
        } else {
            hi = mid;
        }
    }
    if (lo >= static_cast<int>(segs.size())) return static_cast<int>(segs.size()) - 1;
    return lo;
}

static void collect_xors(long double u, const std::vector<Segment>& segs,
                         std::vector<char>& reachable) {
    if (u < -kEps) return;
    if (u <= kEps) {
        reachable[0] = 1;
        return;
    }
    if (segs.empty()) return;
    int i = 0;
    int j = find_segment(segs, u);
    while (i < static_cast<int>(segs.size()) && j >= 0) {
        long double x_l = segs[i].l;
        long double x_r = std::min(segs[i].r, u);
        if (x_l >= u - kEps) break;
        long double y_l = segs[j].l;
        long double y_r = segs[j].r;
        long double left = std::max(x_l, u - y_r);
        long double right = std::min(x_r, u - y_l);
        if (right > left + kEps) {
            int val = segs[i].grundy ^ segs[j].grundy;
            if (val >= static_cast<int>(reachable.size())) {
                reachable.resize(static_cast<size_t>(val + 1), 0);
            }
            reachable[val] = 1;
        }
        if (x_r <= u - y_l + kEps) {
            ++i;
        } else {
            --j;
        }
    }
}

static std::vector<Segment> build_grundy_segments(long double limit, long double sqrt2) {
    std::vector<long double> points = build_breakpoints(limit, sqrt2);
    std::vector<Segment> segs;
    std::vector<char> reachable(128, 0);
    for (size_t i = 0; i + 1 < points.size(); ++i) {
        long double l = points[i];
        long double r = points[i + 1];
        long double mid = (l + r) * 0.5L;
        std::fill(reachable.begin(), reachable.end(), 0);
        if (mid >= 1.0L - kEps) {
            collect_xors(mid - 1.0L, segs, reachable);
        }
        if (mid >= sqrt2 - kEps) {
            collect_xors(mid - sqrt2, segs, reachable);
        }
        int mex = 0;
        while (mex < static_cast<int>(reachable.size()) && reachable[mex]) ++mex;
        if (mex >= static_cast<int>(reachable.size())) {
            reachable.resize(static_cast<size_t>(mex + 1), 0);
        }
        if (!segs.empty() && std::abs(segs.back().r - l) <= kEps && segs.back().grundy == mex) {
            segs.back().r = r;
        } else {
            segs.push_back({l, r, mex});
        }
    }
    return segs;
}

static std::vector<MSpline> build_m_splines(const std::vector<Segment>& gsegs,
                                            long double limit) {
    // m(u) = measure of x where Grundy(x) == Grundy(u-x); build as a piecewise-linear spline.
    int max_g = 0;
    for (const auto& seg : gsegs) max_g = std::max(max_g, seg.grundy);
    std::vector<std::vector<Segment>> per_g(static_cast<size_t>(max_g + 1));
    for (const auto& seg : gsegs) {
        per_g[seg.grundy].push_back(seg);
    }

    std::vector<Event> events;
    for (const auto& lst : per_g) {
        size_t n = lst.size();
        for (size_t i = 0; i < n; ++i) {
            for (size_t j = i; j < n; ++j) {
                long double weight = (i == j) ? 1.0L : 2.0L;
                long double A = lst[i].l;
                long double B = lst[i].r;
                long double C = lst[j].l;
                long double D = lst[j].r;
                long double p1 = A + C;
                long double p2 = A + D;
                long double p3 = B + C;
                long double p4 = B + D;
                if (p1 <= limit + kEps) events.push_back({static_cast<double>(p1), static_cast<double>(weight)});
                if (p2 <= limit + kEps) events.push_back({static_cast<double>(p2), static_cast<double>(-weight)});
                if (p3 <= limit + kEps) events.push_back({static_cast<double>(p3), static_cast<double>(-weight)});
                if (p4 <= limit + kEps) events.push_back({static_cast<double>(p4), static_cast<double>(weight)});
            }
        }
    }

    std::sort(events.begin(), events.end());
    std::vector<MSpline> splines;
    long double prev = 0.0L;
    long double slope = 0.0L;
    long double value = 0.0L;
    size_t idx = 0;
    while (idx < events.size()) {
        long double pos = static_cast<long double>(events[idx].pos);
        if (pos > limit + kEps) break;
        if (pos > prev + kEps) {
            long double c = value - slope * prev;
            splines.push_back({prev, pos, slope, c});
            value += slope * (pos - prev);
            prev = pos;
        }
        long double delta = 0.0L;
        while (idx < events.size() && std::abs(static_cast<long double>(events[idx].pos) - pos) <= kEps) {
            delta += static_cast<long double>(events[idx].delta);
            ++idx;
        }
        slope += delta;
    }
    if (prev < limit - kEps) {
        long double c = value - slope * prev;
        splines.push_back({prev, limit, slope, c});
    }
    if (splines.empty()) {
        splines.push_back({0.0L, limit, 0.0L, 0.0L});
    }
    return splines;
}

static int find_m_segment(const std::vector<MSpline>& splines, long double x) {
    int lo = 0;
    int hi = static_cast<int>(splines.size());
    while (lo < hi) {
        int mid = (lo + hi) >> 1;
        if (splines[mid].r <= x + kEps) {
            lo = mid + 1;
        } else {
            hi = mid;
        }
    }
    if (lo >= static_cast<int>(splines.size())) return static_cast<int>(splines.size()) - 1;
    return lo;
}

static long double eval_m(const std::vector<MSpline>& splines, long double u) {
    if (u <= 0.0L) return 0.0L;
    int idx = find_m_segment(splines, u);
    long double val = splines[idx].slope * u + splines[idx].c;
    if (val < 0.0L && val > -1e-9L) val = 0.0L;
    return val;
}

static long double e_value(const std::vector<MSpline>& splines, long double sqrt2, long double L) {
    long double u1 = L - 1.0L;
    long double u2 = L - sqrt2;
    if (u1 <= 0.0L || u2 <= 0.0L) return 0.0L;
    long double m1 = eval_m(splines, u1);
    long double m2 = eval_m(splines, u2);
    return 0.5L * L * (m1 / u1 + m2 / u2);
}

static long double fprime(long double L, long double sqrt2,
                          long double s1, long double c1,
                          long double s2, long double c2) {
    long double t1 = L - 1.0L;
    long double t2 = L - sqrt2;
    return (s1 + s2) - c1 / (t1 * t1) - c2 * sqrt2 / (t2 * t2);
}

static long double find_root(long double a, long double b, long double sqrt2,
                             long double s1, long double c1,
                             long double s2, long double c2) {
    long double fa = fprime(a, sqrt2, s1, c1, s2, c2);
    long double fb = fprime(b, sqrt2, s1, c1, s2, c2);
    if (std::abs(fa) <= 1e-18L) return a;
    if (std::abs(fb) <= 1e-18L) return b;
    if (fa * fb > 0.0L) return std::numeric_limits<long double>::quiet_NaN();
    for (int it = 0; it < 80; ++it) {
        long double mid = (a + b) * 0.5L;
        long double fm = fprime(mid, sqrt2, s1, c1, s2, c2);
        if (fa * fm <= 0.0L) {
            b = mid;
            fb = fm;
        } else {
            a = mid;
            fa = fm;
        }
    }
    return (a + b) * 0.5L;
}

static long double max_e_in_range(const std::vector<MSpline>& splines,
                                  const std::vector<long double>& breakpoints,
                                  long double sqrt2,
                                  long double Lmin, long double Lmax,
                                  unsigned threads) {
    int start = static_cast<int>(std::lower_bound(breakpoints.begin(), breakpoints.end(), Lmin - kEps) - breakpoints.begin());
    int end = static_cast<int>(std::upper_bound(breakpoints.begin(), breakpoints.end(), Lmax + kEps) - breakpoints.begin());
    int last_interval = static_cast<int>(breakpoints.size()) - 1;
    if (start >= last_interval) return 0.0L;
    if (end > last_interval) end = last_interval;
    if (end <= start) return 0.0L;

    int total = end - start;
    if (threads == 0) threads = 1;
    if (threads > static_cast<unsigned>(total)) threads = static_cast<unsigned>(total);

    std::vector<long double> best_per_thread(threads, 0.0L);
    std::vector<std::thread> pool;
    pool.reserve(threads);
    for (unsigned t = 0; t < threads; ++t) {
        int chunk_start = start + (total * static_cast<int>(t)) / static_cast<int>(threads);
        int chunk_end = start + (total * static_cast<int>(t + 1)) / static_cast<int>(threads);
        pool.emplace_back([&, t, chunk_start, chunk_end]() {
            long double local_best = 0.0L;
            for (int i = chunk_start; i < chunk_end; ++i) {
                long double a = std::max(breakpoints[i], Lmin);
                long double b = std::min(breakpoints[i + 1], Lmax);
                if (b <= a + kEps) continue;

                long double mid = (a + b) * 0.5L;
                long double u1 = mid - 1.0L;
                long double u2 = mid - sqrt2;
                if (u1 <= 0.0L || u2 <= 0.0L) continue;
                int idx1 = find_m_segment(splines, u1);
                int idx2 = find_m_segment(splines, u2);
                long double s1 = splines[idx1].slope;
                long double c1 = splines[idx1].c;
                long double s2 = splines[idx2].slope;
                long double c2 = splines[idx2].c;

                long double ea = e_value(splines, sqrt2, a);
                long double eb = e_value(splines, sqrt2, b);
                local_best = std::max(local_best, std::max(ea, eb));

                long double points[5] = {a, a + (b - a) * 0.25L, mid, a + (b - a) * 0.75L, b};
                long double fp[5];
                for (int p = 0; p < 5; ++p) {
                    fp[p] = fprime(points[p], sqrt2, s1, c1, s2, c2);
                }
                for (int p = 0; p < 4; ++p) {
                    if (fp[p] == 0.0L) {
                        long double ev = e_value(splines, sqrt2, points[p]);
                        local_best = std::max(local_best, ev);
                        continue;
                    }
                    if (fp[p] * fp[p + 1] <= 0.0L) {
                        long double root = find_root(points[p], points[p + 1], sqrt2, s1, c1, s2, c2);
                        if (root == root) {
                            long double ev = e_value(splines, sqrt2, root);
                            local_best = std::max(local_best, ev);
                        }
                    }
                }
            }
            best_per_thread[t] = local_best;
        });
    }
    for (auto& th : pool) th.join();
    long double best = 0.0L;
    for (long double v : best_per_thread) best = std::max(best, v);
    return best;
}

int main(int argc, char** argv) {
    long double Lmin = 200.0L;
    long double Lmax = 500.0L;
    if (argc > 2) {
        Lmin = std::strtold(argv[1], nullptr);
        Lmax = std::strtold(argv[2], nullptr);
    }
    long double sqrt2 = std::sqrt(2.0L);

    long double grundy_limit = Lmax;
    auto gsegs = build_grundy_segments(grundy_limit, sqrt2);

    long double u_limit = Lmax - 1.0L;
    auto splines = build_m_splines(gsegs, u_limit);

    long double global_min = std::min(Lmin, 2.0L);
    std::vector<long double> breakpoints;
    breakpoints.reserve(splines.size() * 4 + 8);
    breakpoints.push_back(global_min);
    breakpoints.push_back(Lmax);
    for (const auto& sp : splines) {
        long double a = sp.l + 1.0L;
        long double b = sp.r + 1.0L;
        long double c = sp.l + sqrt2;
        long double d = sp.r + sqrt2;
        if (a <= Lmax + kEps) breakpoints.push_back(a);
        if (b <= Lmax + kEps) breakpoints.push_back(b);
        if (c <= Lmax + kEps) breakpoints.push_back(c);
        if (d <= Lmax + kEps) breakpoints.push_back(d);
    }
    std::sort(breakpoints.begin(), breakpoints.end());
    std::vector<long double> uniq;
    uniq.reserve(breakpoints.size());
    for (long double v : breakpoints) {
        if (v < global_min - kEps || v > Lmax + kEps) continue;
        if (uniq.empty() || v - uniq.back() > kEps) uniq.push_back(v);
    }
    if (uniq.empty() || global_min - uniq.front() < -kEps) uniq.insert(uniq.begin(), global_min);
    if (uniq.empty() || Lmax - uniq.back() > kEps) uniq.push_back(Lmax);

    unsigned threads = std::max(1u, std::thread::hardware_concurrency());
    long double ans = max_e_in_range(splines, uniq, sqrt2, Lmin, Lmax, threads);

    // Validation checks.
    long double e2 = e_value(splines, sqrt2, 2.0L);
    if (std::abs(e2 - 2.0L) > 1e-7L) {
        std::cerr << "Validation failed: e(2) got " << std::fixed << std::setprecision(12) << e2
                  << " expected 2\n";
        return 1;
    }
    long double e4 = e_value(splines, sqrt2, 4.0L);
    if (std::abs(e4 - 1.11974851L) > 2e-6L) {
        std::cerr << "Validation failed: e(4) got " << std::fixed << std::setprecision(12) << e4
                  << " expected 1.11974851\n";
        return 1;
    }
    long double f2_10 = max_e_in_range(splines, uniq, sqrt2, 2.0L, 10.0L, threads);
    if (std::abs(f2_10 - 2.61969775L) > 2e-6L) {
        std::cerr << "Validation failed: f(2,10) got " << std::fixed << std::setprecision(12) << f2_10
                  << " expected 2.61969775\n";
        return 1;
    }
    long double f10_20 = max_e_in_range(splines, uniq, sqrt2, 10.0L, 20.0L, threads);
    if (std::abs(f10_20 - 5.99374121L) > 2e-6L) {
        std::cerr << "Validation failed: f(10,20) got " << std::fixed << std::setprecision(12) << f10_20
                  << " expected 5.99374121\n";
        return 1;
    }

    std::cout.setf(std::ios::fixed);
    std::cout.precision(8);
    std::cout << static_cast<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.ArrayList;
import java.util.Collections;
import java.util.Locale;

public class Euler644 {

    static final double kEps = 1e-12;

    static class Segment {
        double l, r;
        int grundy;

        Segment(double l, double r, int grundy) {
            this.l = l;
            this.r = r;
            this.grundy = grundy;
        }
    }

    static class MSpline {
        double l, r, slope, c;

        MSpline(double l, double r, double slope, double c) {
            this.l = l;
            this.r = r;
            this.slope = slope;
            this.c = c;
        }
    }

    static class Event implements Comparable<Event> {
        double pos, delta;

        Event(double pos, double delta) {
            this.pos = pos;
            this.delta = delta;
        }

        @Override
        public int compareTo(Event o) {
            return Double.compare(this.pos, o.pos);
        }
    }

    static ArrayList<Double> buildBreakpoints(double limit, double sqrt2) {
        ArrayList<Double> vals = new ArrayList<>();
        int max_b = (int) (limit / sqrt2) + 2;
        for (int b = 0; b <= max_b; ++b) {
            double base = sqrt2 * b;
            if (base > limit + kEps)
                break;
            int a_max = (int) Math.floor(limit - base + kEps);
            for (int a = 0; a <= a_max; ++a) {
                vals.add(base + a);
            }
        }
        Collections.sort(vals);
        ArrayList<Double> uniq = new ArrayList<>();
        for (double v : vals) {
            if (uniq.isEmpty() || v - uniq.get(uniq.size() - 1) > kEps) {
                uniq.add(v);
            }
        }
        if (uniq.isEmpty() || limit - uniq.get(uniq.size() - 1) > kEps) {
            uniq.add(limit);
        }
        return uniq;
    }

    static int findSegment(ArrayList<Segment> segs, double x) {
        int lo = 0;
        int hi = segs.size();
        while (lo < hi) {
            int mid = (lo + hi) >> 1;
            if (segs.get(mid).r <= x + kEps) {
                lo = mid + 1;
            } else {
                hi = mid;
            }
        }
        if (lo >= segs.size())
            return segs.size() - 1;
        return lo;
    }

    static void collectXors(double u, ArrayList<Segment> segs, ArrayList<Integer> reachable) {
        if (u < -kEps)
            return;
        if (u <= kEps) {
            reachable.set(0, 1);
            return;
        }
        if (segs.isEmpty())
            return;
        int i = 0;
        int j = findSegment(segs, u);
        while (i < segs.size() && j >= 0) {
            double x_l = segs.get(i).l;
            double x_r = Math.min(segs.get(i).r, u);
            if (x_l >= u - kEps)
                break;
            double y_l = segs.get(j).l;
            double y_r = segs.get(j).r;
            double left = Math.max(x_l, u - y_r);
            double right = Math.min(x_r, u - y_l);
            if (right > left + kEps) {
                int val = segs.get(i).grundy ^ segs.get(j).grundy;
                while (val >= reachable.size())
                    reachable.add(0);
                reachable.set(val, 1);
            }
            if (x_r <= u - y_l + kEps) {
                ++i;
            } else {
                --j;
            }
        }
    }

    static ArrayList<Segment> buildGrundySegments(double limit, double sqrt2) {
        ArrayList<Double> points = buildBreakpoints(limit, sqrt2);
        ArrayList<Segment> segs = new ArrayList<>();
        ArrayList<Integer> reachable = new ArrayList<>(128);
        for (int i = 0; i < 128; i++)
            reachable.add(0);
        for (int i = 0; i + 1 < points.size(); ++i) {
            double l = points.get(i);
            double r = points.get(i + 1);
            double mid = (l + r) * 0.5;
            for (int k = 0; k < reachable.size(); k++)
                reachable.set(k, 0);
            if (mid >= 1.0 - kEps)
                collectXors(mid - 1.0, segs, reachable);
            if (mid >= sqrt2 - kEps)
                collectXors(mid - sqrt2, segs, reachable);
            int mex = 0;
            while (mex < reachable.size() && reachable.get(mex) == 1)
                ++mex;
            while (mex >= reachable.size())
                reachable.add(0);
            if (!segs.isEmpty() && Math.abs(segs.get(segs.size() - 1).r - l) <= kEps
                    && segs.get(segs.size() - 1).grundy == mex) {
                segs.get(segs.size() - 1).r = r;
            } else {
                segs.add(new Segment(l, r, mex));
            }
        }
        return segs;
    }

    static ArrayList<MSpline> buildMSplines(ArrayList<Segment> gsegs, double limit) {
        int max_g = 0;
        for (Segment seg : gsegs)
            max_g = Math.max(max_g, seg.grundy);
        ArrayList<ArrayList<Segment>> per_g = new ArrayList<>(max_g + 1);
        for (int i = 0; i <= max_g; i++)
            per_g.add(new ArrayList<>());
        for (Segment seg : gsegs)
            per_g.get(seg.grundy).add(seg);

        ArrayList<Event> events = new ArrayList<>();
        for (ArrayList<Segment> lst : per_g) {
            int n = lst.size();
            for (int i = 0; i < n; ++i) {
                for (int j = i; j < n; ++j) {
                    double weight = (i == j) ? 1.0 : 2.0;
                    double A = lst.get(i).l;
                    double B = lst.get(i).r;
                    double C = lst.get(j).l;
                    double D = lst.get(j).r;
                    double p1 = A + C;
                    double p2 = A + D;
                    double p3 = B + C;
                    double p4 = B + D;
                    if (p1 <= limit + kEps)
                        events.add(new Event(p1, weight));
                    if (p2 <= limit + kEps)
                        events.add(new Event(p2, -weight));
                    if (p3 <= limit + kEps)
                        events.add(new Event(p3, -weight));
                    if (p4 <= limit + kEps)
                        events.add(new Event(p4, weight));
                }
            }
        }

        Collections.sort(events);
        ArrayList<MSpline> splines = new ArrayList<>();
        double prev = 0.0;
        double slope = 0.0;
        double value = 0.0;
        int idx = 0;
        while (idx < events.size()) {
            double pos = events.get(idx).pos;
            if (pos > limit + kEps)
                break;
            if (pos > prev + kEps) {
                double c = value - slope * prev;
                splines.add(new MSpline(prev, pos, slope, c));
                value += slope * (pos - prev);
                prev = pos;
            }
            double delta = 0.0;
            while (idx < events.size() && Math.abs(events.get(idx).pos - pos) <= kEps) {
                delta += events.get(idx).delta;
                ++idx;
            }
            slope += delta;
        }
        if (prev < limit - kEps) {
            double c = value - slope * prev;
            splines.add(new MSpline(prev, limit, slope, c));
        }
        if (splines.isEmpty()) {
            splines.add(new MSpline(0.0, limit, 0.0, 0.0));
        }
        return splines;
    }

    static int findMSegment(ArrayList<MSpline> splines, double x) {
        int lo = 0;
        int hi = splines.size();
        while (lo < hi) {
            int mid = (lo + hi) >> 1;
            if (splines.get(mid).r <= x + kEps) {
                lo = mid + 1;
            } else {
                hi = mid;
            }
        }
        if (lo >= splines.size())
            return splines.size() - 1;
        return lo;
    }

    static double evalM(ArrayList<MSpline> splines, double u) {
        if (u <= 0.0)
            return 0.0;
        int idx = findMSegment(splines, u);
        double val = splines.get(idx).slope * u + splines.get(idx).c;
        if (val < 0.0 && val > -1e-9)
            val = 0.0;
        return val;
    }

    static double eValue(ArrayList<MSpline> splines, double sqrt2, double L) {
        double u1 = L - 1.0;
        double u2 = L - sqrt2;
        if (u1 <= 0.0 || u2 <= 0.0)
            return 0.0;
        double m1 = evalM(splines, u1);
        double m2 = evalM(splines, u2);
        return 0.5 * L * (m1 / u1 + m2 / u2);
    }

    static double fPrime(double L, double sqrt2, double s1, double c1, double s2, double c2) {
        double t1 = L - 1.0;
        double t2 = L - sqrt2;
        return (s1 + s2) - c1 / (t1 * t1) - c2 * sqrt2 / (t2 * t2);
    }

    static double findRoot(double a, double b, double sqrt2, double s1, double c1, double s2, double c2) {
        double fa = fPrime(a, sqrt2, s1, c1, s2, c2);
        double fb = fPrime(b, sqrt2, s1, c1, s2, c2);
        if (Math.abs(fa) <= 1e-18)
            return a;
        if (Math.abs(fb) <= 1e-18)
            return b;
        if (fa * fb > 0.0)
            return Double.NaN;
        for (int it = 0; it < 80; ++it) {
            double mid = (a + b) * 0.5;
            double fm = fPrime(mid, sqrt2, s1, c1, s2, c2);
            if (fa * fm <= 0.0) {
                b = mid;
                fb = fm;
            } else {
                a = mid;
                fa = fm;
            }
        }
        return (a + b) * 0.5;
    }

    static double maxEInRange(ArrayList<MSpline> splines, ArrayList<Double> breakpoints, double sqrt2, double Lmin,
            double Lmax) {
        int start = 0;
        while (start < breakpoints.size() && breakpoints.get(start) < Lmin - kEps)
            start++;
        int end = 0;
        while (end < breakpoints.size() && breakpoints.get(end) <= Lmax + kEps)
            end++;
        int last_interval = breakpoints.size() - 1;
        if (start >= last_interval)
            return 0.0;
        if (end > last_interval)
            end = last_interval;
        if (end <= start)
            return 0.0;

        double best = 0.0;
        for (int i = start; i < end; ++i) {
            double a = Math.max(breakpoints.get(i), Lmin);
            double b = Math.min(breakpoints.get(i + 1), Lmax);
            if (b <= a + kEps)
                continue;

            double mid = (a + b) * 0.5;
            double u1 = mid - 1.0;
            double u2 = mid - sqrt2;
            if (u1 <= 0.0 || u2 <= 0.0)
                continue;
            int idx1 = findMSegment(splines, u1);
            int idx2 = findMSegment(splines, u2);
            double s1 = splines.get(idx1).slope;
            double c1 = splines.get(idx1).c;
            double s2 = splines.get(idx2).slope;
            double c2 = splines.get(idx2).c;

            double ea = eValue(splines, sqrt2, a);
            double eb = eValue(splines, sqrt2, b);
            best = Math.max(best, Math.max(ea, eb));

            double[] points = { a, a + (b - a) * 0.25, mid, a + (b - a) * 0.75, b };
            double[] fp = new double[5];
            for (int p = 0; p < 5; ++p)
                fp[p] = fPrime(points[p], sqrt2, s1, c1, s2, c2);
            for (int p = 0; p < 4; ++p) {
                if (fp[p] == 0.0) {
                    double ev = eValue(splines, sqrt2, points[p]);
                    best = Math.max(best, ev);
                    continue;
                }
                if (fp[p] * fp[p + 1] <= 0.0) {
                    double root = findRoot(points[p], points[p + 1], sqrt2, s1, c1, s2, c2);
                    if (!Double.isNaN(root)) {
                        double ev = eValue(splines, sqrt2, root);
                        best = Math.max(best, ev);
                    }
                }
            }
        }
        return best;
    }

    public static String solve() {
        double Lmin = 200.0;
        double Lmax = 500.0;
        double sqrt2 = Math.sqrt(2.0);

        ArrayList<Segment> gsegs = buildGrundySegments(Lmax, sqrt2);
        ArrayList<MSpline> splines = buildMSplines(gsegs, Lmax - 1.0);

        double global_min = Math.min(Lmin, 2.0);
        ArrayList<Double> breakpoints = new ArrayList<>();
        breakpoints.add(global_min);
        breakpoints.add(Lmax);
        for (MSpline sp : splines) {
            double a = sp.l + 1.0;
            double b = sp.r + 1.0;
            double c = sp.l + sqrt2;
            double d = sp.r + sqrt2;
            if (a <= Lmax + kEps)
                breakpoints.add(a);
            if (b <= Lmax + kEps)
                breakpoints.add(b);
            if (c <= Lmax + kEps)
                breakpoints.add(c);
            if (d <= Lmax + kEps)
                breakpoints.add(d);
        }
        Collections.sort(breakpoints);
        ArrayList<Double> uniq = new ArrayList<>();
        for (double v : breakpoints) {
            if (v < global_min - kEps || v > Lmax + kEps)
                continue;
            if (uniq.isEmpty() || v - uniq.get(uniq.size() - 1) > kEps)
                uniq.add(v);
        }
        if (uniq.isEmpty() || global_min - uniq.get(0) < -kEps)
            uniq.add(0, global_min);
        if (uniq.isEmpty() || Lmax - uniq.get(uniq.size() - 1) > kEps)
            uniq.add(Lmax);

        double ans = maxEInRange(splines, uniq, sqrt2, Lmin, Lmax);
        return String.format(Locale.US, "%.8f", ans);
    }

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