Problem 287: Quadtree Encoding (a Simple Compression Algorithm)

View on Project Euler

Project Euler Problem 287 Solution

EulerSolve provides an optimized solution for Project Euler Problem 287, Quadtree Encoding (a Simple Compression Algorithm), with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary The image has size $$2^N\times 2^N,$$ with pixel coordinates \(0\le x,y\le 2^N-1\). A pixel is black exactly when $$ (x-c)^2+(y-c)^2\le R^2,\qquad c=2^{N-1},\quad R^2=2^{2N-2}. $$ So the image is a digital disk centered at \((c,c)\). The quadtree encoding uses: $$\text{uniform block} \Rightarrow 2\text{ bits},\qquad \text{mixed block} \Rightarrow 1+\sum_{i=1}^4 L_i.$$ The goal is to compute the total encoding length for \(N=24\) without expanding all pixels. Mathematical Approach 1) Length depends only on whether a block is uniform. If every pixel in a square block has the same color, the encoder stops immediately and emits a leaf of length \(2\). Otherwise it emits one split bit and recurses into the four children. Since black and white leaves both cost \(2\), we never need to know which color a uniform block has, only whether it is uniform. 2) Reformulate the color test geometrically. For a block $$B=[x_0,x_1]\times [y_0,y_1]$$ with integer pixel coordinates, define the squared distance to the center $$D(x,y)=(x-c)^2+(y-c)^2.$$ The block is: $$\text{all black if } \max_{(x,y)\in B} D(x,y)\le R^2,$$ $$\text{all white if } \min_{(x,y)\in B} D(x,y)>R^2.$$ Only if neither condition holds is the block mixed. 3) Exact formula for the nearest pixel....

Detailed mathematical approach

Problem Summary

For a positive integer \(N\), a black and white image of \(2^N\times2^N\) pixels is described by a sequence of bits, read from left to right. The first bit concerns the whole image. The bit \(0\) splits the current square region of side \(2^n\) into four squares of side \(2^{n-1}\), whose descriptions follow in the order top left, top right, bottom left, bottom right; the pair \(10\) says that the current region is entirely black, and \(11\) that it is entirely white. The \(4\times4\) example of the statement has a description of length \(30\) and a minimal one of length \(16\). The image \(D_N\) consists of the pixels \((x,y)\) with \(0\le x,y\le2^N-1\), where \((0,0)\) is the bottom left pixel, and a pixel is black if

$$(x-2^{N-1})^2+(y-2^{N-1})^2\le2^{2N-2}$$

and white otherwise. We must find the length of the minimal sequence that describes \(D_{24}\). [1]

The whole solution in one sentence: a shortest description splits exactly the squares that contain both colours, and these are recognised by their pixels nearest to and farthest from the centre of the disk, so a recursion counts the splits while subdividing only squares along the circle.

Mathematical Approach

We first find the shortest description of any image, then exploit the shape of the disk.

1. Reading a description

The minimal sequence of the statement, \(0100101111101110\), reads as follows: \(0\), the image is split; \(10\), its top left quadrant is black; \(0\), the top right quadrant is split into four pixels, which \(10\), \(11\), \(11\), \(10\) make black, white, white and black; \(11\) and \(10\), the bottom left quadrant is white and the bottom right one black. The regions form a tree, the region quadtree of the image: [2] a split region has its four quadrants as children, a region described by \(10\) or \(11\) is a leaf, and the bits list the nodes in pre-order, every node before its children. [3] Pre-order alone does not determine a tree in general, but here the first bit of a node says whether it has four children or none. A single pixel cannot be split.

2. The length counts the splits

Claim. A description with \(z\) splits has length \(7z+2\).

Proof. Every node except the whole image is one of the four children of exactly one split node, so the tree has \(4z+1\) nodes: \(z\) split nodes and \(3z+1\) leaves. A split costs one bit and a leaf two bits, together \(z+2(3z+1)=7z+2\).

The minimal sequence of the statement splits twice, and \(7\cdot2+2=16\). The longer one also splits the top left and the bottom right quadrant, although both are entirely black, and \(7\cdot4+2=30\). A shortest description is therefore one with the fewest splits.

3. The fewest splits

The regions that a description can reach are the blocks: for \(0\le k\le N\), the blocks of level \(k\) are the squares of side \(2^k\) whose bottom left pixel has both coordinates divisible by \(2^k\). A block is uniform if all its pixels have the same colour and mixed otherwise; single pixels are uniform. Let \(M_k\) be the number of mixed blocks of level \(k\).

