Detailed mathematical approach

Problem Summary

For each \(n\), the problem sets \(r=2^n-1\) and considers the integer points on the sphere

$S_r=\{(x,y,z)\in \mathbb{Z}^3 : x^2+y^2+z^2=r^2\}.$

The north and south poles are \(N=(0,0,r)\) and \(S=(0,0,-r)\). A move from one lattice point to another has a risk equal to the square of the normalized central angle between them. If \(M(r)\) denotes the minimum total risk of a path from \(N\) to \(S\), then the Project Euler target is

$\sum_{n=1}^{15} M(2^n-1).$

The local source files solve this by combining an exact dense-graph model with a symmetry-reduced fast solver.

Mathematical Approach

Step 1: Edge Risk From Spherical Geometry

If \(a,b\in S_r\), then \(\lVert a\rVert=\lVert b\rVert=r\). Let \(\theta(a,b)\) be the central angle. By the dot-product identity on a sphere,

$\cos\theta(a,b)=\frac{a\cdot b}{r^2},\qquad \theta(a,b)=\arccos\left(\frac{a\cdot b}{r^2}\right).$

The implementations normalize this angle by \(\pi\), clamp the cosine into \([-1,1]\) for numerical safety, and define

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

This is exactly the edge_risk/edgeRisk formula in the C++, Python, and Java solutions.

Step 2: Turn the Moon Into a Dense Weighted Graph

Because the code allows a direct move between any two lattice points, \(S_r\) becomes a complete weighted graph. For a path

$P=(p_0,p_1,\dots,p_k),\qquad p_0=N,\quad p_k=S,$

the total risk is

$R(P)=\sum_{i=0}^{k-1} w(p_i,p_{i+1}).$

Therefore

$M(r)=\min_{P:N\leadsto S} R(P),$

so the mathematical core is a shortest-path problem. The exact C++ checkpoint path uses Dijkstra on the full point set generated by generate_all_points.

Step 3: A Small Worked Example

When \(r=1\), the only lattice points are the six axis points. The direct north-to-south jump has \(\theta=\pi\), so its risk is \(1\). Going through any equatorial axis point uses two quarter-turns, each with normalized angle \(1/2\), hence

$M(1)=2\left(\frac{1}{2}\right)^2=\frac{1}{2}.$

This illustrates why splitting a long move into several shorter moves can reduce total risk: the square function strongly penalizes large angles.

Step 4: Mirror Formula Used by the Fast Solver

The optimized algorithm computes shortest risks from the north pole to a reduced set of representative points in the upper half of the sphere. Suppose one such point is

$p=(x,y,z),\qquad z\ge 0,$

and let its mirror across the equatorial plane be

$p'=(x,y,-z).$

If \(d(p)\) is the minimal risk from \(N\) to \(p\), symmetry gives the same value from \(p'\) to \(S\). So a full north-to-south candidate is

$2d(p)+w(p,p').$

Now

$p\cdot p'=x^2+y^2-z^2=r^2-2z^2,$

hence

$\cos\theta=\frac{r^2-2z^2}{r^2}=1-2\left(\frac{z}{r}\right)^2.$

Writing \(\alpha=\arcsin(z/r)\), we have \(\cos\theta=\cos(2\alpha)\), so for \(0\le z\le r\),

$\theta=2\arcsin\left(\frac{z}{r}\right).$

Therefore the mirror-bridge risk is

$w(p,p')=\left(\frac{2\arcsin(z/r)}{\pi}\right)^2.$

The fast solver finally minimizes

$\boxed{M(r)=\min_{p}\left(2d(p)+\left(\frac{2\arcsin(z_p/r)}{\pi}\right)^2\right).}$

Step 5: Reduced Point Set and Conservative Pruning

The function generate_reduced_points does not enumerate every sign and permutation copy. Instead it scans ordered nonnegative triples with

$0\le x\le y\le z,\qquad x^2+y^2+z^2=r^2,$

then emits only the coordinate permutations needed by the chosen symmetry model. The outer bound \(3x^2<r^2\) comes from \(x\le y\le z\). After sorting and deduplication, this reduced list becomes the vertex set of the fast search.

The code also precomputes

$a_t=\frac{\arcsin(t/r)}{\pi}\qquad (0\le t\le r).$

During Dijkstra, if the current best full north-to-south bound is \(B\) and the current one-sided risk is \(d_c\), the implementation uses the conservative remaining budget \(B-2d_c\). Cheap lower bounds based on differences of the precomputed \(a_t\) values are tested before the expensive acos call. In C++ and Java this test is applied to the \(z\), \(x\), and \(y\) coordinates; the Python version keeps the same overall structure but only uses the \(z\)-based bound.

Step 6: Coarse-to-Fine Search and Checkpoints

The fast solver runs four passes on the reduced point set with strides \(27,9,3,1\). Each coarse pass gives an upper bound that makes the next pass cheaper. This is why minimal_risk_fast repeatedly calls the pruned Dijkstra routine with progressively denser samples.

The C++ file validates the method in two ways. First, it checks the known benchmark

$M(7)=0.1784943998.$

Second, for \(r=31\) it runs both the exact complete-graph Dijkstra and the reduced fast solver and requires agreement to within \(10^{-12}\). Those checks are important because the fast method is an optimization of the exact graph formulation, not a different mathematical model.

How the Code Works

The C++ source contains both the exact validator and the optimized solver. The Python and Java files implement the same reduced-point mathematics directly. All three versions share the same ingredients: the edge-risk formula, precomputed \(\arcsin(z/r)/\pi\) values, representative-point generation, and the mirror closing formula.

Implementation details differ slightly. C++ parallelizes the 15 radii with up to eight pthread workers; Java sums them with a parallel IntStream; Python evaluates them serially. The exact benchmark routines generate_all_points, pole_indices, and dijkstra_complete appear only in the C++ source and are used for correctness checks rather than for the full run.

Complexity Analysis

If \(V_r=|S_r|\), the exact complete-graph Dijkstra uses \(O(V_r^2)\) time and \(O(V_r)\) memory because every relaxation step scans all remaining vertices. That is acceptable for checkpoints such as \(r=31\), but not for the largest radii in the final sum.

The fast solver still performs dense relaxations, so its worst-case cost is quadratic in the reduced vertex count. However, the reduced set is much smaller than the full lattice-point set, the coarse-to-fine schedule reuses increasingly strong upper bounds, and many candidate edges are skipped by the conservative pruning tests. In practice these optimizations are what make the \(n=1,\dots,15\) computation feasible.

Footnotes and References

1. Problem page: https://projecteuler.net/problem=353
2. Great-circle distance and central angle: Wikipedia — Great-circle distance
3. Dijkstra's algorithm: Wikipedia — Dijkstra's algorithm
4. Lattice points on spheres and three-square representations: Wikipedia — Sum of three squares theorem
5. Dot product and spherical law of cosines: Wikipedia — Spherical law of cosines