Problem 262: Mountain Range
View on Project EulerProject Euler Problem 262 Solution
EulerSolve provides an optimized solution for Project Euler Problem 262, Mountain Range, with C++, Python, Java, and a step-by-step mathematical explanation.
Problem Summary The mountain is described by the height function $$ h(x,y)=\Bigl(5000-0.005(x^2+y^2+xy)+12.5(x+y)\Bigr)e^{-\left|10^{-6}(x^2+y^2)-0.0015(x+y)+0.7\right|}. $$ We must travel from $$ A=(200,200)\qquad\text{to}\qquad B=(1400,1400) $$ while obeying the rule that the route may never rise above the smallest possible pass height. Among all such valid routes, we want the shortest one. The implementation therefore solves two problems in sequence: 1. find the minimum feasible ceiling height, 2. under that ceiling, find the shortest route. Mathematical Approach 1. The Minimax Pass Height For any continuous path \(\gamma\) from \(A\) to \(B\), define its ceiling cost by $$ C(\gamma)=\max_{p\in\gamma} h(p). $$ The smallest mountain pass height is then $$ f_*=\min_{\gamma:A\to B} C(\gamma). $$ This is not yet a shortest-path problem in the Euclidean sense. It is a bottleneck problem: we minimize the highest altitude touched by the route. Only after \(f_*\) is known does the actual geometric shortest-path problem begin. 2. Discrete Approximation on a Grid The code samples the square $$ [0,1600]\times[0,1600] $$ on a regular grid with spacing grid-step . With the default setting grid-step = 1.0 , this gives $$ n=1601 $$ grid points on each axis. At each grid point \((i\Delta,j\Delta)\), the program stores the sampled height $$ h(i\Delta,j\Delta)....
Detailed mathematical approach
Problem Summary
The mountain is described by the height function
$$ h(x,y)=\Bigl(5000-0.005(x^2+y^2+xy)+12.5(x+y)\Bigr)e^{-\left|10^{-6}(x^2+y^2)-0.0015(x+y)+0.7\right|}. $$
We must travel from
$$ A=(200,200)\qquad\text{to}\qquad B=(1400,1400) $$
while obeying the rule that the route may never rise above the smallest possible pass height. Among all such valid routes, we want the shortest one.
The implementation therefore solves two problems in sequence:
1. find the minimum feasible ceiling height,
2. under that ceiling, find the shortest route.
Mathematical Approach
1. The Minimax Pass Height
For any continuous path \(\gamma\) from \(A\) to \(B\), define its ceiling cost by
$$ C(\gamma)=\max_{p\in\gamma} h(p). $$
The smallest mountain pass height is then
$$ f_*=\min_{\gamma:A\to B} C(\gamma). $$
This is not yet a shortest-path problem in the Euclidean sense. It is a bottleneck problem: we minimize the highest altitude touched by the route.
Only after \(f_*\) is known does the actual geometric shortest-path problem begin.
2. Discrete Approximation on a Grid
The code samples the square
$$ [0,1600]\times[0,1600] $$
on a regular grid with spacing grid-step. With the default setting grid-step = 1.0, this gives
$$ n=1601 $$
grid points on each axis.
At each grid point \((i\Delta,j\Delta)\), the program stores the sampled height
$$ h(i\Delta,j\Delta). $$
This replaces the continuous terrain by a finite graph whose vertices are grid samples and whose edges connect the four axis-adjacent neighbors.
The four-neighbor choice matters because the code is using the grid only to estimate topological connectivity under a height ceiling; diagonal moves are not needed for the minimax search itself.
3. Why Dijkstra Still Works for a Bottleneck Objective
For a grid path
$$ v_0,v_1,\dots,v_t, $$
the discrete bottleneck cost is
$$ \max\bigl(h(v_0),h(v_1),\dots,h(v_t)\bigr). $$
If the current best ceiling at a vertex is \(c\), and we step to a neighbor \(u\), the new ceiling becomes
$$ c'=\max(c,h(u)). $$
This update is monotone: extending a path never lowers its ceiling. That is exactly the property Dijkstra needs. Once a vertex leaves the priority queue with the smallest currently known ceiling, no later path can improve it.
So the routine minimax_height is a correct Dijkstra-style solver for the minimax problem, with relaxation rule
$$ \text{best}[u]\leftarrow \min\bigl(\text{best}[u], \max(\text{best}[v],h(u))\bigr). $$
The returned value is a grid approximation of the true continuous pass height \(f_*\).
4. From Pass Height to a Forbidden Region
Once \(f_*\) is known, define the forbidden set
$$ \Omega=\{(x,y): h(x,y)>f_*\}. $$
The allowed set is its complement
$$ F=\{(x,y): h(x,y)\le f_*\}. $$
Any valid route must lie entirely in \(F\). Therefore the shortest valid route is simply the shortest Euclidean path from \(A\) to \(B\) in the plane with the obstacle region \(\Omega\) removed.
The critical boundary is the contour
$$ \partial\Omega=\{(x,y): h(x,y)=f_*\}. $$
The whole second phase of the algorithm is about approximating this contour and then solving a shortest-path problem around it.
5. Why the Final Route Reduces to “Segment + Arc + Segment”
For a fixed obstacle with smooth or polygonal boundary, a shortest path in the free plane has a standard structure:
1. whenever the path is strictly inside free space, it is a straight line segment,
2. if it cannot remain straight because the obstacle blocks it, it touches the obstacle boundary,
3. while constrained by the obstacle, it follows the boundary, then leaves it along another straight visible segment.
So once the relevant contour loop is known, the shortest valid route must be one of these forms:
1. the direct segment \(AB\), if it stays below \(f_*\),
2. \(A\to P_i\), then an arc of the contour, then \(P_j\to B\), where \(P_i\) and \(P_j\) are visible contour points.
This is the geometric reason the code only scans visible contour-point pairs instead of searching arbitrary curved paths in the plane.
6. Extracting the Level Set with Marching Squares
The grid already tells us which sampled vertices are above the ceiling and which are below it. The code marks a sample as inside the forbidden region when
$$ h(i\Delta,j\Delta)>f_*. $$
Each grid cell has four corners, so there are \(2^4=16\) possible above/below patterns. The arrays kCaseCount and kCaseEdges are the usual marching-squares lookup tables for those 16 cases.
If an edge crosses the contour \(h=f_*\), the crossing point is approximated by linear interpolation. If an edge endpoint heights are \(h_1\) and \(h_2\), then the crossing fraction is
$$ t=\frac{f_*-h_1}{h_2-h_1}, $$
and the point is
$$ P(t)=P_1+t(P_2-P_1). $$
Cells of ambiguous type produce two contour segments; ordinary crossing cells produce one. Running over all cells produces a polygonal approximation of the level set \(h=f_*\).
7. Quantization, Adjacency, and Closed Loops
The interpolated segment endpoints are floating-point numbers. Two segments that should meet can differ by tiny roundoff noise, so the code quantizes each point with scale
$$ 10^6. $$
That means nearly coincident points are snapped to the same integer key.
The program then builds an adjacency map from these quantized keys. In an ideal simple contour, every contour vertex has degree 2. The code checks exactly that.
Next it walks the adjacency graph to recover closed loops. If several loops appear, the implementation keeps the loop with largest perimeter.
This is an implementation-level robustness choice: numerically extracted contours can split into several pieces, but the outer main obstacle loop is the one relevant for the shortest route around the mountain.
8. Arc Lengths Along the Contour
Once a single contour loop
$$ P_0,P_1,\dots,P_{m-1} $$
has been recovered, the code computes prefix arc lengths
$$ S_0=0,\qquad S_{i+1}=S_i+|P_iP_{i+1}|, $$
with indices taken cyclically. The total contour perimeter is
$$ L=S_m. $$
For two contour vertices \(P_i\) and \(P_j\), one arc has length
$$ s_{ij}=|S_j-S_i|, $$
and the opposite arc has length
$$ L-s_{ij}. $$
The shorter one is the relevant boundary-travel cost for that pair.
9. Visibility Testing from \(A\) and \(B\)
A contour point \(P_i\) is usable from \(A\) only if the straight segment \(AP_i\) stays under the ceiling. Since the code does not solve the exact intersection analytically, it checks visibility by sampling points along the segment with step size vis-step.
If any sampled point satisfies
$$ h(x,y)>f_*+10^{-10}, $$
then the segment is rejected.
This produces two visible-index sets:
$$ V_A=\{i: A\text{ can see }P_i\},\qquad V_B=\{j: B\text{ can see }P_j\}. $$
Only pairs from \(V_A\times V_B\) need to be examined.
10. Final Candidate Formula
For each visible pair \((i,j)\), the candidate route length is
$$ |AP_i|+\min\bigl(s_{ij},\,L-s_{ij}\bigr)+|P_jB|. $$
The code also checks the direct segment \(AB\) first. If it is visible under \(f_*\), then no bent route can beat it, because a straight segment is the shortest Euclidean connection between two points.
Otherwise, the minimum over all visible contour pairs is the program's approximation to the true shortest valid route.
11. Implementation Sanity Checks
The program validates several necessary conditions:
1. \(f_*\) cannot be below the heights of \(A\) or \(B\),
2. the contour must not be empty,
3. every contour node must have degree 2,
4. the final route cannot be shorter than the direct distance \(|AB|\).
With the default parameters, the endpoint heights are approximately
$$ h(A)\approx 7851.540,\qquad h(B)\approx 6964.696. $$
The straight segment from \(A\) to \(B\) reaches heights above
$$ 13275, $$
so it is not feasible at the optimal ceiling.
When run with --verbose, the program reports
$$ f_*\approx 10396.458664, $$
and with default discretization it prints the route length
$$ 2531.205. $$
These are code-level numerical outputs, not symbolic closed forms.
How the Code Works
compute_heights samples the height field on the grid, optionally in parallel. minimax_height then runs the priority-queue minimax search and returns the discrete pass height.
After that, the code marks grid vertices with height above \(f_*\), extracts contour segments cell by cell with marching squares, quantizes segment endpoints, builds an adjacency graph, walks it into loops, chooses the largest loop, and computes prefix arc lengths around that loop.
Next, it evaluates which contour points are visible from \(A\) and \(B\), checks the direct segment \(AB\), and finally scans all visible pairs with the formula
$$ |AP_i|+\min\bigl(s_{ij},L-s_{ij}\bigr)+|P_jB|. $$
The minimum of those candidates is printed to three decimal places.
Complexity Analysis
If the grid has \(n\) points on each axis, then height sampling costs \(O(n^2)\), the minimax Dijkstra pass costs \(O(n^2\log n)\), contour extraction costs \(O(n^2)\), and visibility checks cost
$$ O(|V_A|+|V_B|)\times\text{sampling work} $$
plus the final pair scan
$$ O(|V_A||V_B|). $$
Memory usage is \(O(n^2)\) because the sampled height field dominates storage.
The important tradeoff is numerical: smaller grid-step and vis-step improve the approximation to the continuous problem, but they increase runtime and memory consumption.
Footnotes and References
- Problem page: https://projecteuler.net/problem=262
- Bottleneck / widest path problems: Wikipedia - widest path problem
- Marching squares: Wikipedia - marching squares
- Dijkstra's algorithm: Wikipedia - Dijkstra's algorithm
- Shortest path among obstacles / visibility ideas: Wikipedia - visibility graph
Problem 262 source code
C++
#include <algorithm>
#include <cmath>
#include <cstdint>
#include <cstdlib>
#include <iomanip>
#include <iostream>
#include <limits>
#include <queue>
#include <thread>
#include <unordered_map>
#include <unordered_set>
#include <vector>
namespace {
struct Point {
double x = 0.0;
double y = 0.0;
};
double height(double x, double y) {
const double xx = x * x;
const double yy = y * y;
const double xy = x * y;
const double base = 5000.0 - 0.005 * (xx + yy + xy) + 12.5 * (x + y);
const double expo = -std::abs(0.000001 * (xx + yy) - 0.0015 * (x + y) + 0.7);
return base * std::exp(expo);
}
double dist(const Point &a, const Point &b) {
return std::hypot(a.x - b.x, a.y - b.y);
}
struct Key {
int64_t x = 0;
int64_t y = 0;
bool operator==(const Key &other) const { return x == other.x && y == other.y; }
};
struct KeyHash {
size_t operator()(const Key &k) const {
uint64_t x = static_cast<uint64_t>(k.x);
uint64_t y = static_cast<uint64_t>(k.y);
return static_cast<size_t>((x * 0x9e3779b97f4a7c15ULL) ^ (y + (x << 6) + (x >> 2)));
}
};
Key quantize(const Point &p, double scale) {
return {static_cast<int64_t>(std::llround(p.x * scale)),
static_cast<int64_t>(std::llround(p.y * scale))};
}
Point from_key(const Key &k, double scale) {
return {static_cast<double>(k.x) / scale, static_cast<double>(k.y) / scale};
}
Point interp(const Point &p1, const Point &p2, double v1, double v2, double f) {
double t;
const double denom = v2 - v1;
if (std::abs(denom) < 1e-14) {
t = 0.5;
} else {
t = (f - v1) / denom;
if (t < 0.0) t = 0.0;
if (t > 1.0) t = 1.0;
}
return {p1.x + t * (p2.x - p1.x), p1.y + t * (p2.y - p1.y)};
}
void compute_heights(std::vector<double> &heights, int n, double step, int threads) {
auto worker = [&](int row_start, int row_end) {
for (int i = row_start; i < row_end; ++i) {
const double x = i * step;
const int base = i * n;
for (int j = 0; j < n; ++j) {
const double y = j * step;
heights[base + j] = height(x, y);
}
}
};
if (threads <= 1) {
worker(0, n);
return;
}
int use_threads = std::min(threads, n);
int chunk = (n + use_threads - 1) / use_threads;
std::vector<std::thread> pool;
pool.reserve(use_threads);
for (int t = 0; t < use_threads; ++t) {
int start = t * chunk;
int end = std::min(n, start + chunk);
if (start >= end) break;
pool.emplace_back(worker, start, end);
}
for (auto &th : pool) th.join();
}
double minimax_height(const std::vector<double> &heights, int n, int start_idx, int goal_idx) {
const double inf = std::numeric_limits<double>::infinity();
std::vector<double> best(heights.size(), inf);
using Node = std::pair<double, int>;
std::priority_queue<Node, std::vector<Node>, std::greater<Node>> pq;
best[start_idx] = heights[start_idx];
pq.emplace(best[start_idx], start_idx);
const int dirs[4][2] = {{1, 0}, {-1, 0}, {0, 1}, {0, -1}};
while (!pq.empty()) {
auto [cost, idx] = pq.top();
pq.pop();
if (cost != best[idx]) continue;
if (idx == goal_idx) return cost;
int i = idx / n;
int j = idx - i * n;
for (auto &d : dirs) {
int ni = i + d[0];
int nj = j + d[1];
if (ni < 0 || ni >= n || nj < 0 || nj >= n) continue;
int nidx = ni * n + nj;
double next = std::max(cost, heights[nidx]);
if (next < best[nidx]) {
best[nidx] = next;
pq.emplace(next, nidx);
}
}
}
return inf;
}
bool visible(const Point &a, const Point &b, double f, double step) {
const double dx = b.x - a.x;
const double dy = b.y - a.y;
const double dist_ab = std::hypot(dx, dy);
if (dist_ab == 0.0) return true;
const int samples = std::max(1, static_cast<int>(dist_ab / step));
for (int s = 1; s < samples; ++s) {
const double t = static_cast<double>(s) / samples;
const double x = a.x + dx * t;
const double y = a.y + dy * t;
if (height(x, y) > f + 1e-10) return false;
}
return true;
}
struct EdgePair {
int a = -1;
int b = -1;
};
} // namespace
int main(int argc, char **argv) {
constexpr int kMaxCoord = 1600;
double grid_step = 1.0;
double vis_step = 0.5;
int threads = static_cast<int>(std::thread::hardware_concurrency());
if (threads <= 0) threads = 1;
bool verbose = false;
bool validate = false;
for (int i = 1; i < argc; ++i) {
std::string arg = argv[i];
if (arg == "--threads" && i + 1 < argc) {
threads = std::max(1, std::atoi(argv[++i]));
} else if (arg == "--grid-step" && i + 1 < argc) {
grid_step = std::max(0.25, std::atof(argv[++i]));
} else if (arg == "--vis-step" && i + 1 < argc) {
vis_step = std::max(0.05, std::atof(argv[++i]));
} else if (arg == "--verbose") {
verbose = true;
} else if (arg == "--validate") {
validate = true;
}
}
const int n = static_cast<int>(kMaxCoord / grid_step) + 1;
std::vector<double> heights(static_cast<size_t>(n) * n, 0.0);
compute_heights(heights, n, grid_step, threads);
const Point A{200.0, 200.0};
const Point B{1400.0, 1400.0};
const int start_idx = static_cast<int>(A.x / grid_step) * n + static_cast<int>(A.y / grid_step);
const int goal_idx = static_cast<int>(B.x / grid_step) * n + static_cast<int>(B.y / grid_step);
const double f_min = minimax_height(heights, n, start_idx, goal_idx);
if (!std::isfinite(f_min)) {
std::cerr << "Failed to find a connecting elevation.\n";
return 1;
}
if (f_min + 1e-9 < height(A.x, A.y) || f_min + 1e-9 < height(B.x, B.y)) {
std::cerr << "Validation failed: f_min below endpoint height.\n";
return 1;
}
if (verbose) {
std::cerr << std::fixed << std::setprecision(6) << "f_min=" << f_min << "\n";
}
std::vector<uint8_t> inside(static_cast<size_t>(n) * n, 0);
for (size_t idx = 0; idx < heights.size(); ++idx) {
inside[idx] = heights[idx] > f_min ? 1 : 0;
}
static const int kCaseCount[16] = {
0, 1, 1, 1, 1, 2, 1, 1,
1, 1, 2, 1, 1, 1, 1, 0
};
static const EdgePair kCaseEdges[16][2] = {
{{-1, -1}, {-1, -1}}, // 0
{{3, 0}, {-1, -1}}, // 1
{{0, 1}, {-1, -1}}, // 2
{{3, 1}, {-1, -1}}, // 3
{{1, 2}, {-1, -1}}, // 4
{{3, 2}, {0, 1}}, // 5
{{0, 2}, {-1, -1}}, // 6
{{3, 2}, {-1, -1}}, // 7
{{2, 3}, {-1, -1}}, // 8
{{0, 2}, {-1, -1}}, // 9
{{0, 1}, {2, 3}}, // 10
{{1, 2}, {-1, -1}}, // 11
{{3, 1}, {-1, -1}}, // 12
{{0, 1}, {-1, -1}}, // 13
{{3, 0}, {-1, -1}}, // 14
{{-1, -1}, {-1, -1}} // 15
};
std::vector<std::pair<Point, Point>> segments;
segments.reserve(9000);
for (int i = 0; i < n - 1; ++i) {
const double x = i * grid_step;
for (int j = 0; j < n - 1; ++j) {
const double y = j * grid_step;
const int idx_ll = i * n + j;
const int idx_lr = (i + 1) * n + j;
const int idx_ur = (i + 1) * n + (j + 1);
const int idx_ul = i * n + (j + 1);
const bool c0 = inside[idx_ll];
const bool c1 = inside[idx_lr];
const bool c2 = inside[idx_ur];
const bool c3 = inside[idx_ul];
const int mask = (c0 ? 1 : 0) | (c1 ? 2 : 0) | (c2 ? 4 : 0) | (c3 ? 8 : 0);
if (mask == 0 || mask == 15) continue;
const Point ll{x, y};
const Point lr{x + grid_step, y};
const Point ur{x + grid_step, y + grid_step};
const Point ul{x, y + grid_step};
const double v_ll = heights[idx_ll];
const double v_lr = heights[idx_lr];
const double v_ur = heights[idx_ur];
const double v_ul = heights[idx_ul];
Point edges[4];
edges[0] = interp(ll, lr, v_ll, v_lr, f_min);
edges[1] = interp(lr, ur, v_lr, v_ur, f_min);
edges[2] = interp(ul, ur, v_ul, v_ur, f_min);
edges[3] = interp(ll, ul, v_ll, v_ul, f_min);
for (int s = 0; s < kCaseCount[mask]; ++s) {
int e0 = kCaseEdges[mask][s].a;
int e1 = kCaseEdges[mask][s].b;
if (e0 < 0 || e1 < 0) continue;
segments.emplace_back(edges[e0], edges[e1]);
}
}
}
if (segments.empty()) {
std::cerr << "Validation failed: no contour segments found.\n";
return 1;
}
const double key_scale = 1e6;
std::unordered_map<Key, std::vector<Key>, KeyHash> adj;
adj.reserve(segments.size() * 2);
for (const auto &seg : segments) {
Key k1 = quantize(seg.first, key_scale);
Key k2 = quantize(seg.second, key_scale);
adj[k1].push_back(k2);
adj[k2].push_back(k1);
}
for (const auto &kv : adj) {
if (kv.second.size() != 2) {
std::cerr << "Validation failed: contour node degree != 2.\n";
return 1;
}
}
std::unordered_set<Key, KeyHash> visited;
std::vector<Key> contour_keys;
contour_keys.reserve(adj.size());
double best_perimeter = -1.0;
for (const auto &kv : adj) {
const Key start_key = kv.first;
if (visited.count(start_key)) continue;
std::vector<Key> loop;
loop.reserve(1024);
Key prev{std::numeric_limits<int64_t>::min(), std::numeric_limits<int64_t>::min()};
Key cur = start_key;
while (true) {
if (visited.count(cur)) {
if (cur == start_key) break;
std::cerr << "Validation failed: contour loop revisited unexpectedly.\n";
return 1;
}
visited.insert(cur);
loop.push_back(cur);
const auto &nbrs = adj[cur];
Key next = nbrs[0];
if (prev == next) next = nbrs[1];
prev = cur;
cur = next;
if (loop.size() > adj.size() + 1) {
std::cerr << "Validation failed: contour walk overflow.\n";
return 1;
}
}
double per = 0.0;
for (size_t i = 0; i < loop.size(); ++i) {
Point a = from_key(loop[i], key_scale);
Point b = from_key(loop[(i + 1) % loop.size()], key_scale);
per += dist(a, b);
}
if (per > best_perimeter) {
best_perimeter = per;
contour_keys = std::move(loop);
}
}
if (contour_keys.empty()) {
std::cerr << "Validation failed: contour loop not found.\n";
return 1;
}
if (validate && visited.size() != adj.size()) {
std::cerr << "Validation warning: multiple contour loops detected; using largest.\n";
}
std::vector<Point> contour;
contour.reserve(contour_keys.size());
for (const auto &k : contour_keys) {
contour.push_back(from_key(k, key_scale));
}
const int m = static_cast<int>(contour.size());
std::vector<double> prefix(static_cast<size_t>(m) + 1, 0.0);
for (int i = 0; i < m; ++i) {
prefix[i + 1] = prefix[i] + dist(contour[i], contour[(i + 1) % m]);
}
const double perimeter = prefix[m];
std::vector<int> visA;
std::vector<int> visB;
std::vector<double> distA(m, 0.0);
std::vector<double> distB(m, 0.0);
visA.reserve(m);
visB.reserve(m);
for (int i = 0; i < m; ++i) {
distA[i] = dist(A, contour[i]);
distB[i] = dist(B, contour[i]);
if (visible(A, contour[i], f_min, vis_step)) visA.push_back(i);
if (visible(B, contour[i], f_min, vis_step)) visB.push_back(i);
}
if (visA.empty() || visB.empty()) {
std::cerr << "Validation failed: no visible contour points from A or B.\n";
return 1;
}
double best = std::numeric_limits<double>::infinity();
if (visible(A, B, f_min, vis_step)) best = dist(A, B);
for (int i : visA) {
for (int j : visB) {
double arc = std::abs(prefix[j] - prefix[i]);
arc = std::min(arc, perimeter - arc);
double total = distA[i] + distB[j] + arc;
if (total < best) best = total;
}
}
if (!std::isfinite(best)) {
std::cerr << "Validation failed: shortest path not found.\n";
return 1;
}
if (validate) {
const double direct = dist(A, B);
if (best + 1e-7 < direct) {
std::cerr << "Validation failed: shortest path shorter than straight line.\n";
return 1;
}
if (contour.size() < 1000) {
std::cerr << "Validation warning: contour resolution seems too low.\n";
}
}
std::cout << std::fixed << std::setprecision(3) << best << "\n";
return 0;
}
Python
import sys
import math
import heapq
from concurrent.futures import ThreadPoolExecutor
def compute_height(x, y):
xx = x * x
yy = y * y
xy = x * y
base = 5000.0 - 0.005 * (xx + yy + xy) + 12.5 * (x + y)
expo = -abs(0.000001 * (xx + yy) - 0.0015 * (x + y) + 0.7)
return base * math.exp(expo)
def dist(a, b):
return math.hypot(a[0] - b[0], a[1] - b[1])
def quantize(p, scale):
return (round(p[0] * scale), round(p[1] * scale))
def from_key(k, scale):
return (k[0] / scale, k[1] / scale)
def interp(p1, p2, v1, v2, f):
denom = v2 - v1
if abs(denom) < 1e-14:
t = 0.5
else:
t = (f - v1) / denom
if t < 0.0: t = 0.0
if t > 1.0: t = 1.0
return (p1[0] + t * (p2[0] - p1[0]), p1[1] + t * (p2[1] - p1[1]))
def worker_heights(row_start, row_end, n, step):
local_heights = {}
for i in range(row_start, row_end):
x = i * step
base = i * n
for j in range(n):
y = j * step
local_heights[base + j] = compute_height(x, y)
return local_heights
def minimax_height(heights, n, start_idx, goal_idx):
inf = float('inf')
best = [inf] * len(heights)
pq = []
best[start_idx] = heights[start_idx]
heapq.heappush(pq, (best[start_idx], start_idx))
dirs = [(1, 0), (-1, 0), (0, 1), (0, -1)]
while pq:
cost, idx = heapq.heappop(pq)
if cost != best[idx]:
continue
if idx == goal_idx:
return cost
i = idx // n
j = idx % n
for di, dj in dirs:
ni = i + di
nj = j + dj
if 0 <= ni < n and 0 <= nj < n:
nidx = ni * n + nj
nxt = max(cost, heights[nidx])
if nxt < best[nidx]:
best[nidx] = nxt
heapq.heappush(pq, (nxt, nidx))
return inf
def visible(a, b, f, step):
dx = b[0] - a[0]
dy = b[1] - a[1]
dist_ab = math.hypot(dx, dy)
if dist_ab == 0.0:
return True
samples = max(1, int(dist_ab / step))
for s in range(1, samples):
t = s / samples
x = a[0] + dx * t
y = a[1] + dy * t
if compute_height(x, y) > f + 1e-10:
return False
return True
def solve():
kMaxCoord = 1600
grid_step = 1.0
vis_step = 0.5
threads = 4 # Adjust based on logic, Python threading is GIL limited, but computation is isolated if we use ProcessPool, however ThreadPool is easier and math is somewhat fast here. Wait, actually we can just sequentially do it or use ProcessPool.
import multiprocessing
threads = multiprocessing.cpu_count() or 1
n = int(kMaxCoord / grid_step) + 1
heights = [0.0] * (n * n)
# Compute heights
import os
if threads > 1 and os.name == 'posix':
with multiprocessing.Pool(threads) as pool:
chunk = (n + threads - 1) // threads
args = []
for t in range(threads):
start = t * chunk
end = min(n, start + chunk)
if start < end:
args.append((start, end, n, grid_step))
results = pool.starmap(worker_heights, args)
for res_dict in results:
for idx, h in res_dict.items():
heights[idx] = h
else:
for idx, h in worker_heights(0, n, n, grid_step).items():
heights[idx] = h
A = (200.0, 200.0)
B = (1400.0, 1400.0)
start_idx = int(A[0] / grid_step) * n + int(A[1] / grid_step)
goal_idx = int(B[0] / grid_step) * n + int(B[1] / grid_step)
f_min = minimax_height(heights, n, start_idx, goal_idx)
inside = [0] * (n * n)
for idx in range(len(heights)):
if heights[idx] > f_min:
inside[idx] = 1
kCaseCount = [0, 1, 1, 1, 1, 2, 1, 1, 1, 1, 2, 1, 1, 1, 1, 0]
kCaseEdges = [
[(-1, -1), (-1, -1)],
[(3, 0), (-1, -1)],
[(0, 1), (-1, -1)],
[(3, 1), (-1, -1)],
[(1, 2), (-1, -1)],
[(3, 2), (0, 1)],
[(0, 2), (-1, -1)],
[(3, 2), (-1, -1)],
[(2, 3), (-1, -1)],
[(0, 2), (-1, -1)],
[(0, 1), (2, 3)],
[(1, 2), (-1, -1)],
[(3, 1), (-1, -1)],
[(0, 1), (-1, -1)],
[(3, 0), (-1, -1)],
[(-1, -1), (-1, -1)]
]
segments = []
for i in range(n - 1):
x = i * grid_step
for j in range(n - 1):
y = j * grid_step
idx_ll = i * n + j
idx_lr = (i + 1) * n + j
idx_ur = (i + 1) * n + (j + 1)
idx_ul = i * n + (j + 1)
mask = (inside[idx_ll] | (inside[idx_lr] << 1) | (inside[idx_ur] << 2) | (inside[idx_ul] << 3))
if mask == 0 or mask == 15:
continue
ll = (x, y)
lr = (x + grid_step, y)
ur = (x + grid_step, y + grid_step)
ul = (x, y + grid_step)
v_ll = heights[idx_ll]
v_lr = heights[idx_lr]
v_ur = heights[idx_ur]
v_ul = heights[idx_ul]
edges = [
interp(ll, lr, v_ll, v_lr, f_min),
interp(lr, ur, v_lr, v_ur, f_min),
interp(ul, ur, v_ul, v_ur, f_min),
interp(ll, ul, v_ll, v_ul, f_min)
]
for s in range(kCaseCount[mask]):
e0, e1 = kCaseEdges[mask][s]
if e0 < 0 or e1 < 0:
continue
segments.append((edges[e0], edges[e1]))
key_scale = 1e6
adj = {}
for seg in segments:
k1 = quantize(seg[0], key_scale)
k2 = quantize(seg[1], key_scale)
if k1 not in adj: adj[k1] = []
if k2 not in adj: adj[k2] = []
adj[k1].append(k2)
adj[k2].append(k1)
visited = set()
contour_keys = []
best_perimeter = -1.0
for start_key in adj:
if start_key in visited:
continue
loop = []
prev = (-float('inf'), -float('inf'))
cur = start_key
while True:
if cur in visited:
break
visited.add(cur)
loop.append(cur)
nbrs = adj[cur]
nxt = nbrs[0]
if prev == nxt:
if len(nbrs) > 1:
nxt = nbrs[1]
prev = cur
cur = nxt
per = 0.0
for i in range(len(loop)):
a = from_key(loop[i], key_scale)
b = from_key(loop[(i + 1) % len(loop)], key_scale)
per += dist(a, b)
if per > best_perimeter:
best_perimeter = per
contour_keys = loop
contour = [from_key(k, key_scale) for k in contour_keys]
m = len(contour)
prefix = [0.0] * (m + 1)
for i in range(m):
prefix[i + 1] = prefix[i] + dist(contour[i], contour[(i + 1) % m])
perimeter = prefix[m]
visA = []
visB = []
distA = [0.0] * m
distB = [0.0] * m
for i in range(m):
distA[i] = dist(A, contour[i])
distB[i] = dist(B, contour[i])
if visible(A, contour[i], f_min, vis_step): visA.append(i)
if visible(B, contour[i], f_min, vis_step): visB.append(i)
best = float('inf')
if visible(A, B, f_min, vis_step):
best = dist(A, B)
for i in visA:
for j in visB:
arc = abs(prefix[j] - prefix[i])
arc = min(arc, perimeter - arc)
total = distA[i] + distB[j] + arc
if total < best:
best = total
return f"{best:.3f}"
if __name__ == '__main__':
print(solve())
Java
import java.util.*;
import java.util.concurrent.*;
public class Euler262 {
static class Point {
double x, y;
Point(double x, double y) {
this.x = x;
this.y = y;
}
}
static double height(double x, double y) {
double xx = x * x;
double yy = y * y;
double xy = x * y;
double base = 5000.0 - 0.005 * (xx + yy + xy) + 12.5 * (x + y);
double expo = -Math.abs(0.000001 * (xx + yy) - 0.0015 * (x + y) + 0.7);
return base * Math.exp(expo);
}
static double dist(Point a, Point b) {
return Math.hypot(a.x - b.x, a.y - b.y);
}
static class Key {
long x, y;
Key(long x, long y) {
this.x = x;
this.y = y;
}
@Override
public boolean equals(Object o) {
if (this == o)
return true;
if (!(o instanceof Key))
return false;
Key other = (Key) o;
return x == other.x && y == other.y;
}
@Override
public int hashCode() {
long hashX = x;
return (int) ((hashX * 0x9e3779b97f4a7c15L) ^ (y + (hashX << 6) + (hashX >> 2)));
}
}
static Key quantize(Point p, double scale) {
return new Key(Math.round(p.x * scale), Math.round(p.y * scale));
}
static Point fromKey(Key k, double scale) {
return new Point(k.x / scale, k.y / scale);
}
static Point interp(Point p1, Point p2, double v1, double v2, double f) {
double t;
double denom = v2 - v1;
if (Math.abs(denom) < 1e-14) {
t = 0.5;
} else {
t = (f - v1) / denom;
if (t < 0.0)
t = 0.0;
if (t > 1.0)
t = 1.0;
}
return new Point(p1.x + t * (p2.x - p1.x), p1.y + t * (p2.y - p1.y));
}
static class Node implements Comparable<Node> {
double cost;
int idx;
Node(double cost, int idx) {
this.cost = cost;
this.idx = idx;
}
public int compareTo(Node o) {
return Double.compare(this.cost, o.cost);
}
}
public static String solve() {
int kMaxCoord = 1600;
double gridStep = 1.0;
double visStep = 0.5;
int threads = Math.max(1, Runtime.getRuntime().availableProcessors());
int n = (int) (kMaxCoord / gridStep) + 1;
double[] heights = new double[n * n];
int useThreads = Math.min(threads, n);
int chunk = (n + useThreads - 1) / useThreads;
ExecutorService executor = Executors.newFixedThreadPool(useThreads);
List<Future<?>> futures = new ArrayList<>();
for (int t = 0; t < useThreads; ++t) {
final int start = t * chunk;
final int end = Math.min(n, start + chunk);
if (start >= end)
break;
futures.add(executor.submit(() -> {
for (int i = start; i < end; ++i) {
double x = i * gridStep;
int base = i * n;
for (int j = 0; j < n; ++j) {
double y = j * gridStep;
heights[base + j] = height(x, y);
}
}
}));
}
for (Future<?> f : futures) {
try {
f.get();
} catch (Exception e) {
}
}
executor.shutdown();
Point A = new Point(200.0, 200.0);
Point B = new Point(1400.0, 1400.0);
int startIdx = (int) (A.x / gridStep) * n + (int) (A.y / gridStep);
int goalIdx = (int) (B.x / gridStep) * n + (int) (B.y / gridStep);
double[] bestHeights = new double[n * n];
Arrays.fill(bestHeights, Double.POSITIVE_INFINITY);
PriorityQueue<Node> pq = new PriorityQueue<>();
bestHeights[startIdx] = heights[startIdx];
pq.offer(new Node(bestHeights[startIdx], startIdx));
int[][] dirs = { { 1, 0 }, { -1, 0 }, { 0, 1 }, { 0, -1 } };
double fMin = Double.POSITIVE_INFINITY;
while (!pq.isEmpty()) {
Node curr = pq.poll();
if (curr.cost != bestHeights[curr.idx])
continue;
if (curr.idx == goalIdx) {
fMin = curr.cost;
break;
}
int i = curr.idx / n;
int j = curr.idx % n;
for (int[] d : dirs) {
int ni = i + d[0];
int nj = j + d[1];
if (ni < 0 || ni >= n || nj < 0 || nj >= n)
continue;
int nidx = ni * n + nj;
double next = Math.max(curr.cost, heights[nidx]);
if (next < bestHeights[nidx]) {
bestHeights[nidx] = next;
pq.offer(new Node(next, nidx));
}
}
}
byte[] inside = new byte[n * n];
for (int idx = 0; idx < heights.length; ++idx) {
if (heights[idx] > fMin)
inside[idx] = 1;
}
int[] kCaseCount = { 0, 1, 1, 1, 1, 2, 1, 1, 1, 1, 2, 1, 1, 1, 1, 0 };
int[][][] kCaseEdges = {
{ { -1, -1 }, { -1, -1 } },
{ { 3, 0 }, { -1, -1 } },
{ { 0, 1 }, { -1, -1 } },
{ { 3, 1 }, { -1, -1 } },
{ { 1, 2 }, { -1, -1 } },
{ { 3, 2 }, { 0, 1 } },
{ { 0, 2 }, { -1, -1 } },
{ { 3, 2 }, { -1, -1 } },
{ { 2, 3 }, { -1, -1 } },
{ { 0, 2 }, { -1, -1 } },
{ { 0, 1 }, { 2, 3 } },
{ { 1, 2 }, { -1, -1 } },
{ { 3, 1 }, { -1, -1 } },
{ { 0, 1 }, { -1, -1 } },
{ { 3, 0 }, { -1, -1 } },
{ { -1, -1 }, { -1, -1 } }
};
List<Point[]> segments = new ArrayList<>();
for (int i = 0; i < n - 1; ++i) {
double x = i * gridStep;
for (int j = 0; j < n - 1; ++j) {
double y = j * gridStep;
int idxLL = i * n + j;
int idxLR = (i + 1) * n + j;
int idxUR = (i + 1) * n + (j + 1);
int idxUL = i * n + (j + 1);
int mask = (inside[idxLL] == 1 ? 1 : 0) |
(inside[idxLR] == 1 ? 2 : 0) |
(inside[idxUR] == 1 ? 4 : 0) |
(inside[idxUL] == 1 ? 8 : 0);
if (mask == 0 || mask == 15)
continue;
Point ll = new Point(x, y);
Point lr = new Point(x + gridStep, y);
Point ur = new Point(x + gridStep, y + gridStep);
Point ul = new Point(x, y + gridStep);
Point[] edges = {
interp(ll, lr, heights[idxLL], heights[idxLR], fMin),
interp(lr, ur, heights[idxLR], heights[idxUR], fMin),
interp(ul, ur, heights[idxUL], heights[idxUR], fMin),
interp(ll, ul, heights[idxLL], heights[idxUL], fMin)
};
for (int s = 0; s < kCaseCount[mask]; ++s) {
int e0 = kCaseEdges[mask][s][0];
int e1 = kCaseEdges[mask][s][1];
if (e0 < 0 || e1 < 0)
continue;
segments.add(new Point[] { edges[e0], edges[e1] });
}
}
}
double keyScale = 1e6;
Map<Key, List<Key>> adj = new HashMap<>();
for (Point[] seg : segments) {
Key k1 = quantize(seg[0], keyScale);
Key k2 = quantize(seg[1], keyScale);
adj.computeIfAbsent(k1, k -> new ArrayList<>()).add(k2);
adj.computeIfAbsent(k2, k -> new ArrayList<>()).add(k1);
}
Set<Key> visited = new HashSet<>();
List<Key> contourKeys = new ArrayList<>();
double bestPerimeter = -1.0;
Key prevDummy = new Key(Long.MIN_VALUE, Long.MIN_VALUE);
for (Key startKey : adj.keySet()) {
if (visited.contains(startKey))
continue;
List<Key> loop = new ArrayList<>();
Key prev = prevDummy;
Key cur = startKey;
while (true) {
if (visited.contains(cur))
break;
visited.add(cur);
loop.add(cur);
List<Key> nbrs = adj.get(cur);
Key nxt = nbrs.get(0);
if (prev.equals(nxt) && nbrs.size() > 1) {
nxt = nbrs.get(1);
}
prev = cur;
cur = nxt;
}
double per = 0.0;
for (int i = 0; i < loop.size(); ++i) {
Point a = fromKey(loop.get(i), keyScale);
Point b = fromKey(loop.get((i + 1) % loop.size()), keyScale);
per += dist(a, b);
}
if (per > bestPerimeter) {
bestPerimeter = per;
contourKeys = loop;
}
}
List<Point> contour = new ArrayList<>();
for (Key k : contourKeys) {
contour.add(fromKey(k, keyScale));
}
int m = contour.size();
double[] prefix = new double[m + 1];
for (int i = 0; i < m; ++i) {
prefix[i + 1] = prefix[i] + dist(contour.get(i), contour.get((i + 1) % m));
}
double perimeter = prefix[m];
List<Integer> visA = new ArrayList<>();
List<Integer> visB = new ArrayList<>();
double[] distA = new double[m];
double[] distB = new double[m];
for (int i = 0; i < m; ++i) {
distA[i] = dist(A, contour.get(i));
distB[i] = dist(B, contour.get(i));
if (visible(A, contour.get(i), fMin, visStep))
visA.add(i);
if (visible(B, contour.get(i), fMin, visStep))
visB.add(i);
}
double bestDist = Double.POSITIVE_INFINITY;
if (visible(A, B, fMin, visStep)) {
bestDist = dist(A, B);
}
for (int i : visA) {
for (int j : visB) {
double arc = Math.abs(prefix[j] - prefix[i]);
arc = Math.min(arc, perimeter - arc);
double total = distA[i] + distB[j] + arc;
if (total < bestDist) {
bestDist = total;
}
}
}
return String.format(Locale.US, "%.3f", bestDist);
}
static boolean visible(Point a, Point b, double f, double step) {
double dx = b.x - a.x;
double dy = b.y - a.y;
double distAb = Math.hypot(dx, dy);
if (distAb == 0.0)
return true;
int samples = Math.max(1, (int) (distAb / step));
for (int s = 1; s < samples; ++s) {
double t = (double) s / samples;
if (height(a.x + dx * t, a.y + dy * t) > f + 1e-10)
return false;
}
return true;
}
public static void main(String[] args) {
System.out.println(solve());
}
}