Claim. Every description splits every mixed block, and splitting exactly the mixed blocks gives a description. Hence the minimal sequence is unique, and its length is

$$L_{\min}=7M+2,\qquad M=M_1+M_2+\dots+M_N.$$

Proof. Take a description and a mixed block \(B\), and follow the description from the whole image into the quadrants that contain \(B\). If this walk stopped at a leaf larger than \(B\), that leaf would be uniform, and so would \(B\). Hence the walk reaches \(B\), and \(B\) is split, because a mixed region cannot be a leaf. Conversely, split a region exactly when it is mixed. A block containing a mixed block is mixed itself, so this process reaches and splits every mixed block and nothing else, and it stops at uniform blocks, at the latest at single pixels. It makes \(M\) splits; every other description splits these \(M\) blocks and at least one more, so by Section 2 it is longer.

The programs compute this length recursively: a uniform block costs \(2\) bits, a mixed block \(1\) bit plus the costs of its four quadrants. Since a reflection of the image maps blocks to blocks, the orientation conventions do not affect the length.

4. The disk

Write \(c=2^{N-1}\) for both coordinates of the centre and \(R=2^{N-1}\) for the radius, so that \(R^2=2^{2N-2}\); a pixel is black exactly when

$$D(x,y)=(x-c)^2+(y-c)^2\le R^2.$$

The middle of the image lies between pixels, at \(c-\tfrac12\) in both coordinates, so the disk is shifted half a pixel up and to the right. It meets the left column and the bottom row in single pixels, \((0,c)\) and \((c,0)\), but the right column and the top row in several, five each for \(N=3\).

5. Worked example: \(D_2\)

For \(N=2\) we have \(c=2\) and \(R^2=4\). Row by row from the top row \(y=3\) down to \(y=0\), with \(\blacksquare\) for black and \(\square\) for white, the image is

$$\begin{matrix}\square&\blacksquare&\blacksquare&\blacksquare\\\blacksquare&\blacksquare&\blacksquare&\blacksquare\\\square&\blacksquare&\blacksquare&\blacksquare\\\square&\square&\blacksquare&\square\end{matrix}$$

The top right quadrant is black, and the other three quadrants contain both colours, so \(M_2=1\), \(M_1=3\) and \(M=4\). The minimal sequence is \(0\), followed by \(011101010\) for the top left quadrant, \(10\) for the top right one, \(011101111\) for the bottom left one and \(010101011\) for the bottom right one. Its length is \(1+9+2+9+9=30=7\cdot4+2\).

6. Two pixels decide a block

Let \(B\) be the block of side \(s\) with the columns \(x_0\le x\le x_1\) and the rows \(y_0\le y\le y_1\), where \(x_1=x_0+s-1\) and \(y_1=y_0+s-1\). Put \(x^*=x_0\) if \(c\lt x_0\), \(x^*=x_1\) if \(c\gt x_1\) and \(x^*=c\) otherwise, that is, \(c\) clamped to the interval from \(x_0\) to \(x_1\), [4] and let \(\hat x\) be the one of \(x_0\) and \(x_1\) farther from \(c\). Define \(y^*\) and \(\hat y\) in the same way.

Claim. Over the pixels of \(B\), the smallest value of \(D\) is \(D_{\min}=D(x^*,y^*)\) and the largest is \(D_{\max}=D(\hat x,\hat y)\). Consequently \(B\) is black if \(D_{\max}\le R^2\), white if \(D_{\min}\gt R^2\), and mixed if \(D_{\min}\le R^2\lt D_{\max}\).

Proof. The function \(t\mapsto(t-c)^2\) decreases for \(t\le c\) and increases for \(t\ge c\). On the integers from \(x_0\) to \(x_1\) it is therefore smallest at the integer closest to \(c\), which is \(x^*\) because \(c\) is an integer, and largest at the end \(\hat x\). Now \(D\) is the sum of a function of \(x\) alone and a function of \(y\) alone, and \(B\) contains every combination of its columns and rows, so the minimum of \(D\) over \(B\) is attained at \((x^*,y^*)\) and the maximum at \((\hat x,\hat y)\). If \(D_{\max}\le R^2\), every pixel is black; if \(D_{\min}\gt R^2\), every pixel is white; otherwise \((x^*,y^*)\) is black and \((\hat x,\hat y)\) is white.

