Problem 314: The Mouse on the Moon
View on Project EulerProject Euler Problem 314 Solution
EulerSolve provides an optimized solution for Project Euler Problem 314, The Mouse on the Moon, with C++, Python, Java, and a step-by-step mathematical explanation.
Problem Summary The wall is built on lattice posts inside a \(500\times 500\) square. The quantity to maximize is $$\frac{\text{enclosed area}}{\text{wall length}}.$$ By symmetry, the full wall can be reconstructed from one quarter of the boundary, namely a monotone chain from the top midpoint to the right midpoint. With $$G=250,$$ we only optimize a quarter-chain from $$ (0,G)\quad\text{to}\quad (G,0).$$ If this quarter-chain has area \(A\) under it and length \(L\), then the full symmetric wall has area \(4A\) and length \(4L\), so the global ratio is exactly $$\frac{4A}{4L}=\frac{A}{L}.$$ Mathematical Approach 1) Why one quarter is enough. The square and the objective are symmetric with respect to both coordinate axes through the center. Therefore an optimal wall can be assumed to be symmetric too. In the first quadrant, its boundary is a chain starting at the top midpoint \((0,G)\) and ending at the right midpoint \((G,0)\). We may also assume this chain is monotone: \(x\) never decreases and \(y\) never increases. Any backtracking would increase perimeter without helping the enclosed area. 2) Why the chain can be taken convex. If a local part of the chain bends the wrong way, replacing that dent by its chord increases or preserves the enclosed area while shortening the boundary. So an optimum must be convex in the quarter-square....
Detailed mathematical approach
Problem Summary
The wall is built on lattice posts inside a \(500\times 500\) square. The quantity to maximize is
$$\frac{\text{enclosed area}}{\text{wall length}}.$$
By symmetry, the full wall can be reconstructed from one quarter of the boundary, namely a monotone chain from the top midpoint to the right midpoint.
With
$$G=250,$$
we only optimize a quarter-chain from
$$ (0,G)\quad\text{to}\quad (G,0).$$
If this quarter-chain has area \(A\) under it and length \(L\), then the full symmetric wall has area \(4A\) and length \(4L\), so the global ratio is exactly
$$\frac{4A}{4L}=\frac{A}{L}.$$
Mathematical Approach
1) Why one quarter is enough.
The square and the objective are symmetric with respect to both coordinate axes through the center. Therefore an optimal wall can be assumed to be symmetric too. In the first quadrant, its boundary is a chain starting at the top midpoint \((0,G)\) and ending at the right midpoint \((G,0)\).
We may also assume this chain is monotone: \(x\) never decreases and \(y\) never increases. Any backtracking would increase perimeter without helping the enclosed area.
2) Why the chain can be taken convex.
If a local part of the chain bends the wrong way, replacing that dent by its chord increases or preserves the enclosed area while shortening the boundary. So an optimum must be convex in the quarter-square. Convexity means the slopes of successive segments are nondecreasing as we move from the top midpoint to the right midpoint.
3) Primitive step vectors are enough.
A segment from one lattice post to another with displacement
$$ (u,v),\qquad u,v\ge 0$$
passes through intermediate lattice posts whenever \(\gcd(u,v)>1\). Splitting such a segment at those intermediate posts does not change either its geometric length or the area under it. Therefore we only need primitive vectors
$$\gcd(u,v)=1.$$
The code includes the axis vectors \((1,0)\) and \((0,1)\) as well, so horizontal and vertical runs are still possible.
4) Coordinate transform used by the DP.
The actual quarter-chain lives in ordinary coordinates \((x,y)\) from \((0,G)\) to \((G,0)\). The implementation stores instead
$$j=G-y,$$
so the start point becomes \((0,0)\) and the endpoint becomes \((G,G)\). A geometric segment
$$ (x,y)\to(x+u,y-v)$$
therefore appears in the DP as
$$ (x,j)\to(x+u,j+v).$$
This is why the code fills `dp[x+u][j+v]` from `dp[x][j]`.
5) Exact area increment of one segment.
The area under one segment is the area of a trapezoid. If the segment goes from height \(y\) down to height \(y-v\) over horizontal distance \(u\), then
$$\Delta A=u\cdot\frac{y+(y-v)}{2}=u\left(y-\frac v2\right).$$
Since \(y=G-j\), this becomes
$$\Delta A=u\left(G-j-\frac v2\right).$$
This is exactly what the implementation computes as
$$u\left(G-\frac v2\right)-ju,$$
which explains the code lines `base_area = u * (GRID - 0.5 * v)` and the later `area_gain -= u` as \(j\) increases.
6) Exact length increment.
The boundary contribution of the segment is simply
$$\Delta L=\sqrt{u^2+v^2}.$$
So if a quarter-chain is described by a sequence of primitive vectors \((u_i,v_i)\), then
$$A=\sum_i \Delta A_i,\qquad L=\sum_i \sqrt{u_i^2+v_i^2}.$$
7) Fractional objective and Dinkelbach transform.
The target ratio is
$$\rho^*=\max \frac{A}{L}.$$
Instead of maximizing a ratio directly, fix a real parameter \(\lambda\) and maximize
$$F_\lambda=A-\lambda L.$$
For the optimal ratio \(\rho^*\), the maximum value of \(F_{\rho^*}\) is exactly \(0\). This is the standard Dinkelbach principle for fractional programming.
So the algorithm repeatedly solves the easier DP problem for a fixed \(\lambda\), obtains the best chain \((A,L)\), and updates
$$\lambda\leftarrow \frac{A}{L}.$$
When the update stops changing, we have reached the optimal ratio.
8) DP state and transition.
Let `score[x][j]` be the best value of
$$A-\lambda L$$
for a convex monotone chain ending at transformed point \((x,j)\), using only directions processed so far.
The vectors are sorted by increasing slope
$$\frac vu,$$
with \((1,0)\) first and \((0,1)\) last. This enforces nondecreasing slope and therefore convexity.
For one vector \((u,v)\), the transition is
$$dp[x+u,j+v]=\max\Bigl(dp[x+u,j+v],\ dp[x,j]+\Delta A-\lambda\Delta L\Bigr).$$
Because the update is performed in-place while the current vector is active, the same direction may be used repeatedly. That is exactly what we want for long straight portions of the chain.
9) Why sorting by slope avoids self-intersections.
All chosen segments move rightward and downward in actual coordinates, and their slopes are processed in nondecreasing order. This means the tangent direction of the chain only turns one way. Combined with monotonicity, the quarter-boundary stays convex and cannot self-intersect.
Worked Geometric Interpretation
Suppose the quarter-chain contains a single diagonal step from \((175,250)\) to \((250,175)\), together with the horizontal and vertical parts needed to connect \((0,250)\) to \((250,0)\). This is exactly the kind of “cut off a corner” example described in the problem statement. The DP can represent that chain using repeated \((1,0)\), then one copy of \((75,75)\) split into primitive slope-\(1\) steps, then repeated \((0,1)\).
The area under that chain is the quarter-area of the final enclosed wall, and the quarter-length is the corresponding quarter of the full boundary. The optimization simply searches all such convex monotone quarter-chains and chooses the one with maximum \(A/L\).
Algorithm
1) Build all primitive vectors \((u,v)\) with \(0\le u,v\le G\), plus \((1,0)\) and \((0,1)\).
2) Sort them by slope \(v/u\).
3) For a fixed \(\lambda\), run a DP on the \((G+1)\times(G+1)\) grid of transformed endpoints.
4) Recover the optimal quarter-area \(A\) and quarter-length \(L\) at state \((G,G)\).
5) Update \(\lambda\leftarrow A/L\) and repeat until convergence.
Complexity Analysis
Let \(K\) be the number of primitive directions. One DP solve costs
$$O(KG^2)$$
time and
$$O(G^2)$$
memory. The outer Dinkelbach iteration converges in only a small number of rounds, so the full computation is practical.
Checks And Final Result
The C++ program converges to
$$132.52756426,$$
which is the required maximum area-to-perimeter ratio rounded to eight decimal places.
Further Reading
- Problem page: https://projecteuler.net/problem=314
- Fractional programming / Dinkelbach method: https://en.wikipedia.org/wiki/Fractional_programming
- Primitive lattice vectors: https://en.wikipedia.org/wiki/Coprime_integers
Problem 314 source code
C++
#include <algorithm>
#include <cmath>
#include <iomanip>
#include <iostream>
#include <numeric>
#include <vector>
namespace {
constexpr int GRID = 250;
constexpr int DIM = GRID + 1;
constexpr double kNeg = -1e100;
constexpr double kEps = 1e-12;
struct Vec {
int u;
int v;
double len;
double slope;
};
int idx(int x, int y) { return x * DIM + y; }
std::vector<Vec> build_vectors() {
std::vector<Vec> vecs;
vecs.reserve(40000);
vecs.push_back({1, 0, 1.0, 0.0});
for (int u = 1; u <= GRID; ++u) {
for (int v = 1; v <= GRID; ++v) {
if (std::gcd(u, v) == 1) {
double len = std::sqrt(static_cast<double>(u) * u +
static_cast<double>(v) * v);
double slope = static_cast<double>(v) / u;
vecs.push_back({u, v, len, slope});
}
}
}
vecs.push_back({0, 1, 1.0, 1e100});
std::sort(vecs.begin(), vecs.end(),
[](const Vec &a, const Vec &b) { return a.slope < b.slope; });
return vecs;
}
struct Result {
double area;
double len;
double score;
};
Result solve_for_lambda(double lambda, const std::vector<Vec> &vecs) {
static std::vector<double> score(DIM * DIM);
static std::vector<double> area(DIM * DIM);
static std::vector<double> length(DIM * DIM);
std::fill(score.begin(), score.end(), kNeg);
std::fill(area.begin(), area.end(), 0.0);
std::fill(length.begin(), length.end(), 0.0);
score[idx(0, 0)] = 0.0;
for (const auto &vec : vecs) {
int u = vec.u;
int v = vec.v;
double len = vec.len;
double base_area = static_cast<double>(u) * (GRID - 0.5 * v);
double base_score = base_area - lambda * len;
for (int i = 0; i <= GRID - u; ++i) {
double *score_row = &score[idx(i, 0)];
double *area_row = &area[idx(i, 0)];
double *len_row = &length[idx(i, 0)];
double *score_row_to = &score[idx(i + u, 0)];
double *area_row_to = &area[idx(i + u, 0)];
double *len_row_to = &length[idx(i + u, 0)];
double area_gain = base_area;
double score_gain = base_score;
for (int j = 0; j <= GRID - v; ++j) {
double cur_score = score_row[j];
if (cur_score > kNeg / 2) {
int tj = j + v;
double cand_score = cur_score + score_gain;
if (cand_score > score_row_to[tj] + kEps) {
score_row_to[tj] = cand_score;
area_row_to[tj] = area_row[j] + area_gain;
len_row_to[tj] = len_row[j] + len;
}
}
area_gain -= u;
score_gain -= u;
}
}
}
int end = idx(GRID, GRID);
return {area[end], length[end], score[end]};
}
} // namespace
int main() {
auto vecs = build_vectors();
double lambda = 132.0;
for (int iter = 0; iter < 40; ++iter) {
Result res = solve_for_lambda(lambda, vecs);
double new_lambda = res.area / res.len;
if (std::abs(new_lambda - lambda) < 1e-12 ||
std::abs(res.score) < 1e-11) {
lambda = new_lambda;
break;
}
lambda = new_lambda;
}
std::cout << std::fixed << std::setprecision(8) << lambda << '\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 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
import java.util.*;
public class Euler314 {
static final int GRID = 250;
static final int DIM = GRID + 1;
static final double K_NEG = -1e100;
static final double K_EPS = 1e-12;
static class Vec implements Comparable<Vec> {
int u, v;
double len, slope;
Vec(int u, int v, double len, double slope) {
this.u = u;
this.v = v;
this.len = len;
this.slope = slope;
}
public int compareTo(Vec o) {
return Double.compare(this.slope, o.slope);
}
}
static int gcd(int a, int b) {
while (b != 0) {
int temp = b;
b = a % b;
a = temp;
}
return a;
}
static List<Vec> buildVectors() {
List<Vec> vecs = new ArrayList<>();
vecs.add(new Vec(1, 0, 1.0, 0.0));
for (int u = 1; u <= GRID; ++u) {
for (int v = 1; v <= GRID; ++v) {
if (gcd(u, v) == 1) {
double len = Math.sqrt((double) u * u + (double) v * v);
double slope = (double) v / u;
vecs.add(new Vec(u, v, len, slope));
}
}
}
vecs.add(new Vec(0, 1, 1.0, 1e100));
Collections.sort(vecs);
return vecs;
}
static class Result {
double area, len, score;
Result(double a, double l, double s) {
area = a;
len = l;
score = s;
}
}
static double[] score = new double[DIM * DIM];
static double[] area = new double[DIM * DIM];
static double[] length = new double[DIM * DIM];
static double[] rowMax = new double[DIM];
static Result solveForLambda(double lambda, List<Vec> vecs) {
Arrays.fill(score, K_NEG);
Arrays.fill(area, 0.0);
Arrays.fill(length, 0.0);
Arrays.fill(rowMax, K_NEG);
score[0] = 0.0;
rowMax[0] = 0.0;
for (int vi = 0; vi < vecs.size(); vi++) {
Vec vec = vecs.get(vi);
int u = vec.u;
int v = vec.v;
double len = vec.len;
double baseArea = (double) u * (GRID - 0.5 * v);
double baseScore = baseArea - lambda * len;
for (int i = 0; i <= GRID - u; ++i) {
if (rowMax[i] <= K_NEG / 2)
continue;
int idxFrom = i * DIM;
int idxTo = (i + u) * DIM;
double areaGain = baseArea;
double scoreGain = baseScore;
for (int j = 0; j <= GRID - v; ++j) {
double curScore = score[idxFrom + j];
if (curScore > K_NEG / 2) {
int tj = idxTo + j + v;
double candScore = curScore + scoreGain;
if (candScore > score[tj] + K_EPS) {
score[tj] = candScore;
area[tj] = area[idxFrom + j] + areaGain;
length[tj] = length[idxFrom + j] + len;
if (candScore > rowMax[i + u]) {
rowMax[i + u] = candScore;
}
}
}
areaGain -= u;
scoreGain -= u;
}
}
}
int end = GRID * DIM + GRID;
return new Result(area[end], length[end], score[end]);
}
public static String solve() {
List<Vec> vecs = buildVectors();
double lambda = 132.0;
for (int iter = 0; iter < 40; ++iter) {
Result res = solveForLambda(lambda, vecs);
double newLambda = res.area / res.len;
if (Math.abs(newLambda - lambda) < 1e-12 || Math.abs(res.score) < 1e-11) {
lambda = newLambda;
break;
}
lambda = newLambda;
}
return String.format(Locale.US, "%.8f", lambda);
}
public static void main(String[] args) {
System.out.println(solve());
}
}