Problem 262: Mountain Range

Project 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 height of a mountainous landscape above the point $(x,y)$ is $h=P\,e^{-|Q|}$, where

$$\begin{aligned}P(x,y)&=5000-\frac{x^2+y^2+xy}{200}\\&\quad+\frac{25(x+y)}{2},\end{aligned}$$

$$\begin{aligned}Q(x,y)&=\frac{x^2+y^2}{10^6}-\frac{3(x+y)}{2000}\\&\quad+\frac{7}{10}.\end{aligned}$$

A mosquito at $A=(200,200)$ rises vertically to an elevation $f$ and then flies at this constant elevation, without leaving the square $0\le x,y\le1600$, until it is directly above $B=(1400,1400)$; it can only pass above points where $h\le f$. With $f_{\min}$ the smallest elevation that allows such a trip, we must find the length of the shortest path between the points $A'$ and $B'$ above $A$ and $B$ at the elevation $f_{\min}$, rounded to three decimal places. ^[1] The flight is horizontal, so this is the length of a plane curve from $A$ to $B$ in the square along which $h\le f_{\min}$. The statement gives no intermediate values.

The whole solution in one sentence: the mountains form a ring whose crest is too high to cross, so the mosquito must fly round it through a pass at the highest point of the left edge of the square (or of the bottom edge, in mirror image), whose height is $f_{\min}$, and at this level the shortest flight is a straight line, an arc of the contour line $h=f_{\min}$ and another straight line, which the program approximates on a grid of heights.

Mathematical Approach

We first describe the landscape exactly, find the pass and prove that its height is $f_{\min}$, determine the shape of the shortest flight, and then follow the grid methods of the program.

1. A ring of mountains

Completing the squares in $Q$ gives

$$\begin{aligned}10^6\,Q&=(x-750)^2+(y-750)^2\\&\quad-425000.\end{aligned}$$

So $Q=0$ exactly on the circle $K$ with centre $M=(750,750)$ and radius $\rho=\sqrt{425000}\approx651.92$, $Q\lt0$ inside it and $Q\gt0$ outside. Expanding $P$ around $M$ with $u=x-750$, $v=y-750$ and $w=u+v=x+y-1500$ gives

$$P=15312.5-\frac{u^2+v^2+uv}{200}+\frac{5w}{4}.$$

On $K$ we have $u^2+v^2=\rho^2$ and $uv=(w^2-\rho^2)/2$, so there the height $h=P$ is

$$14250-\frac{w^2}{400}+\frac{5w}{4},$$

where $|w|\le\rho\sqrt2\approx921.95$ because $(u+v)^2\le2(u^2+v^2)$. This parabola in $w$ is highest, $14406.25$, at $w=250$, and lowest at $w=-\rho\sqrt2$, the point $(289.02,289.02)$ of $K$ nearest to $A$, where it equals $12125-1.25\rho\sqrt2$${}\approx10972.557$. So every point of $K$ is at least $10972.55$ high. It is a sharp crest: with $n$ the outward normal, the slope of $h$ is $\partial P/\partial n+P|\nabla Q|$ just inside $K$ and $\partial P/\partial n-P|\nabla Q|$ just outside, where $|\nabla Q|=2\rho/10^6$, $P|\nabla Q|\ge14.3$ and $\partial P/\partial n$ lies between $-11.55$ and $-3.14$ on $K$, so $h$ rises up to $K$ and falls after it. The centre lies in a basin, $h(M)=15312.5\,e^{-0.425}$${}\approx10010.85$. Both $A$ and $B$ lie outside $K$, at the distances $777.82$ and $919.24$ from $M$, with $h(A)\approx7851.54$ and $h(B)\approx6964.70$. The segment $AB$ crosses $K$ twice and is highest at the second crossing $(1210.98,1210.98)$, where $h\approx13277.44$.

2. The lowest elevation is a minimax value

A trip at the elevation $f$ exists if and only if some path $\gamma$ from $A$ to $B$ in the square has its highest point at most $f$. Hence

$$f_{\min}=\min_{\gamma}\,\max_{p\in\gamma}h(p)$$

over all such paths. This is a minimax, or bottleneck, problem: only the highest point of $\gamma$ matters, not its length. ^[2] Picture a flood: as the water level $f$ rises, the lakes around $A$ and $B$ grow, and $f_{\min}$ is the level at which they first touch. In a smooth landscape this usually happens at a saddle point, the mountain pass; ^[4] here it happens on the edge of the square.

