Problem 353: Risky Moon

Detailed mathematical approach

Problem Summary

A moon is the sphere \(C(r)\) of radius \(r\) centred at \((0,0,0)\); its stations are the points of \(C(r)\) with integer coordinates, \((0,0,r)\) the North Pole station and \((0,0,-r)\) the South Pole station. Every two stations are joined by a road, the shortest great arc through them. A road of length \(d\) has risk \(\left(\frac{d}{\pi r}\right)^2\), and a journey's risk is the sum of its roads' risks. The direct pole-to-pole journey has length \(\pi r\) and risk \(1\); the one through \((0,r,0)\) has the same length but risk \(0.5\), each half contributing \(\left(\frac12\right)^2\). The minimal risk between the poles is \(M(r)\). The statement gives \(M(7)=0.1784943998\) and asks for \(\sum_{n=1}^{15}M(2^n-1)\), both to ten decimal places. ^[1]

The whole solution in one sentence: a road's risk is the square of its central angle over \(\pi\), so \(M(r)\) is a shortest-path distance on the integer points of \(C(r)\), and a reflection in the equator plus the eight symmetries fixing the polar axis cut that graph — provably without loss — to a sixteenth of the sphere, small enough for a dense Dijkstra with a cutoff on latitude.

Mathematical Approach

1. The risk of a road

Let \(\theta(a,b)\in[0,\pi]\) be the angle two stations subtend at the centre. Both have length \(r\), so the dot product gives

$$\cos\theta(a,b)=\frac{a\cdot b}{r^2},$$

and the arc between them has length \(d=r\,\theta(a,b)\). ^[2] Since \(\pi r\) is the length of a half great circle,

$$w(a,b)=\left(\frac{d}{\pi r}\right)^2=\left(\frac{\theta(a,b)}{\pi}\right)^2 .$$

The radius has vanished: only angles matter, and as \(a\cdot b\) is an integer every cosine is rational with denominator \(r^2\). On \(C(3)\), \((0,0,3)\) and \((1,2,2)\) give \(\cos\theta=\frac69=\frac23\), \(\theta=0.841068670568\) and \(w=0.071674251557\); antipodal stations give risk \(1\), the direct pole-to-pole journey.

2. Why many short roads beat one long one

Let a pole-to-pole journey use \(k\) roads with angles \(\theta_1,\dots,\theta_k\). The angular distance is a metric (Section 9) and the poles are \(\pi\) apart, so \(\sum_i\theta_i\ge\pi\). Cauchy–Schwarz on \((\theta_1,\dots,\theta_k)\) and \((1,\dots,1)\) gives ^[5]

$$\Bigl(\sum_{i=1}^{k}\theta_i\Bigr)^{2}\le k\sum_{i=1}^{k}\theta_i^{2},$$

so the risk of the journey is

$$\sum_{i=1}^{k}\Bigl(\frac{\theta_i}{\pi}\Bigr)^{2}\ge\frac{1}{k\pi^{2}}\Bigl(\sum_i\theta_i\Bigr)^{2}\ge\frac1k,$$

with equality only when all angles equal \(\pi/k\). A journey of \(k\) roads is never safer than \(1/k\), and the quarter-arcs of the statement attain that limit at \(k=2\). Risk thus rewards splitting a journey into many equal pieces, but the pieces must land on stations, which are finite in number. Minimal risk is not minimal length: the journeys below have total angle \(1.0547\pi\) at \(r=3\) and \(1.1855\pi\) at \(r=15\).

3. How many stations a moon has