In \(D_2\), the top right quadrant has \(D=0\) at its nearest pixel \((2,2)\) and \(D=2\le4\) at its farthest pixel \((3,3)\), so it is black. The other quadrants have \(D=1\) at \((1,2)\) and \(D=5\) at \((0,3)\), \(D=2\) at \((1,1)\) and \(D=8\) at \((0,0)\), and \(D=1\) at \((2,1)\) and \(D=5\) at \((3,0)\), so they are mixed. The test costs the same few multiplications for a block of any size.

7. A larger example: \(D_3\)

For \(N=3\), with \(c=4\) and \(R^2=16\), the image is, again from the top row down,

$$\begin{matrix}\square&\square&\blacksquare&\blacksquare&\blacksquare&\blacksquare&\blacksquare&\square\\\square&\blacksquare&\blacksquare&\blacksquare&\blacksquare&\blacksquare&\blacksquare&\blacksquare\\\square&\blacksquare&\blacksquare&\blacksquare&\blacksquare&\blacksquare&\blacksquare&\blacksquare\\\blacksquare&\blacksquare&\blacksquare&\blacksquare&\blacksquare&\blacksquare&\blacksquare&\blacksquare\\\square&\blacksquare&\blacksquare&\blacksquare&\blacksquare&\blacksquare&\blacksquare&\blacksquare\\\square&\blacksquare&\blacksquare&\blacksquare&\blacksquare&\blacksquare&\blacksquare&\blacksquare\\\square&\square&\blacksquare&\blacksquare&\blacksquare&\blacksquare&\blacksquare&\square\\\square&\square&\square&\square&\blacksquare&\square&\square&\square\end{matrix}$$

The whole image and its four quadrants are mixed. Of the sixteen blocks of side \(2\), seven are mixed: those with the bottom left pixels \((2,0)\), \((4,0)\), \((6,0)\), \((0,2)\), \((0,4)\), \((0,6)\) and \((6,6)\). Hence \(M=1+4+7=12\), and the minimal length is \(7\cdot12+2=86\). The shift is visible: the bottom left block of side \(2\) is white, the top right one is mixed. The quadrants differ: in the order of the statement their minimal lengths are \(23\), \(16\), \(23\) and \(23\), and \(4685\), \(4503\), \(4916\) and \(4685\) for \(N=10\). Only the symmetry exchanging \(x\) and \(y\) survives, which swaps the top left and the bottom right quadrant; multiplying one quadrant by four gives a wrong length.

8. Only the blocks along the circle are mixed

For \(N=10\) the numbers of mixed blocks of the levels \(k=1,\dots,10\) are \(1023\), \(762\), \(440\), \(232\), \(122\), \(60\), \(28\), \(12\), \(4\) and \(1\), in total \(M=2684\), and the minimal length of \(D_{10}\) is \(18790\). The counts roughly halve from level to level, because a mixed block lies close to the circle.

Claim. For every level \(k\), \(M_k\lt2\sqrt2\,\pi\cdot2^{N-k}\); hence \(M\lt2\sqrt2\,\pi\,(2^N-1)\), less than \(8.9\cdot2^N\).

Proof. A mixed block of side \(s\) contains a black pixel \(p\) and a white pixel \(q\) at distance \(1\): a walk inside the block from a black pixel to a white one, one column or one row at a time, changes colour at some step. With \(C=(c,c)\) and \(|\cdot|\) the Euclidean length we have \(|p-C|\le R\lt|q-C|\le|p-C|+1\), so \(|p-C|\gt R-1\). Draw every pixel as the unit square centred at it: the block becomes a square of side \(s\) whose points all lie within the distance \(a=(s-\tfrac12)\sqrt2\) of \(p\), hence in the ring between the circles of radii \(R-1-a\) and \(R+a\) around \(C\). For \(k\le N-2\) we have \(R\ge2s\), so the inner radius is \(R-1-a\ge(2-\sqrt2)(s-\tfrac12)\gt0\), and the ring has the area

$$\pi\bigl((R+a)^2-(R-1-a)^2\bigr)=\pi(2a+1)(2R-1).$$

The squares of the blocks of one level do not overlap, so \(M_k s^2\le\pi(2a+1)(2R-1)\), and the claim follows from \((2a+1)/s\lt2\sqrt2\) and \((2R-1)/s\lt2^{N-k}\). The levels \(N-1\) and \(N\) have only \(4\) blocks and \(1\) block, which satisfy the bound anyway.