3. The pass on the edge of the square

On the left edge $x=0$ we have $Q\ge0.1375\gt0$, so the absolute value can be dropped, and differentiating gives

$$\frac{d}{dy}\,h(0,y)=\frac{e^{-Q(0,y)}}{10^8}\,c(y),$$

$$\begin{aligned}c(y)&=y^3-3250y^2-125000y\\&\quad+2\cdot10^9.\end{aligned}$$

For $0\le y\le1600$ the derivative $c'(y)=y(3y-6500)-125000$ is negative, and $c(0)\gt0\gt c(1600)$, so $c$ has exactly one root there, $y_0\approx895.4834$. The height increases along the edge up to $y_0$ and decreases after it, so the top of the edge is $S=(0,y_0)$ with $h(S)\approx10396.4622$; since $h(x,y)=h(y,x)$, the point $S'=(y_0,0)$ on the bottom edge is equally high. We show that $f_{\min}=h(S)$, using two facts: $P$ is concave, so on the square it is at least its smallest corner value $5000$; and where $Q\ge0$, $h_x=e^{-Q}(P_x-PQ_x)$ with $P_x=12.5-(2x+y)/200$ and $Q_x=(x-750)/500000$, and similarly for $y$.

A route at this height. Fly from $A$ straight to $(0,200)$, up the left edge to $(0,1600)$, along the top edge to $(1400,1600)$ and down to $B$. On the first leg $P_x\gt0\gt Q_x$, so $h$ grows towards $A$; on the left edge $h\le h(S)$; on the top edge $P\le13212.5$ and $Q\ge0.2975$, so $h\le13212.5\,e^{-0.2975}\lt9813$; on the last leg $P_y\lt0\lt Q_y$, so $h\le h(B)$. Hence $f_{\min}\le h(S)$.
No lower route. Let $C_1=(a,y_0)$, $a\approx114.52$, be the left point of $K$ at the height $y_0$, $L_1$ the segment from $S$ to $C_1$, and $L_2$ its mirror image from $S'$ to $C_2=(y_0,a)$. On $L_1$ we have $Q\ge0$, $Q_x\lt0$ and, as $x\lt115$ and $y_0\lt896$, $P_x\gt6.87$, so $h$ increases strictly from $S$ to $C_1$, and $h\ge h(S)$ on $L_1$ and $L_2$, with equality only at $S$ and $S'$. The edges from $S$ to the corner $(0,0)$ and on to $S'$, the segment $L_2$, the arc of $K$ from $C_2$ to $C_1$ through $(289.02,289.02)$ and $L_1$ form a simple closed curve around $A$ that leaves $B$ outside. By the Jordan curve theorem every path from $A$ to $B$ meets this curve, ^[5] so it has a last point on it. A path inside the square cannot leave the enclosed region across the two edge pieces, so that point lies on $L_1$, $L_2$ (with $S$ and $S'$) or $K$, where $h\ge h(S)$. Hence $f_{\min}\ge h(S)$.

So $f_{\min}=h(S)$. The pass is not a saddle point of $h$: at $S$ the height still increases into the square, $\partial h/\partial x\approx22.44$. It is made by the boundary of the area; outside the square the mosquito could fly much lower, for example $h(-200,895)\approx6303.49$. At the level $f_{\min}$ the strict increase along $L_1$ and $L_2$ forces every flight through $S$ or $S'$.

4. The same level on a grid

The program samples $h$ at the $1601^2=2\,563\,201$ points with integer coordinates and looks for a chain of grid points from $A$ to $B$, each a horizontal or vertical neighbour of the previous one, whose highest point is as low as possible: Dijkstra's algorithm with the length of a chain replaced by its highest point. ^[2] Every point $v$ carries a label $b(v)$, the lowest highest point of the chains from $A$ to $v$ found so far, starting with $b(A)=h(A)$; the algorithm repeatedly finishes the unfinished point $v$ with the smallest label and updates each neighbour $u$ by

$$b(u)\leftarrow\min\bigl(b(u),\,\max(b(v),h(u))\bigr).$$

Dijkstra's argument works because extending a chain never lowers its highest point and $\max(b,h(u))$ grows with $b$. ^[3] Suppose all finished points, starting with $A$, carry optimal labels, and $v$ is chosen next. Any chain from $A$ to $v$ has a first unfinished point $u$; the update from its finished predecessor made $b(u)$ at most the highest point of the chain up to $u$, and $b(v)\le b(u)$, so no chain beats $b(v)$, which itself comes from a chain. A tiny example, with $A$ at the bottom left (height $1$) and $B$ at the top right (height $3$):