The stations are the integer solutions of \(x^2+y^2+z^2=r^2\) counted with signs and order, the sequence \(A016725(r)\). ^[6] It is \(6\) times a multiplicative function that is \(1\) on powers of \(2\), so it depends only on the odd part of \(r\); for an odd prime \(p\) it is \(6p\) when \(p\equiv1\pmod 4\) and \(6(p+2)\) when \(p\equiv3\pmod 4\). As \(2^n-1\) is prime only for prime \(n\), and \(2^{11}-1=23\cdot89\), among \(1\le n\le15\) exactly \(n=2,3,5,7,13\) give primes, ^[7] all \(3\) modulo \(4\), so those radii carry \(30\), \(54\), \(198\), \(774\) and \(49158\) stations. The counts for \(n=1,\dots,15\) are

$$\begin{aligned}&6,\ 30,\ 54,\ 150,\ 198,\ 918,\\ &774,\ 2550,\ 3942,\ 12870,\\ &13350,\ 59670,\ 49158,\\ &174150,\ 272646 .\end{aligned}$$

They need not grow with \(r\): \(C(127)\) carries \(774\) stations against \(918\) on \(C(63)\), and \(M(127)=0.046342185232\) is accordingly worse than \(M(63)=0.041515684149\). Nor do more stations always help, since different radii offer different directions: \(C(13)\) has \(78\) stations against \(54\) on \(C(7)\), yet \(M(13)=0.181756601531\) exceeds \(M(7)=0.178494399755\). Both the number of stations and where they sit matter.

4. A shortest-path problem

Let the stations be the vertices of a complete graph whose edge \(\{a,b\}\) has weight \(w(a,b)\ge0\). A journey's risk is the weight of the corresponding walk, so

$$M(r)=\operatorname{dist}\bigl((0,0,r),(0,0,-r)\bigr),$$

which Dijkstra's algorithm computes because no weight is negative. ^[4] It gives the start risk \(0\) and the others a tentative value, then repeatedly settles the unsettled station of least tentative risk and relaxes the roads out of it, lowering a neighbour's value whenever the route through the settled station is cheaper. Settling the least value is safe because any other route to that station begins with an unsettled value at least as large and then adds non-negative risk. With a plain array this costs \(\Theta(V^2)\), right for a graph with \(\Theta(V^2)\) roads.

At \(r=1\) the stations are the six axis points. Two distinct ones are perpendicular or antipodal, so every road costs \(\frac14\) or \(1\) and \(k\) roads cost at least \(k/4\): one road is the direct journey of risk \(1\), three or more cost at least \(\frac34\), and two perpendicular quarter-circles cost exactly \(\frac12\). Hence \(M(1)=\frac12\).

5. The reflection in the equator

Let \(\sigma(x,y,z)=(x,y,-z)\). It maps stations to stations, preserves angles and swaps the poles, so if \(d(p)\) is the minimal risk from the North Pole station to \(p\), the minimal risk from \(p\) to the South Pole station is \(d(\sigma p)\).

Every pole-to-pole journey has a first road leaving the closed northern hemisphere, from some \(u\) with \(z_u\ge0\) to some \(v\) with \(z_v\lt0\); cutting there and minimising each part gives the exact identity

$$M(r)=\min_{u,v}\ \bigl(d(u)+w(u,v)+d(\sigma v)\bigr).$$

The published programs evaluate only the symmetric case \(v=\sigma u\),

$$\widehat M(r)=\min_{p}\ \bigl(2\,d(p)+w(p,\sigma p)\bigr),$$

the risk of a real journey — travel optimally to \(p\), cross to \(\sigma p\), mirror the first half — so \(\widehat M(r)\ge M(r)\). The crossing road has a closed form: from \(p\cdot\sigma p=x^2+y^2-z^2=r^2-2z^2\),

$$\cos\theta(p,\sigma p)=1-2\Bigl(\frac zr\Bigr)^{2},$$

which is \(\cos2\alpha\) for \(\alpha=\arcsin(z/r)\), so \(\theta(p,\sigma p)=2\alpha\) when \(0\le z\le r\) and

$$w(p,\sigma p)=\Bigl(\frac{2}{\pi}\arcsin\frac zr\Bigr)^{2}.$$