The same comparison of unit squares with an area shows that a disk of radius \(r\) contains \(\pi r^2\) lattice points up to an error of the order of \(r\). [5] Our bound is generous, about \(9090\) for \(N=10\), but it grows like \(2^N\), the side of the image, not like \(4^N\), the number of pixels, about \(2.8\cdot10^{14}\) for \(N=24\).

How the Code Works

The C++ program [6] starts in main. parse_arguments accepts --skip-checkpoints and --n= followed by digits, the last one counting, with the default \(N=24\); any other argument, or a final \(N\) outside \(1\le N\le30\), ends the program with exit code \(1\). The digit loop of parse_int_after_prefix has no overflow check, so values from \(2^{31}\) on overflow its int first; a test build read \(2^{32}+3\) as \(3\). Unless --skip-checkpoints is given, run_checkpoints compares encode_length_fast(n) with encode_length_bruteforce(n) for \(n=1,\dots,5\) and ends the program with exit code \(2\) on a difference. Then main prints encode_length_fast of the chosen \(N\).

encode_length_fast sets size to \(2^N\), c to \(2^{N-1}\) and the unsigned 64-bit r2 to \(R^2\) and calls the recursive lambda rec, which receives itself as self. A call rec(self, x0, y0, side) is Section 6 for one block: nearest_x and nearest_y are \(x^*\) and \(y^*\), computed with std::clamp, and min_d2 is \(D_{\min}\); dx_max is the larger of std::llabs(x0 - c) and std::llabs(x1 - c), likewise dy_max, and max_d2 is \(D_{\max}\). If max_d2 <= r2 or min_d2 > r2, the block is uniform and the call returns \(2\). Otherwise it returns \(1\) plus the lengths of the four quadrants, written in the order of the statement, top left at (x0, y0 + half) since \(y\) grows upwards; the order in which C++ evaluates the four calls is unspecified and does not matter for the sum. By Section 3 the result is \(7M+2\). The coordinates are int, the differences are squared as long long, and the squared distances, at most \(2^{2N-1}\), fit in 64 bits.

encode_length_bruteforce stores the image in a vector of vectors black and builds the summed-area table ps, in which ps[y][x] counts the black pixels with column less than \(x\) and row less than \(y\); [7] sum_black obtains the number of black pixels of a block from four entries of the table: the count below and to the left of its top right corner, minus the counts below the block and to its left, plus the count below and to the left of the block, which was subtracted twice. For the top right quadrant of \(D_2\) this gives \(11-4-4+1=4\). Its own rec calls a block uniform when this number is \(0\) or \(s^2\), so it checks the test of Section 6 against the actual pixels. The bitmap and the table grow like \(4^N\), and run_checkpoints uses them only for \(N\le5\).

The Java program has the same recursion in the static method rec(x0, y0, side, c, r2), with all values of type long and the clamp written as Math.max(x0, Math.min(c, x1)); solve fixes \(N=24\), without options or checkpoints. The Python file is only a bridge to the C++ program: solve compiles solutionsCpp/Euler287.cpp with clang++ or g++ into the hidden file .euler287_py_bridge when that is missing or older than the source, runs it without arguments, checkpoints included, and returns the last non-empty output line via parse_output. A failed checkpoint makes subprocess.check_output raise CalledProcessError. [8]

Complexity and Verification

The recursion makes \(4M+1\) calls, one for the whole image and four for every mixed block, each with a constant number of operations, so by Section 8 the time is \(O(2^N)\). The memory is the recursion stack of at most \(N+1\) nested calls. The brute force needs \(O(4^N)\) time and memory. On an Apple M4 Max the C++ program, built with clang++ -O2, takes about \(0.4\) seconds for \(N=24\) including the checkpoints, and the Java program about \(0.3\) seconds.

For this article we also ran these checks, outside the published solutions:

  • An independent Python brute force builds the minimal sequence from the pixels. For \(N=1,\dots,12\) its lengths \(9\), \(30\), \(86\), \(212\), \(499\), \(1052\), \(2242\), \(4552\), \(9382\), \(18790\), \(37459\) and \(75658\) agree with the C++ program, and each is \(7\) times its number of splits plus \(2\), the splits being exactly the mixed blocks.
  • A Python version of the recursion makes exactly \(4M+1\) calls for every \(N\le14\).
  • A method without recursion and without the test of Section 6 counts the mixed blocks strip by strip: the black pixels of every row form one interval, and a block is black if its columns lie in the intervals of all its rows and white if they meet none of them. It agrees with the C++ program for every \(N\le24\), and for \(N\le22\) every \(M_k\) stays below \(45\) percent of the bound of Section 8.
  • For \(N=24\), the C++ program built with clang++ and with g++, the Java program and the Python bridge print the same value.