$$\begin{matrix}4&9&3\\5&8&6\\1&7&2\end{matrix}$$

The points are finished in the order of the heights $1,5,4$, with the labels $1,5,5$, and then $7,2,6,3$, all with the label $7$. Every chain passes the middle column, with the heights $7$, $8$ and $9$, so it climbs to at least $7$, and the chain along the bottom row and up the right column attains $7$. On the real grid the C++ and Python searches take $454\,302$ points out of the queue, $B$ last, and return $h(0,895)\approx10396.4587$, $0.0035$ below $f_{\min}$: the grid points $(0,895)$ and $(0,896)$ lie on both sides of the true top at $y_0\approx895.48$, and both are lower. At this level the continuous problem has no flight, since the obstacle reaches over the edge for $895\lt y\lt895.97$, but no grid point there is above the level.

5. The obstacle is convex

At the level $f_{\min}$ the obstacle, the region above the level, is a ring around $K$. A flight can never enter the disk inside $K$, because $K$ is higher than $f_{\min}$, so we may fill the basin. Let $G=P\,e^{-Q}$; it equals $h$ outside $K$ and exceeds it inside. $P$ is concave and positive on the square and $Q$ is convex, so $P$ and $e^{-Q}$ are log-concave, hence so is their product $G$, and its superlevel sets are convex. ^[10] The set $D$ of points of the square with $G\ge f_{\min}$ is therefore convex. It contains $K$, where $G=P\gt f_{\min}$, hence the whole disk, and outside the disk it is exactly the obstacle. So the outer contour of the obstacle, the boundary of $D$, is a convex curve, and by the bounds of Section 3 it touches the edges only at $S$ and $S'$. The reflection $(x,y)\mapsto(y,x)$ fixes $h$, $A$ and $B$ and exchanges $S$ and $S'$, so the shortest flights through $S$ and $S'$ are mirror images of equal length; we follow the one through $S$.

6. Shortest paths round a convex obstacle

Picture a string from $A$ to $B$ pulled tight round the obstacle. Where it is free it is straight, since a curved piece could be replaced by a shorter chord; where it touches the obstacle it follows the contour; and it meets and leaves the contour along tangents, since a corner at the junction could be cut off. For polygonal obstacles this is the rule that a shortest path bends only at their corners. ^[7] A worked example: round the disk of radius $3$ about the origin from $(-5,0)$ to $(5,0)$. The tangent from $(-5,0)$ has the length $\sqrt{5^2-3^2}=4$, and the radius to its upper tangent point makes the angle $\alpha=\arccos(3/5)\approx0.9273$ with the negative $x$-axis. Between the two tangent points the string follows an arc of the angle $\pi-2\alpha\approx1.2870$, of length $3(\pi-2\alpha)\approx3.8610$, so the shortest path has the length $8+3(\pi-2\alpha)\approx11.8610$.

7. The shortest flight

The flight through $S$ is a shortest path from $A$ to $S$ followed by one from $S$ to $B$, both round the convex set $D$. By Section 6 it goes straight from $A$ to a tangent point $T_A$, follows the contour upwards through $S$ and round the upper left part of the mountain to a tangent point $T_B$, and goes straight to $B$. Since $\nabla h$ is perpendicular to the contour, a tangent point $T$ seen from $X$ solves $h(T)=f_{\min}$ and $(T-X)\cdot\nabla h(T)=0$; this gives $T_A\approx(48.3,624.7)$ and $T_B\approx(873.0,1536.0)$. With $s$ for the length of an arc of the contour, the flight has the length

$$|AT_A|+s(T_A,T_B)+|T_BB|.$$

The straight pieces stay in the square because the square is convex.

8. Tracing the contour with marching squares

A grid point is inside the obstacle when its height exceeds the level $f$ of the search; each unit square of the grid, a cell, is treated on its own. Its corners give a $4$-bit number from $0$ to $15$, with the bits $1$, $2$, $4$, $8$ for the lower left, lower right, upper right and upper left corner. Cells with the number $0$ or $15$ are skipped; otherwise a table of the $16$ cases names the two sides on which the contour enters and leaves the cell, and the crossing point on a side with the corner heights $v_1,v_2$ lies at the fraction

