Problem 287: Quadtree Encoding (a Simple Compression Algorithm)
View on Project EulerProject 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
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. The minimum of \(D(x,y)\) over the block is obtained by clamping the center coordinate into the interval of the block:
$$x_{\min}=\operatorname{clamp}(c,x_0,x_1),\qquad y_{\min}=\operatorname{clamp}(c,y_0,y_1).$$
Then
$$d_{\min}^2=(x_{\min}-c)^2+(y_{\min}-c)^2.$$
This works because \(D(x,y)\) is the sum of two independent convex quadratic functions, one in \(x\) and one in \(y\).
4) Exact formula for the farthest pixel. The maximum of \(D(x,y)\) over an axis-aligned block is always attained at a corner. Therefore we can compute
$$dx_{\max}=\max(|x_0-c|,|x_1-c|),\qquad dy_{\max}=\max(|y_0-c|,|y_1-c|),$$
and then
$$d_{\max}^2=dx_{\max}^2+dy_{\max}^2.$$
This is exactly what the fast solver uses.
5) Uniformity test. With these two numbers, the recursion decision is immediate:
$$d_{\max}^2\le R^2 \Rightarrow \text{all black},$$
$$d_{\min}^2>R^2 \Rightarrow \text{all white},$$
otherwise the circle boundary crosses the block and the block must be subdivided.
6) Recursive length formula. If \(B\) is uniform,
$$L(B)=2.$$
If \(B\) is mixed and split into \(B_{NW},B_{NE},B_{SW},B_{SE}\), then
$$L(B)=1+L(B_{NW})+L(B_{NE})+L(B_{SW})+L(B_{SE}).$$
The child order matters for the bitstream but not for the total length. The code uses the order NW, NE, SW, SE.
Worked Examples
Example 1: \(N=1\). The image is \(2\times2\), with \(c=1\) and \(R^2=1\). Among the four pixels, \((1,1)\), \((1,0)\), and \((0,1)\) are black, while \((0,0)\) is white. So the root block is mixed and must split into four \(1\times1\) leaves:
$$L=1+4\cdot 2=9.$$
Example 2: small checkpoints. The fast method agrees with the brute-force quadtree on
$$N=1,2,3,4,5$$
and the corresponding lengths are
$$9,\ 30,\ 86,\ 212,\ 499.$$
These values are a strong sanity check for the geometric pruning logic.
Why the Fast Algorithm Is Correct
The brute-force method explicitly colors every pixel and then compresses. The fast method never builds the bitmap: it replaces “is this block monochrome?” by the exact min/max distance test above. Because those tests are necessary and sufficient, every block is classified exactly as in the brute-force version, so both traversals visit the same logical quadtree and produce the same total length.
Complexity Analysis
The running time is proportional to the number of visited quadtree nodes, not to the total number of pixels \(4^N\). Large regions fully inside or fully outside the disk terminate immediately; only blocks near the circle boundary recurse deeply. Memory usage is just the recursion depth:
$$O(N).$$
Further Reading
- Problem page: https://projecteuler.net/problem=287
- Quadtree: https://en.wikipedia.org/wiki/Quadtree
- Distance to an axis-aligned box: https://en.wikipedia.org/wiki/Euclidean_distance
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());
}
}