Footnotes and References

The references supply the background; the length formula \(7M+2\), the two-pixel test, the bound on the mixed blocks and all counts are derived or computed in this article.

  1. Project Euler 287 — Quadtree Encoding (a Simple Compression Algorithm). The official statement: 0 splits a region into its top left, top right, bottom left and bottom right sub-regions, 10 and 11 mark all-black and all-white regions, a 4 × 4 example has descriptions of lengths 30 and 16, and the task asks for the length of the minimal sequence for the disk image D_24.
  2. Wikipedia — Quadtree. A region quadtree divides a region into four equal quadrants, subquadrants and so on; in a binary image a region that is not entirely 0s or 1s is subdivided, so every leaf is a block of one value.
  3. Wikipedia — Tree traversal. Pre-order visits a node before traversing its subtrees, so every parent comes before its children; a pre-order sequence alone does not describe the underlying tree uniquely.
  4. cppreference — std::clamp. Since C++17, std::clamp(v, lo, hi) returns lo if v is less than lo, hi if hi is less than v, and v otherwise.
  5. MathWorld — Gauss's Circle Problem. The number N(r) of lattice points on or inside a circle of radius r centred at the origin is πr² + E(r), and Gauss showed that |E(r)| ≤ 2√2 πr for r ≥ 1.
  6. C++ — Euler287.cpp. The linked immutable revision contains parse_arguments, encode_length_fast, encode_length_bruteforce and run_checkpoints; it is the implementation record for this article.
  7. Wikipedia — Summed-area table. Every entry holds the sum of all values above and to the left of it, inclusive, and the sum over any rectangle follows from four entries.
  8. Python — subprocess. check_output runs a command, waits for it and returns its output, and it raises CalledProcessError if the return code is non-zero.

Problem 287 source code

C++

#include <algorithm>
#include <cstdint>
#include <iostream>
#include <string>
#include <vector>

namespace {

using u64 = std::uint64_t;

struct Options {
    int n = 24;
    bool run_checkpoints = true;
};

bool parse_int_after_prefix(const std::string& arg, const std::string& prefix, int& value) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }
    int parsed = 0;
    for (char c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }
        parsed = parsed * 10 + static_cast<int>(c - '0');
    }
    value = parsed;
    return true;
}

bool parse_arguments(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 (parse_int_after_prefix(arg, "--n=", options.n)) {
            continue;
        }
        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return options.n >= 1 && options.n <= 30;
}

u64 encode_length_fast(const int n) {
    const int size = 1 << n;
    const int c = 1 << (n - 1);
    const u64 r2 = 1ULL << (2 * n - 2);

    const auto rec = [&](auto&& self, const int x0, const int y0, const int side) -> u64 {
        const int x1 = x0 + side - 1;
        const int y1 = y0 + side - 1;

        const int nearest_x = std::clamp(c, x0, x1);
        const int nearest_y = std::clamp(c, y0, y1);
        const long long dx_min = static_cast<long long>(nearest_x) - c;
        const long long dy_min = static_cast<long long>(nearest_y) - c;
        const u64 min_d2 = static_cast<u64>(dx_min * dx_min + dy_min * dy_min);

        const long long dx_max = std::max(std::llabs(static_cast<long long>(x0) - c),
                                          std::llabs(static_cast<long long>(x1) - c));
        const long long dy_max = std::max(std::llabs(static_cast<long long>(y0) - c),
                                          std::llabs(static_cast<long long>(y1) - c));
        const u64 max_d2 = static_cast<u64>(dx_max * dx_max + dy_max * dy_max);

        if (max_d2 <= r2 || min_d2 > r2) {
            return 2ULL;
        }

        const int half = side / 2;
        return 1ULL + self(self, x0, y0 + half, half) + self(self, x0 + half, y0 + half, half) +
               self(self, x0, y0, half) + self(self, x0 + half, y0, half);
    };

    return rec(rec, 0, 0, size);
}