$$t=\frac{f-v_1}{v_2-v_1}$$

of the way from the first corner, by linear interpolation. ^[6] For example, with $f=10396$ and the corner heights $10390$, $10400$, $10405$, $10395$ in this order, the right corners are inside, the number is $2+4=6$, the table gives the bottom and top sides, and with $t=0.6$ and $t=0.1$ the contour crosses the cell from $(0.6,0)$ to $(0.1,1)$. The numbers $5$ and $10$, with two opposite corners inside, need two segments, but no such cell occurs here: $8\,444$ cells are crossed, each by one segment.

9. From pieces to closed curves

Every crossing point is computed in both cells that share its side. To merge the copies safely, the program multiplies the coordinates by $10^6$, rounds them to integers and records the neighbours of every rounded point; all $8\,444$ points have exactly two. Walking from neighbour to neighbour traces closed curves, here two: the outer contour with $6\,244$ points and a length of about $4895.10$, and the inner one with $2\,200$ points and about $1719.81$; the program keeps the longer. For its points $P_0,\dots,P_{m-1}$ and $P_m=P_0$ it stores the sums

$$S_0=0,\qquad S_{k+1}=S_k+|P_kP_{k+1}|,$$

so the total length is $L=S_m$ and the two arcs between $P_i$ and $P_j$ have the lengths $s_{ij}=|S_j-S_i|$ and $L-s_{ij}$. The traced points lie within $2\cdot10^{-4}$ of the true contour, the error of the linear interpolation, so these sums approximate arc lengths of the contour. ^[8]

10. Visible points and the minimum

A contour point counts as visible from $A$ if the segment to it passes a sampled test: points about $0.5$ apart on it must not be higher than $f+10^{-10}$. Here $1\,152$ contour points pass from $A$ and $1\,328$ from $B$, and $AB$ fails. The program evaluates

$$|AP_i|+\min\bigl(s_{ij},\,L-s_{ij}\bigr)+|P_jB|$$

for all $1\,152\cdot1\,328=1\,529\,856$ pairs and outputs the smallest value. Each candidate is the length of a flight along the traced contour, and for a convex polygon with exact visibility the minimum would be the shortest flight round it: such a flight bends only at corners, ^[7] and a chord between non-neighbouring corners of a convex chain runs through the obstacle, so between its first and last bend the flight steps from corner to corner. The traced chain is convex only up to small wiggles from the interpolation, and visibility is sampled, so the result is an approximation; the last section compares it with an exact computation. The minimum is attained at $P_i\approx(48.49,624.00)$ and $P_j\approx(873.00,1536.05)$, next to $T_A$ and $T_B$.

How the Code Works

The C++ program ^[9] does its work in main. compute_heights fills heights with $h$ at $(i\Delta,j\Delta)$ under the index $i\,n+j$, where $n=\lfloor1600/\Delta\rfloor+1$ and $\Delta=1$ is the default of --grid-step (another step should divide $200$, as coordinates are truncated to grid indices). minimax_height is the search of Section 4: a std::priority_queue ordered by std::greater holds pairs of a label and an index, best holds the labels, outdated entries are skipped, and the function returns when $B$ comes out; the level is then checked against $h(A)$ and $h(B)$.

The byte array inside marks the heights above the level, the marching squares loop reads the tables kCaseCount and kCaseEdges, and interp places the crossing points, clamping $t$ to $[0,1]$. quantize makes a Key of two 64-bit integers with std::llround, and the std::unordered_map called adj lists the neighbours of each key. The program stops if a key does not have exactly two neighbours or a walk meets a visited key other than its start; each walk steps to the neighbour it did not come from, and the longest loop is kept. The warning about several loops under --validate compares the number of visited keys with the size of adj, which are always equal after all walks, so it never appears, although this level has two loops.

visible tests the points $a+(s/k)(b-a)$ for $s=1,\dots,k-1$, with $k=\max(1,\lfloor|ab|/\delta\rfloor)$ and $\delta=0.5$ the default of --vis-step; the ends are skipped, since $A$ and $B$ lie below the level and a traced point may lie slightly above it. The double loop over visA and visB evaluates the formula of Section 10 with std::min(arc, perimeter - arc), which does not depend on where the walk started or its direction. Further checks stop the program on an empty visible set or a non-finite result, and --validate compares the result with $|AB|\approx1697.06$; these are consistency checks, not error bounds.