For an equatorial \(p\) the crossing road degenerates and costs nothing. This is why the programs tabulate \(a_t=\arcsin(t/r)/\pi\) for \(0\le t\le r\).

6. Why the symmetric crossing loses nothing

Let \(u,q\) lie in the closed northern hemisphere with latitudes \(0\le\alpha\le\beta\le\frac\pi2\), let \(\lambda\) be their difference in longitude and put \(t=\cos\alpha\cos\beta\cos\lambda\) and \(s=\sin\alpha\sin\beta\), so that \(\cos\theta(u,q)=t+s\) and \(\cos\theta(u,\sigma q)=t-s\). Write \(g(c)=(\arccos c)^2\), so \(\pi^2w=g(\cos\theta)\). With \(\theta=\arccos c\),

$$\begin{aligned}g'(c)&=-\frac{2\theta}{\sin\theta},\\ g''(c)&=\frac{2(\sin\theta-\theta\cos\theta)}{\sin^3\theta}.\end{aligned}$$

That numerator vanishes at \(\theta=0\) and has derivative \(\theta\sin\theta\ge0\), so \(g''\ge0\): \(g\) is convex and \(g'\) increasing. Hence \(c\mapsto g(c-s)-g(c+s)\) does not increase, and \(t\le\cos\alpha\cos\beta\) because \(\cos\lambda\le1\), so

$$\begin{aligned}&g(t-s)-g(t+s)\\ &\ \ge g\bigl(\cos(\alpha+\beta)\bigr)-g\bigl(\cos(\beta-\alpha)\bigr)\\ &\ =(\alpha+\beta)^2-(\beta-\alpha)^2=4\alpha\beta,\end{aligned}$$

which is at least \(4\alpha^2=\pi^2w(u,\sigma u)\). The minimal-risk distance also obeys \(d(u)\le d(q)+w(u,q)\), and adding the two inequalities gives

$$\begin{aligned}2d(u)+w(u,\sigma u)\\ \le d(u)+d(q)+w(u,\sigma q).\end{aligned}$$

In the identity of Section 5 put \(q=\sigma v\), which is northern; the term becomes \(d(u)+w(u,\sigma q)+d(q)\), symmetric in \(u\) and \(q\) because \(\sigma\) is an isometry. Naming the station of smaller latitude \(u\), the display bounds it below by a symmetric candidate, so \(\widehat M(r)\le M(r)\) and therefore \(\widehat M(r)=M(r)\).

7. A worked example: \(r=3\)

\(C(3)\) has \(30\) stations: six axis points and the \(24\) sign-and-order variants of \((1,2,2)\). Four satisfy \(0\le x\le y\) and \(z\ge0\), and they are the only ones the program keeps:

• \((0,0,3)\), the North Pole station, \(d=0\);
• \((1,2,2)\), reached directly, \(d=0.071674251557\);
• \((2,2,1)\), where the direct road costs \(0.153528046876\) but the detour through \((1,2,2)\) costs \(0.071674251557\)\({}+0.022945591975\)\({}=0.094619843532\);
• \((0,3,0)\), where the direct road costs \(0.25\) but two roads through \((1,2,2)\) cost \(2\cdot0.071674251557\)\({}=0.143348503114\).

Twice each distance plus its crossing road gives the candidates \(1\), \(0.359163618137\), \(0.236045666446\) and \(0.286697006228\). The winner is \((2,2,1)\), whose crossing road has \(\cos\theta=\frac79\) and risk \(0.046805979382\):

$$\begin{aligned}M(3)&=2\cdot0.094619843532\\ &\quad+0.046805979382\\ &=0.236045666446,\end{aligned}$$

realised by

$$\begin{aligned}&(0,0,3)\to(1,2,2)\to(2,2,1)\\ &\to(2,2,-1)\to(1,2,-2)\\ &\to(0,0,-3),\end{aligned}$$

five roads, mirror-symmetric about the equator; an unreduced Dijkstra over all \(30\) stations returns the same value. The candidate \(1\) is the direct journey, and the journey of the statement through \((0,3,0)\) has risk \(0.5\), which the detour through \((1,2,2)\) improves to \(0.286697006228\).

8. One sixteenth of the sphere

The eight signed swaps of the first two coordinates, \((x,y)\mapsto(\pm x,\pm y)\) and \((x,y)\mapsto(\pm y,\pm x)\), map stations to stations and fix both poles. With \(\sigma\) they generate a group of order \(16\) whose fundamental domain is the wedge

$$W=\{(x,y,z): 0\le x\le y,\ z\ge0\}.$$

Off its boundary each orbit meets \(W\) once, so \(W\) holds about a sixteenth of the stations; the wedge sizes for \(n=1,\dots,15\) are

$$\begin{aligned}&2,\ 4,\ 5,\ 13,\ 14,\ 60,\ 50,\\ &169,\ 250,\ 808,\ 838,\\ &3738,\ 3074,\ 10888,\ 17042,\end{aligned}$$

and the ratio of stations to wedge size climbs to \(15.998\) at \(r=32767\), so the largest instance shrinks from \(272646\) vertices to \(17042\) and the work drops by about \(256\).

Computing \(d\) inside \(W\) alone also costs nothing. Fold a station \(p\) by taking absolute values and then, if needed, swapping the first two coordinates so they increase; call the result \(F(p)\in W\). Absolute values cannot decrease a dot product, nor can sorting the first two coordinates of both points the same way, since \(ac+bd-ad-bc=(b-a)(d-c)\ge0\) for \(a\le b\) and \(c\le d\). So \(F(u)\cdot F(v)\ge u\cdot v\) and \(w(F(u),F(v))\le w(u,v)\). As \(F\) keeps the norm and fixes the pole and every station of \(W\), folding a journey from the pole to a \(p\in W\) and deleting repeats gives a journey inside \(W\) to the same \(p\) of no greater risk. So \(d\) computed in \(W\) is the true \(d\), and with Section 6 the programs compute \(M(r)\) exactly.

9. A lower bound for an angle, for free

The spherical law of cosines reads \(\cos c=\cos a\cos b+\sin a\sin b\cos C\) for a spherical triangle with sides \(a,b,c\) and the angle \(C\) opposite \(c\). ^[3] As \(\cos C\ge-1\) it gives \(\cos c\ge\cos(a+b)\), hence \(c\le a+b\) whenever \(a+b\le\pi\); and if \(a+b\gt\pi\) then \(c\le\pi\lt a+b\) anyway. So the angular distance obeys the triangle inequality, as Section 2 assumed.

Call \(\varphi_p=\arcsin(z_p/r)\) the latitude of \(p\). Its angular distance to the North Pole station is \(\frac\pi2-\varphi_p\), so the triangle inequality through that station gives \(|\varphi_a-\varphi_b|\le\theta(a,b)\) and hence, with the table of Section 5,

$$(a_{z_a}-a_{z_b})^2\le w(a,b).$$

The \(x\)- and \(y\)-axes serve as poles equally well, so the same holds with \(x\) or \(y\) for \(z\), at a cost of two lookups and a multiplication each, with no \(\arccos\).

10. The cutoff and the coarse-to-fine passes

Suppose a journey of risk \(B\) is already known, so the answer is at most \(B\). Let the search have settled a station \(c\) at risk \(\rho\) and be about to relax a station \(j\). Any candidate reaching \(j\) along this road costs at least \(2(\rho+w(c,j))\ge2(\rho+f^2)\), with \(f\) one of the three differences of Section 9, so the road cannot improve \(B\) once \(f^2\ge B/2-\rho\). The programs discard it when \(f^2\gt B-2\rho\); for \(\rho\le B/2\) that threshold is the larger of the two, so they prune less than they safely could and never lose an improving road. When \(\rho\gt B/2\) the test discards every road out of \(c\), rightly, since any journey through \(c\) already costs more than \(B\).

As the wedge list is sorted by decreasing \(z\), the quantity \(a_{z_c}-a_{z_j}\) grows with \(j\) once \(j\) is past \(c\): after the first failure of the \(z\)-test there every later unsettled station fails it too, so the loop can stop.

The bound \(B\) comes from the passes themselves, which take every \(27\)th station of the wedge, then every \(9\)th, every \(3\)rd, and finally all; each sample starts at index \(0\) and contains the previous one, and a sample of fewer than two stations is skipped. The first pass that runs has the cutoff disabled and returns the risk of a real journey, a legitimate \(B\), and the nesting keeps that journey available later. At \(r=32767\) the samples hold \(632\), \(1894\), \(5681\) and \(17042\) stations, leave bounds \(5.190\), \(2.988\), \(1.696\) and \(1.000\) times the final value, and spend \(199\,396\), \(433\,933\), \(2\,461\,330\) and \(13\,046\,627\) evaluations of \(\arccos\) — \(16.1\) million against \(145.2\) million for one unpruned pass.

How the Code Works

edge_risk (edgeRisk in Java) implements Section 1: ^[8] dot product, division by \(r^2\), a clamp into \([-1,1]\), acos, division by \(\pi\), square. The clamp is purely defensive, as \(|a\cdot b|\le r^2\) and both numbers are exact in floating point for \(r\le32767\).

generate_reduced_points builds the table \(a_t\) and enumerates the sorted triples \(0\le x\le y\le z\) with \(x^2+y^2+z^2=r^2\) without a square root, keeping the residual h\({}=x^2+y^2+z^2-r^2\): raising \(y\) adds \(2y+1\), lowering \(z\) subtracts \(2z-1\). For each \(x\) it starts at \(y=x\), lowers \(z\) until h is not positive, keeps that \(z\) in z0 for the next \(x\), then walks \(y\) up while \(y\le z\), lowering \(z\) again when h turns positive and recording a hit when h is \(0\). Its condition \(3x^2\lt r^2\) loses nothing: \(x\le y\le z\) forces \(3x^2\le r^2\), and \(3x^2=r^2\) is impossible in positive integers. From each triple it emits the wedge stations of Section 8 — every distinct value once as the third coordinate, the other two increasing — the five branches (y == 0, x == 0, y == z, x == y, general) being that rule written out for the coincidences. A sort by decreasing \(z\), then \(y\), then \(x\), with duplicates dropped, puts the North Pole station at index \(0\).

dijkstra_pruned is the dense Dijkstra of Section 4 on that list with the mirror closing of Section 5. The arrays risk and done hold tentative values and settled flags, and \(2.0\) serves as "unreached" — safe, since every wedge station is northern, so the direct road from the pole costs at most \(\frac14\). When riskmax is positive the relaxation loop applies the cutoff of Section 10 with dmax\({}=\)riskmax\({}-2\rho\), testing the \(z\), \(x\) and \(y\) differences in that order and breaking out on the \(z\)-test once \(j\) is past the current station; when it is \(0\) every road is relaxed. A linear scan picks the unsettled station of least risk, and the routine returns the smallest \(2\cdot\)risk[i]\({}+(2a_{z_i})^2\) over its own sample — on the final, full pass that is \(M(r)\).

minimal_risk_fast generates the wedge and tries the strides \(27\), \(9\), \(3\) and \(1\), skipping a sample of fewer than two stations; each call to dijkstra_pruned takes the previous result as its bound, the first taking \(0\), which disables the cutoff. solve adds the fifteen independent values.

The C++ file also holds the reference implementation: generate_all_points lists every station by scanning \(x\) and \(y\) and testing whether \(r^2-x^2-y^2\) is a square, pole_indices locates the poles, dijkstra_complete runs an unreduced Dijkstra between them, and run_checkpoints requires \(M(7)\) within \(0.5\cdot10^{-10}\) of \(0.1784943998\) and the fast method within \(10^{-12}\) of dijkstra_complete at \(r=31\). main returns status \(2\) on failure, --skip-checkpoints omits both, and the Python and Java files carry no such checks.

The integer work differs: C++ keeps coordinates and \(r^2\) in int and h in long long, Java pairs int coordinates with long for \(r^2\) and h, Python integers are unbounded; the largest \(r^2=1073676289\) still fits a signed 32-bit integer. Risks are long double in C++, ordinary doubles in Java and Python. The three dedupe differently — HashSet, std::unique, a seen set — yet end with the same ordered list, and the Python cutoff keeps only the \(z\)-test, costing angle evaluations but changing no value.

Complexity and Verification

Generating the wedge is \(\Theta(r^2)\): the outer coordinate takes order \(r\) values and the inner walk order \(r\) steps for a constant fraction of them — measured \(0.218\,r^2\) steps at every radius, \(233\,658\,285\) at \(r=32767\), or \(0.31\) s. The arcsine table costs \(\Theta(r)\) and the sort \(O(V_W\log V_W)\) for a wedge of \(V_W\) stations. Each Dijkstra pass settles stations at one scan of the sample apiece, so a pass is \(O(V_W^2)\), and the final pass settles nearly all of the wedge, so the four are \(\Theta(V_W^2)\). The method is thus \(\Theta(r^2+V_W^2)\) in time and \(\Theta(V_W+r)\) in memory, road risks being computed when needed rather than stored; at \(r=32767\) the table holds \(32768\) entries against \(17042\) stations. The cutoff removes \(\arccos\) calls, not selection scans: \(16\,141\,286\) instead of \(145\,206\,361\) at \(r=32767\), \(27.4\) instead of \(216.9\) million over all fifteen radii. Those are the three-axis counts of C++ and Java; Python's \(z\)-only test needs \(38\,575\,270\) at \(r=32767\) for the same value.

Squaring the angle makes a short road's risk well conditioned. With \(w=(\theta/\pi)^2\) and \(\theta=\arccos c\),

$$\frac{dw}{dc}=-\frac{2\theta}{\pi^2\sin\theta},$$

which tends to \(-2/\pi^2\) as \(\theta\to0\): the sensitivity of \(\arccos\) near \(c=1\) ^[3] is cancelled by the vanishing derivative of the square. Both ends of a wedge road have non-negative coordinates, so \(c\ge0\), \(\theta\le\frac\pi2\) and \(|dw/dc|\le1/\pi\) throughout: one rounding of \(c\) moves a risk by at most \(7\cdot10^{-17}\). Against \(60\)-digit arithmetic the worst error over \(50\,000\) random wedge roads at \(r=32767\) was \(6.9\cdot10^{-17}\), and the closest pair there, \((8574,22362,22363)\) and \((8574,22363,22362)\) at the least squared chord length two distinct stations can have, \(2\), loses nothing: its risk \(1.9\cdot10^{-10}\) comes out right to \(5\cdot10^{-19}\). With at most \(827\) roads in a journey the accumulated error is of order \(10^{-12}\), far inside the \(5\cdot10^{-11}\) ten decimals demand. This is evidence rather than a certificate: it bounds neither the library's own error in \(\arccos\) nor the effect of near-ties in the comparisons.

Besides the two checkpoints in the C++ file, the following was verified for this article by computation.

• \(M(1)=\frac12\), \(M(3)=0.236045666446\), \(M(7)=0.178494399755\) — which rounds to the \(0.1784943998\) of the statement — \(M(15)=0.107147947962\) and \(M(31)=0.087445300344\), each from an unreduced Dijkstra over all stations as well as from the wedge-and-mirror method.
• For every radius \(r\le100\) an unreduced pole-to-pole Dijkstra agreed with the wedge-and-mirror value to within \(2.8\cdot10^{-17}\), and \(d\) inside the wedge agreed with \(d\) on the whole sphere at every wedge station — corroborating the proofs of Sections 6 and 8.
• For all fifteen radii the four pruned passes return exactly what one unpruned pass over the wedge returns, and the Python cutoff agrees with the C++ and Java one. Brute force confirms that generate_reduced_points emits exactly the stations of \(W\) for every \(r=2^n-1\) up to \(511\).
• Direct enumeration reproduced the station counts and agreed with \(6p\) and \(6(p+2)\) for every odd prime below \(100\), and the winning journeys use \(2\), \(5\), \(7\), \(15\), \(15\), \(35\), \(35\), \(68\), \(91\), \(151\), \(157\), \(331\), \(313\), \(585\) and \(827\) roads. The inequalities of Sections 6 and 9 held on \(280\,000\) random pairs of stations each, the latter on all three axes.

On one machine the C++ program needs \(0.95\) s with both checkpoints and \(0.67\) s without, Java \(3.1\) s and Python \(45\) s; all three print the same ten decimal places.

Footnotes and References

The references supply the background facts; the shortest-path reduction, the bound \(1/k\), the two reduction proofs, the latitude cutoff and the error analysis are derived here.

1. Project Euler 353 — Risky Moon. The official statement places stations at the integer points of the sphere C(r), joins every two by the shortest great arc, charges a road of length d the risk (d/(πr))² and a journey the sum of its roads' risks, notes a risk of 1 for the direct pole-to-pole journey against 0.5 through (0, r, 0) at the same length, gives M(7) = 0.1784943998, and asks for the sum of M(2ⁿ − 1), n = 1 to 15, as a.bcdefghijk.
2. Wikipedia — Great-circle distance. It calls Δσ the central angle between two points of a sphere, gives Δσ = arccos(n₁·n₂) for the unit normals and d = rΔσ for the arc length, and notes that only the arctangent form is well conditioned at all angles — the first two facts are Section 1, the last motivates the error analysis.
3. Wikipedia — Spherical law of cosines. It states cos c = cos a cos b + sin a sin b cos C for a spherical triangle with sides a, b, c and C opposite c, and warns that recovering a small c by inverting the cosine magnifies rounding errors. Section 9 derives the triangle inequality from it.
4. Wikipedia — Dijkstra's algorithm. It finds shortest paths from one source to every other node, requires non-negative edge weights, and runs in Θ(|V|²) when the unsettled vertices sit in an array and the minimum is taken by linear search — the version used here.
5. Wikipedia — Cauchy–Schwarz inequality. For real numbers it gives (Σ uᵢvᵢ)² ≤ (Σ uᵢ²)(Σ vᵢ²), with equality exactly when the sequences are linearly dependent. Section 2 takes v all ones, which yields the bound 1/k.
6. OEIS A016725 — integer solutions of x² + y² + z² = n². It counts them with zeros allowed and signs and order distinguished — 6, 6, 30, 6, 30, 30, 54, 6, 102, 30, 78, … — and gives a(n) = 6·b(n) with b multiplicative, b(2ᵉ) = 1, b(pᵉ) = pᵉ for p ≡ 1 (mod 4) and b(pᵉ) = pᵉ + 2(pᵉ − 1)/(p − 1) for p ≡ 3 (mod 4): the station counts of Section 3.
7. Wikipedia — Mersenne prime. It proves from the factorisation of 2^(ab) − 1 that 2ⁿ − 1 is prime only when n is prime, and gives 2¹¹ − 1 = 2047 = 23 × 89 as the smallest prime exponent with a composite Mersenne number: hence only n = 2, 3, 5, 7, 13 below 16.
8. C++ — Euler353.cpp. The linked immutable revision contains generate_reduced_points, dijkstra_pruned, minimal_risk_fast and the exact routines run by run_checkpoints. It is the implementation record for this article; the wider comparisons above are extra validation, not tests in that file.