Problem 695: Random Rectangles
View on Project EulerProject Euler Problem 695 Solution
EulerSolve provides an optimized solution for Project Euler Problem 695, Random Rectangles, with C++, Python, Java, and a step-by-step mathematical explanation.
Problem Summary The implementations evaluate an expected median area arising from a random rectangle construction. After normalizing two complementary horizontal gaps and two complementary vertical gaps, the geometry depends only on two ratios \(x,y\in[0,1]\). The target value becomes a two-dimensional integral of a normalized kernel, followed by a constant rescaling by \(1/4\). Mathematical Approach Introduce the complementary fractions $$x_0=1-x,\qquad y_0=1-y.$$ The normalized geometry is completely determined by these four numbers. Each admissible relative arrangement contributes the median of three candidate areas, and the kernel is the average over six equally weighted arrangements. Step 1: Express the Median of Three Areas For any triple \((a,b,c)\), the median can be written without sorting as $$\operatorname{med}(a,b,c)=a+b+c-\min(a,b,c)-\max(a,b,c).$$ This identity is exactly what the implementation uses. It avoids a separate ordering step and gives the middle value directly....
Detailed mathematical approach
Problem Summary
The implementations evaluate an expected median area arising from a random rectangle construction. After normalizing two complementary horizontal gaps and two complementary vertical gaps, the geometry depends only on two ratios \(x,y\in[0,1]\). The target value becomes a two-dimensional integral of a normalized kernel, followed by a constant rescaling by \(1/4\).
Mathematical Approach
Introduce the complementary fractions
$$x_0=1-x,\qquad y_0=1-y.$$
The normalized geometry is completely determined by these four numbers. Each admissible relative arrangement contributes the median of three candidate areas, and the kernel is the average over six equally weighted arrangements.
Step 1: Express the Median of Three Areas
For any triple \((a,b,c)\), the median can be written without sorting as
$$\operatorname{med}(a,b,c)=a+b+c-\min(a,b,c)-\max(a,b,c).$$
This identity is exactly what the implementation uses. It avoids a separate ordering step and gives the middle value directly.
Step 2: Build the Six Normalized Configurations
From the normalized side fractions, the six median terms are
$$\begin{aligned} T_1(x,y)&=\operatorname{med}(xy,\ x_0y_0,\ 1),\\ T_2(x,y)&=\operatorname{med}(xy,\ x_0,\ y_0),\\ T_3(x,y)&=\operatorname{med}(xy_0,\ x_0y,\ 1),\\ T_4(x,y)&=\operatorname{med}(xy_0,\ x_0,\ y),\\ T_5(x,y)&=\operatorname{med}(x,\ x_0y,\ y_0),\\ T_6(x,y)&=\operatorname{med}(x,\ x_0y_0,\ y). \end{aligned}$$
The normalized kernel is their average:
$$M(x,y)=\frac{T_1(x,y)+T_2(x,y)+T_3(x,y)+T_4(x,y)+T_5(x,y)+T_6(x,y)}{6}.$$
This formula is the mathematical core of the solution: once \(M(x,y)\) is available, the rest of the task is numerical integration.
Step 3: Exploit the Symmetries
Replacing \(x\) with \(1-x\), replacing \(y\) with \(1-y\), or swapping \(x\) and \(y\) only permutes the six triples above. Therefore
$$M(x,y)=M(1-x,y)=M(x,1-y)=M(y,x).$$
These identities are important for two reasons. First, they confirm that the six-term average has the expected geometric symmetry. Second, they justify computing the double quadrature with symmetry in the node indices.
Step 4: Turn the Expectation into an Integral
The normalized mean contribution is
$$J=\int_0^1\int_0^1 M(x,y)\,dx\,dy.$$
After the geometric normalization, the original expected area is obtained by multiplying by the scale factor left outside the ratio variables. In the implemented derivation that factor is \(1/4\), so the final quantity is
$$E=\frac{J}{4}.$$
The entire computational problem is therefore reduced to evaluating \(J\) very accurately.
Step 5: Approximate the Integral with Gauss-Legendre Quadrature
Let \((\xi_i,w_i)_{i=1}^n\) be the \(n\)-point Gauss-Legendre nodes and weights on \([0,1]\). Then
$$J_n=\sum_{i=1}^{n}\sum_{j=1}^{n} w_i w_j M(\xi_i,\xi_j)$$
approximates \(J\). Because the kernel is symmetric in its two arguments, the implementations reorganize the sum as
$$J_n=\sum_{i=1}^{n} w_i^2 M(\xi_i,\xi_i)+2\sum_{1\le i\lt j\le n} w_i w_j M(\xi_i,\xi_j).$$
This keeps the same quadrature rule while cutting the off-diagonal work almost in half.
Step 6: Cancel the Leading Error with Richardson Extrapolation
The kernel is smooth inside regions separated by median-switching boundaries, so the global quadrature error behaves like a piecewise-smooth integral rather than a fully analytic one. The implementation models the leading term as
$$E_n=E+\frac{C}{n^2}+O\left(\frac{1}{n^4}\right).$$
Computing once with \(n\) nodes and once with \(2n\) nodes gives
$$E_{2n}=E+\frac{C}{4n^2}+O\left(\frac{1}{n^4}\right).$$
Eliminating the unknown constant \(C\) yields
$$E_{\text{rich}}=E_{2n}+\frac{E_{2n}-E_n}{3}.$$
This Richardson-refined value is the reported answer.
Worked Example: Evaluate the Kernel at \(x=0\), \(y=\frac{1}{4}\)
Here \(x_0=1\) and \(y_0=\frac{3}{4}\). The six median terms become
$$\begin{aligned} T_1&=\operatorname{med}\left(0,\frac{3}{4},1\right)=\frac{3}{4},\\ T_2&=\operatorname{med}\left(0,1,\frac{3}{4}\right)=\frac{3}{4},\\ T_3&=\operatorname{med}\left(0,\frac{1}{4},1\right)=\frac{1}{4},\\ T_4&=\operatorname{med}\left(0,1,\frac{1}{4}\right)=\frac{1}{4},\\ T_5&=\operatorname{med}\left(0,\frac{1}{4},\frac{3}{4}\right)=\frac{1}{4},\\ T_6&=\operatorname{med}\left(0,\frac{3}{4},\frac{1}{4}\right)=\frac{1}{4}. \end{aligned}$$
Therefore
$$M\left(0,\frac{1}{4}\right)=\frac{\frac{3}{4}+\frac{3}{4}+\frac{1}{4}+\frac{1}{4}+\frac{1}{4}+\frac{1}{4}}{6}=\frac{5}{12},$$
which is one of the checkpoint values verified by the implementation.
How the Code Works
The C++, Python, and Java implementations all follow the same numerical plan. They evaluate the six median expressions for any requested pair \((x,y)\), generate Gauss-Legendre nodes and weights on \([0,1]\) from Legendre roots, and accumulate the symmetric double sum for two resolutions, \(n\) and \(2n\).
Each integral approximation is divided by \(4\), then the two results are combined with Richardson extrapolation. The C++ and Java implementations parallelize the outer loop over quadrature rows, while the Python implementation delegates to the compiled numerical solver and returns the parsed final value.
Before reporting the answer, the workflow checks fixed kernel values, symmetry identities, and benchmark quadrature outputs at moderate orders. Those checks guard against both algebraic mistakes in the kernel and numerical mistakes in the quadrature generator.
Complexity Analysis
For order \(n\), generating the Gauss-Legendre nodes and weights costs \(O(n^2)\) arithmetic overall because each root evaluation uses a degree-\(n\) recurrence and there are \(O(n)\) roots. The dominant cost is the double quadrature itself, which requires \(O(n^2)\) kernel evaluations and \(O(n)\) memory for the node and weight arrays. Because the sum is split by rows, the wall-clock time benefits directly from multi-core parallelism.
Footnotes and References
- Problem page: https://projecteuler.net/problem=695
- Gaussian quadrature: Wikipedia — Gaussian quadrature
- Legendre polynomials: Wikipedia — Legendre polynomials
- Richardson extrapolation: Wikipedia — Richardson extrapolation
- Median: Wikipedia — Median
Problem 695 source code
C++
#include <algorithm>
#include <array>
#include <cerrno>
#include <cmath>
#include <cstdint>
#include <cstdlib>
#include <iomanip>
#include <iostream>
#include <limits>
#include <string>
#include <thread>
#include <vector>
namespace {
using u32 = std::uint32_t;
constexpr int kDefaultBaseN = 8192;
struct Options {
int base_n = kDefaultBaseN;
bool run_checkpoints = true;
bool allow_multithreading = true;
unsigned requested_threads = 0U;
};
inline double median3(const double a, const double b, const double c) {
return a + b + c - std::min(a, std::min(b, c)) - std::max(a, std::max(b, c));
}
// M(x, y): expectation over the 6 relative label permutations of the median area,
// after normalizing gaps so that a+b=1 and c+d=1.
inline double mean_median_normalized(const double x, const double y) {
const double x0 = 1.0 - x;
const double y0 = 1.0 - y;
const double m1 = median3(x * y, x0 * y0, 1.0);
const double m2 = median3(x * y, x0, y0);
const double m3 = median3(x * y0, x0 * y, 1.0);
const double m4 = median3(x * y0, x0, y);
const double m5 = median3(x, x0 * y, y0);
const double m6 = median3(x, x0 * y0, y);
return (m1 + m2 + m3 + m4 + m5 + m6) * (1.0 / 6.0);
}
bool parse_u32_after_prefix(const std::string& arg, const char* prefix, u32& value) {
const std::string p(prefix);
if (arg.rfind(p, 0) != 0U) {
return false;
}
const std::string tail = arg.substr(p.size());
if (tail.empty()) {
return false;
}
std::uint64_t parsed = 0ULL;
for (const char c : tail) {
if (c < '0' || c > '9') {
return false;
}
parsed = parsed * 10ULL + static_cast<std::uint64_t>(c - '0');
if (parsed > static_cast<std::uint64_t>(std::numeric_limits<u32>::max())) {
return false;
}
}
value = static_cast<u32>(parsed);
return true;
}
bool parse_unsigned_after_prefix(const std::string& arg,
const char* prefix,
unsigned& value) {
u32 parsed = 0U;
if (!parse_u32_after_prefix(arg, prefix, parsed)) {
return false;
}
value = static_cast<unsigned>(parsed);
return true;
}
bool parse_arguments(const 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 (arg == "--single-thread") {
options.allow_multithreading = false;
continue;
}
u32 parsed_u32 = 0U;
if (parse_u32_after_prefix(arg, "--n=", parsed_u32)) {
if (parsed_u32 > static_cast<u32>(std::numeric_limits<int>::max())) {
std::cerr << "--n is too large.\n";
return false;
}
options.base_n = static_cast<int>(parsed_u32);
continue;
}
unsigned parsed_unsigned = 0U;
if (parse_unsigned_after_prefix(arg, "--threads=", parsed_unsigned)) {
options.requested_threads = parsed_unsigned;
continue;
}
std::cerr << "Unknown argument: " << arg << '\n';
return false;
}
if (options.base_n < 2) {
std::cerr << "--n must be at least 2.\n";
return false;
}
return true;
}
unsigned choose_thread_count(const bool allow_multithreading,
const unsigned requested_threads,
const int workload) {
if (!allow_multithreading || workload < 2) {
return 1U;
}
unsigned threads = requested_threads;
if (threads == 0U) {
threads = std::thread::hardware_concurrency();
if (threads == 0U) {
threads = 1U;
}
}
return std::max(1U, std::min<unsigned>(threads, static_cast<unsigned>(workload)));
}
struct NodesWeights {
std::vector<double> x;
std::vector<double> w;
};
NodesWeights gauss_legendre_unit_interval(const int n) {
NodesWeights nw;
nw.x.assign(static_cast<std::size_t>(n), 0.0);
nw.w.assign(static_cast<std::size_t>(n), 0.0);
const double pi = std::acos(-1.0);
const double eps = 1e-15;
const int half = (n + 1) / 2;
for (int i = 0; i < half; ++i) {
double z = std::cos(pi * (static_cast<double>(i) + 0.75) / (static_cast<double>(n) + 0.5));
double z_prev = 0.0;
double p1 = 0.0;
double p2 = 0.0;
double pp = 0.0;
do {
p1 = 1.0;
p2 = 0.0;
for (int j = 1; j <= n; ++j) {
const double p3 = p2;
p2 = p1;
p1 = ((2.0 * static_cast<double>(j) - 1.0) * z * p2 -
(static_cast<double>(j) - 1.0) * p3) /
static_cast<double>(j);
}
pp = static_cast<double>(n) * (z * p1 - p2) / (z * z - 1.0);
z_prev = z;
z = z_prev - p1 / pp;
} while (std::abs(z - z_prev) > eps);
const double weight = 2.0 / ((1.0 - z * z) * pp * pp);
const int i_left = i;
const int i_right = n - 1 - i;
nw.x[static_cast<std::size_t>(i_left)] = 0.5 * (-z + 1.0);
nw.x[static_cast<std::size_t>(i_right)] = 0.5 * (z + 1.0);
nw.w[static_cast<std::size_t>(i_left)] = 0.5 * weight;
nw.w[static_cast<std::size_t>(i_right)] = 0.5 * weight;
}
return nw;
}
double integrate_j(const int n,
const bool allow_multithreading,
const unsigned requested_threads) {
const NodesWeights nw = gauss_legendre_unit_interval(n);
const unsigned threads = choose_thread_count(allow_multithreading, requested_threads, n);
std::vector<std::thread> pool;
pool.reserve(threads);
std::vector<double> partial(threads, 0.0);
std::atomic<int> next_i(0);
for (unsigned t = 0; t < threads; ++t) {
pool.emplace_back([&, t]() {
double local_sum = 0.0;
while (true) {
const int i = next_i.fetch_add(1, std::memory_order_relaxed);
if (i >= n) {
break;
}
const double xi = nw.x[static_cast<std::size_t>(i)];
const double wi = nw.w[static_cast<std::size_t>(i)];
double row = wi * wi * mean_median_normalized(xi, xi);
for (int j = i + 1; j < n; ++j) {
row += 2.0 * wi * nw.w[static_cast<std::size_t>(j)] *
mean_median_normalized(xi, nw.x[static_cast<std::size_t>(j)]);
}
local_sum += row;
}
partial[static_cast<std::size_t>(t)] = local_sum;
});
}
for (auto& th : pool) {
th.join();
}
double total = 0.0;
for (const double s : partial) {
total += s;
}
return total;
}
bool nearly_equal(const double a, const double b, const double eps) {
return std::abs(a - b) <= eps;
}
bool run_checkpoints() {
{
const double c1 = mean_median_normalized(0.0, 0.0);
const double c2 = mean_median_normalized(0.0, 0.5);
const double c3 = mean_median_normalized(0.25, 0.25);
const double c4 = mean_median_normalized(0.0, 0.25);
if (!nearly_equal(c1, 1.0 / 3.0, 1e-15)) {
std::cerr << "Checkpoint failed: M(0,0) != 1/3.\n";
return false;
}
if (!nearly_equal(c2, 0.5, 1e-15)) {
std::cerr << "Checkpoint failed: M(0,0.5) != 1/2.\n";
return false;
}
if (!nearly_equal(c3, 0.375, 1e-15)) {
std::cerr << "Checkpoint failed: M(0.25,0.25) != 3/8.\n";
return false;
}
if (!nearly_equal(c4, 5.0 / 12.0, 1e-15)) {
std::cerr << "Checkpoint failed: M(0,0.25) != 5/12.\n";
return false;
}
}
{
constexpr std::array<std::pair<double, double>, 4> points = {
std::pair<double, double>{0.19, 0.73},
std::pair<double, double>{0.41, 0.22},
std::pair<double, double>{0.8, 0.3},
std::pair<double, double>{0.63, 0.63},
};
for (const auto& [x, y] : points) {
const double f = mean_median_normalized(x, y);
const double fx = mean_median_normalized(1.0 - x, y);
const double fy = mean_median_normalized(x, 1.0 - y);
const double fs = mean_median_normalized(y, x);
if (!nearly_equal(f, fx, 2e-15) || !nearly_equal(f, fy, 2e-15) ||
!nearly_equal(f, fs, 2e-15)) {
std::cerr << "Checkpoint failed: symmetry mismatch at x=" << x << ", y=" << y
<< ".\n";
return false;
}
}
}
{
const double e256 = integrate_j(256, false, 1U) / 4.0;
const double e512 = integrate_j(512, false, 1U) / 4.0;
const double e1024 = integrate_j(1024, false, 1U) / 4.0;
if (!nearly_equal(e256, 0.10177862574564227, 5e-14)) {
std::cerr << "Checkpoint failed: n=256 quadrature mismatch.\n";
return false;
}
if (!nearly_equal(e512, 0.10177867340598332, 5e-14)) {
std::cerr << "Checkpoint failed: n=512 quadrature mismatch.\n";
return false;
}
if (!nearly_equal(e1024, 0.10177868278503224, 5e-14)) {
std::cerr << "Checkpoint failed: n=1024 quadrature mismatch.\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 1;
}
const int base_n = options.base_n;
const int fine_n = base_n * 2;
const auto t0 = std::chrono::steady_clock::now();
const double e_base = integrate_j(base_n, options.allow_multithreading, options.requested_threads) / 4.0;
const double e_fine = integrate_j(fine_n, options.allow_multithreading, options.requested_threads) / 4.0;
// Leading quadrature error behaves like O(n^-2) because of piecewise-smooth boundaries.
const double e_richardson = e_fine + (e_fine - e_base) / 3.0;
const auto t1 = std::chrono::steady_clock::now();
const double elapsed_sec = std::chrono::duration<double>(t1 - t0).count();
std::cout << std::fixed << std::setprecision(15);
std::cout << "E_base(n=" << base_n << ") = " << e_base << '\n';
std::cout << "E_fine(n=" << fine_n << ") = " << e_fine << '\n';
std::cout << "E_richardson = " << e_richardson << '\n';
std::cout << "Answer (10 d.p.) = " << std::setprecision(10) << e_richardson << '\n';
std::cout << std::setprecision(6) << "Elapsed seconds = " << elapsed_sec << '\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_like = []
answers = []
equals = []
for line in lines:
lower = line.lower()
if "answer" in lower:
if ":" in line:
answer_like.append(line.rsplit(":", 1)[1].strip())
elif "=" in line:
answer_like.append(line.rsplit("=", 1)[1].strip())
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 answer_like:
return answer_like[-1]
if answers:
return answers[-1]
if equals:
return equals[-1]
if ":" in lines[-1]:
tail = lines[-1].rsplit(":", 1)[1].strip()
if tail:
return tail
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.List;
import java.util.concurrent.*;
import java.util.Locale;
public class Euler695 {
static double median3(double a, double b, double c) {
return a + b + c - Math.min(a, Math.min(b, c)) - Math.max(a, Math.max(b, c));
}
static double meanMedianNormalized(double x, double y) {
double x0 = 1.0 - x;
double y0 = 1.0 - y;
double m1 = median3(x * y, x0 * y0, 1.0);
double m2 = median3(x * y, x0, y0);
double m3 = median3(x * y0, x0 * y, 1.0);
double m4 = median3(x * y0, x0, y);
double m5 = median3(x, x0 * y, y0);
double m6 = median3(x, x0 * y0, y);
return (m1 + m2 + m3 + m4 + m5 + m6) / 6.0;
}
static class NodesWeights {
double[] x;
double[] w;
}
static NodesWeights gaussLegendreUnitInterval(int n) {
NodesWeights nw = new NodesWeights();
nw.x = new double[n];
nw.w = new double[n];
double pi = Math.PI;
double eps = 1e-15;
int half = (n + 1) / 2;
for (int i = 0; i < half; ++i) {
double z = Math.cos(pi * (i + 0.75) / (n + 0.5));
double zPrev = 0.0;
double p1 = 0.0, p2 = 0.0, pp = 0.0;
do {
p1 = 1.0;
p2 = 0.0;
for (int j = 1; j <= n; ++j) {
double p3 = p2;
p2 = p1;
p1 = ((2.0 * j - 1.0) * z * p2 - (j - 1.0) * p3) / j;
}
pp = n * (z * p1 - p2) / (z * z - 1.0);
zPrev = z;
z = zPrev - p1 / pp;
} while (Math.abs(z - zPrev) > eps);
double weight = 2.0 / ((1.0 - z * z) * pp * pp);
int iLeft = i;
int iRight = n - 1 - i;
nw.x[iLeft] = 0.5 * (-z + 1.0);
nw.x[iRight] = 0.5 * (z + 1.0);
nw.w[iLeft] = 0.5 * weight;
nw.w[iRight] = 0.5 * weight;
}
return nw;
}
static double integrateJ(int n) {
NodesWeights nw = gaussLegendreUnitInterval(n);
int threads = Runtime.getRuntime().availableProcessors();
if (threads < 1)
threads = 1;
ExecutorService pool = Executors.newFixedThreadPool(threads);
List<Future<Double>> futures = new ArrayList<>();
for (int t = 0; t < threads; ++t) {
int start = (n * t) / threads;
int end = (n * (t + 1)) / threads;
futures.add(pool.submit(() -> {
double localSum = 0.0;
for (int i = start; i < end; ++i) {
double xi = nw.x[i];
double wi = nw.w[i];
double row = wi * wi * meanMedianNormalized(xi, xi);
for (int j = i + 1; j < n; ++j) {
row += 2.0 * wi * nw.w[j] * meanMedianNormalized(xi, nw.x[j]);
}
localSum += row;
}
return localSum;
}));
}
double total = 0.0;
for (Future<Double> f : futures) {
try {
total += f.get();
} catch (Exception e) {
}
}
pool.shutdown();
return total;
}
public static String solve() {
int baseN = 8192;
int fineN = baseN * 2;
double eBase = integrateJ(baseN) / 4.0;
double eFine = integrateJ(fineN) / 4.0;
double eRichardson = eFine + (eFine - eBase) / 3.0;
return String.format(Locale.US, "%.10f", eRichardson);
}
public static void main(String[] args) {
System.out.println(solve());
}
}