The Python and Java versions follow the same steps with fixed parameters and without the C++ error checks; their walks stop at the first visited key. All three use 64-bit floating point. Python uses tuples as keys, a dictionary and heapq; Java runs the search inside solve with a PriorityQueue of Node objects compared by label only, so it breaks ties differently from the label-and-index order of C++ and Python, and uses a HashMap with its own Key class. Their key rounding treats exact halves differently, which does not matter here.

Complexity and Verification

With $n=1601$ the grid has $N=n^2=2\,563\,201$ points. The heights, marching squares and the hash map cost $O(N)$, the last in expectation, and the search $O(N\log N)$, since a point enters the binary heap at most once per neighbour. ^[3] A visibility test makes at most $1+D/\delta$ height evaluations, $D=1600\sqrt2$, so for $m$ contour points and the visible sets $V_A$, $V_B$ the time is

$$O\bigl(N\log N+m\,D/\delta+|V_A|\,|V_B|\bigr).$$

Here the visibility tests make $5\,409\,298$ height evaluations, most blocked segments being rejected early, and the pair scan makes one step per pair. The measured peak memory is about $49$ MB in C++, mostly heights and best, $120$ MB in Java and $340$ MB in Python; on one desktop computer the runs take about $0.05$ s, $0.13$ s and $2$ s.

All three programs print the same value. Independent numerical checks for this article, not part of the published solutions and not error bounds:

The cubic $c$, solved with 40-digit arithmetic, gives $y_0\approx895.483409899$ and $f_{\min}\approx10396.462193284$; grid searches with the steps $1$, $2$ and $4$ and four or eight neighbours all find the pass next to $S$ and $S'$.
A continuous computation finds the outer contour as a function $r(\varphi)$ of the angle around $M$ by bisection on $2^{18}$ rays, the tangent points, and the arc as a sum of chords. ^[8] With $2^{16}$, $2^{18}$ and $2^{20}$ rays the lengths agree to within $10^{-6}$, and the rounded result agrees with the programs.
The unrounded grid result is about $2.8\cdot10^{-4}$ longer. The exact contour at the program's level, cut off at the edge like the traced one, adds $3.5\cdot10^{-4}$; the traced chain is $0.8\cdot10^{-4}$ shorter than that exact arc, and its $94$ small wiggles among $1\,884$ corners matter less than $10^{-7}$. Such an error could change the third decimal only for an exact length just below a rounding boundary, which is not the case here.

Footnotes and References

The references below supply the background facts used above; the proofs and computations specific to this problem are given in the article.

Project Euler 262 — Mountain Range. The official statement gives h, A(200, 200), B(1400, 1400) and the area 0 ≤ x, y ≤ 1600, defines f_min as the minimum constant elevation allowing the trip, and asks for the length of the shortest path between A′ and B′ at that elevation, rounded to three decimal places.
Wikipedia — Widest path problem. The minimax or bottleneck shortest path problem asks for the path that minimises the maximum weight of its edges, and a modified Dijkstra's algorithm computes the bottlenecks from a start vertex.
Wikipedia — Dijkstra's algorithm. The algorithm repeatedly finalises the unvisited node with the smallest tentative label; nodes may be queued unconditionally and outdated entries skipped, and a binary heap gives the time O((V + E) log V).
Wikipedia — Mountain pass theorem. The theorem is explained by two low points separated by a ridge, whose lowest crossing, the mountain pass, is almost always a saddle point.
Wikipedia — Jordan curve theorem. Every simple closed curve divides the plane into an interior and an exterior, and every continuous path from one to the other meets the curve.
Wikipedia — Marching squares. The corners of a cell, compared with the level, give a 4-bit index into a table of 16 cases; linear interpolation places the contour on the sides, and the saddle cases are ambiguous.
Wikipedia — Visibility graph. Among polygonal obstacles a shortest path is straight except at obstacle vertices, so it is a shortest path in the graph of mutually visible vertices.
Wikipedia — Arc length. The length of a curve is the least upper bound of the lengths of inscribed polygons; in polar coordinates it is the integral of √(r² + (dr/dθ)²) dθ.
C++ — Euler262.cpp. The linked immutable C++ revision contains compute_heights, minimax_height, the marching squares tables, the loop walk, visible and the pair scan. It is the implementation record for this article.
Wikipedia — Logarithmically concave function. Every concave nonnegative function is log-concave, products of log-concave functions are log-concave, and log-concave functions have convex superlevel sets.

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());
    }
}