u64 encode_length_bruteforce(const int n) {
    const int size = 1 << n;
    const int c = 1 << (n - 1);
    const u64 r2 = 1ULL << (2 * n - 2);

    std::vector<std::vector<int>> black(static_cast<std::size_t>(size),
                                        std::vector<int>(static_cast<std::size_t>(size), 0));
    for (int y = 0; y < size; ++y) {
        for (int x = 0; x < size; ++x) {
            const long long dx = static_cast<long long>(x) - c;
            const long long dy = static_cast<long long>(y) - c;
            black[static_cast<std::size_t>(y)][static_cast<std::size_t>(x)] =
                (static_cast<u64>(dx * dx + dy * dy) <= r2) ? 1 : 0;
        }
    }

    std::vector<std::vector<int>> ps(static_cast<std::size_t>(size + 1),
                                     std::vector<int>(static_cast<std::size_t>(size + 1), 0));
    for (int y = 0; y < size; ++y) {
        int row_sum = 0;
        for (int x = 0; x < size; ++x) {
            row_sum += black[static_cast<std::size_t>(y)][static_cast<std::size_t>(x)];
            ps[static_cast<std::size_t>(y + 1)][static_cast<std::size_t>(x + 1)] =
                ps[static_cast<std::size_t>(y)][static_cast<std::size_t>(x + 1)] + row_sum;
        }
    }

    const auto sum_black = [&](const int x0, const int y0, const int side) {
        const int x1 = x0 + side;
        const int y1 = y0 + side;
        return ps[static_cast<std::size_t>(y1)][static_cast<std::size_t>(x1)] -
               ps[static_cast<std::size_t>(y0)][static_cast<std::size_t>(x1)] -
               ps[static_cast<std::size_t>(y1)][static_cast<std::size_t>(x0)] +
               ps[static_cast<std::size_t>(y0)][static_cast<std::size_t>(x0)];
    };

    const auto rec = [&](auto&& self, const int x0, const int y0, const int side) -> u64 {
        const int total = side * side;
        const int blacks = sum_black(x0, y0, side);
        if (blacks == 0 || blacks == total) {
            return 2ULL;
        }
        const int half = side / 2;
        return 1ULL + self(self, x0, y0 + half, half) + self(self, x0 + half, y0 + half, half) +
               self(self, x0, y0, half) + self(self, x0 + half, y0, half);
    };

    return rec(rec, 0, 0, size);
}

bool run_checkpoints() {
    for (int n = 1; n <= 5; ++n) {
        if (encode_length_fast(n) != encode_length_bruteforce(n)) {
            std::cerr << "Checkpoint failed for N=" << n << '\n';
            return false;
        }
    }
    return true;
}

}  // namespace

int main(int argc, char** argv) {
    Options options;
    if (!parse_arguments(argc, argv, options)) {
        return 1;
    }
    if (options.run_checkpoints && !run_checkpoints()) {
        return 2;
    }
    std::cout << encode_length_fast(options.n) << '\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 ""

    answer_candidates = []
    equal_candidates = []
    for line in lines:
        m1 = ANSWER_RE.search(line)
        if m1:
            answer_candidates.append(m1.group(1).strip())
        m2 = EQUAL_RE.search(line)
        if m2:
            equal_candidates.append(m2.group(1).strip())

    if answer_candidates:
        return answer_candidates[-1]
    if equal_candidates:
        return equal_candidates[-1]
    return lines[-1]


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 = subprocess.check_output([str(binary)], text=True)
    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

public class Euler287 {
    static long encodeLengthFast(int n) {
        int size = 1 << n;
        long c = 1L << (n - 1);
        long r2 = 1L << (2 * n - 2);

        return rec(0, 0, size, c, r2);
    }

    static long rec(long x0, long y0, long side, long c, long r2) {
        long x1 = x0 + side - 1;
        long y1 = y0 + side - 1;

        long nearestX = Math.max(x0, Math.min(c, x1));
        long nearestY = Math.max(y0, Math.min(c, y1));

        long dxMin = nearestX - c;
        long dyMin = nearestY - c;
        long minD2 = dxMin * dxMin + dyMin * dyMin;

        long dxMax = Math.max(Math.abs(x0 - c), Math.abs(x1 - c));
        long dyMax = Math.max(Math.abs(y0 - c), Math.abs(y1 - c));
        long maxD2 = dxMax * dxMax + dyMax * dyMax;

        if (maxD2 <= r2 || minD2 > r2) {
            return 2L;
        }

        long half = side / 2;
        return 1L + rec(x0, y0 + half, half, c, r2) +
                rec(x0 + half, y0 + half, half, c, r2) +
                rec(x0, y0, half, c, r2) +
                rec(x0 + half, y0, half, c, r2);
    }

    public static String solve() {
        int n = 24;
        return String.valueOf(encodeLengthFast(n));
    }

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