Problem 353: Risky Moon

Project Euler Problem 353 Solution

EulerSolve provides an optimized solution for Project Euler Problem 353, Risky Moon, with C++, Python, Java, and a step-by-step mathematical explanation.

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....

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

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

Problem 353 source code

C++

#include <algorithm>
#include <array>
#include <atomic>
#include <cmath>
#include <cstdint>
#include <iomanip>
#include <iostream>
#include <limits>
#include <pthread.h>
#include <thread>
#include <vector>

namespace {

struct Point {
    int x;
    int y;
    int z;
};

constexpr long double kPi = 3.141592653589793238462643383279502884L;

long double edge_risk(const Point& a, const Point& b, int r) {
    const long double dot = static_cast<long double>(a.x) * b.x +
                            static_cast<long double>(a.y) * b.y +
                            static_cast<long double>(a.z) * b.z;
    long double c = dot / (static_cast<long double>(r) * static_cast<long double>(r));
    if (c > 1.0L) c = 1.0L;
    if (c < -1.0L) c = -1.0L;
    const long double t = std::acos(c) / kPi;
    return t * t;
}

std::vector<Point> generate_all_points(int r) {
    const int rr = r * r;
    std::vector<Point> pts;
    pts.reserve(220000);
    for (int x = -r; x <= r; ++x) {
        const int x2 = x * x;
        for (int y = -r; y <= r; ++y) {
            const int z2 = rr - x2 - y * y;
            if (z2 < 0) continue;
            const int z = static_cast<int>(std::llround(std::sqrt(static_cast<long double>(z2))));
            if (z * z != z2) continue;
            pts.push_back({x, y, z});
            if (z != 0) pts.push_back({x, y, -z});
        }
    }
    return pts;
}

std::pair<int, int> pole_indices(const std::vector<Point>& pts, int r) {
    int north = -1;
    int south = -1;
    for (int i = 0; i < static_cast<int>(pts.size()); ++i) {
        if (pts[static_cast<std::size_t>(i)].x == 0 && pts[static_cast<std::size_t>(i)].y == 0 &&
            pts[static_cast<std::size_t>(i)].z == r) {
            north = i;
        }
        if (pts[static_cast<std::size_t>(i)].x == 0 && pts[static_cast<std::size_t>(i)].y == 0 &&
            pts[static_cast<std::size_t>(i)].z == -r) {
            south = i;
        }
    }
    return {north, south};
}

long double dijkstra_complete(const std::vector<Point>& pts, int src, int dst, int r) {
    const int n = static_cast<int>(pts.size());
    std::vector<long double> dist(static_cast<std::size_t>(n),
                                  std::numeric_limits<long double>::infinity());
    std::vector<unsigned char> done(static_cast<std::size_t>(n), 0);
    dist[static_cast<std::size_t>(src)] = 0.0L;
    int u = src;
    while (true) {
        done[static_cast<std::size_t>(u)] = 1;
        if (u == dst) break;
        const long double du = dist[static_cast<std::size_t>(u)];
        for (int v = 0; v < n; ++v) {
            if (done[static_cast<std::size_t>(v)] || v == u) continue;
            const long double nd = du + edge_risk(pts[static_cast<std::size_t>(u)],
                                                  pts[static_cast<std::size_t>(v)], r);
            if (nd < dist[static_cast<std::size_t>(v)]) dist[static_cast<std::size_t>(v)] = nd;
        }
        int best = -1;
        long double bd = std::numeric_limits<long double>::infinity();
        for (int v = 0; v < n; ++v) {
            if (done[static_cast<std::size_t>(v)]) continue;
            if (dist[static_cast<std::size_t>(v)] < bd) {
                bd = dist[static_cast<std::size_t>(v)];
                best = v;
            }
        }
        if (best < 0) break;
        u = best;
    }
    return dist[static_cast<std::size_t>(dst)];
}

std::vector<Point> generate_reduced_points(int r, std::vector<long double>& asin_pi) {
    asin_pi.resize(static_cast<std::size_t>(r + 1));
    const long double inv_r = 1.0L / static_cast<long double>(r);
    for (int z = 0; z <= r; ++z) {
        asin_pi[static_cast<std::size_t>(z)] = std::asin(static_cast<long double>(z) * inv_r) / kPi;
    }

    std::vector<Point> pts;
    pts.reserve(30000);
    const int rr = r * r;

    int z0 = r;
    for (int x = 0; 3LL * x * x < rr; ++x) {
        int y = x;
        int z = z0;
        long long h = 1LL * x * x + 1LL * y * y + 1LL * z * z - rr;
        while (h > 0) {
            h -= 2LL * z - 1LL;
            --z;
        }
        z0 = z;

        while (y <= z) {
            if (h == 0) {
                if (y == 0) {
                    pts.push_back({0, 0, r});
                    pts.push_back({0, r, 0});
                } else if (x == 0) {
                    pts.push_back({0, y, z});
                    pts.push_back({0, z, y});
                    pts.push_back({y, z, 0});
                } else if (y == z) {
                    pts.push_back({x, y, y});
                    pts.push_back({y, y, x});
                } else if (x == y) {
                    pts.push_back({x, x, z});
                    pts.push_back({x, z, x});
                } else {
                    pts.push_back({x, y, z});
                    pts.push_back({x, z, y});
                    pts.push_back({y, z, x});
                }
            }

            h += 2LL * y + 1LL;
            ++y;
            if (h > 0) {
                h -= 2LL * z - 1LL;
                --z;
            }
        }
    }

    std::sort(pts.begin(), pts.end(), [](const Point& a, const Point& b) {
        if (a.z != b.z) return a.z > b.z;
        if (a.y != b.y) return a.y > b.y;
        return a.x > b.x;
    });
    pts.erase(std::unique(pts.begin(), pts.end(),
                          [](const Point& a, const Point& b) {
                              return a.x == b.x && a.y == b.y && a.z == b.z;
                          }),
              pts.end());
    return pts;
}

long double dijkstra_pruned(const std::vector<Point>& pts,
                            const std::vector<long double>& asin_pi,
                            int r,
                            long double riskmax) {
    const int n = static_cast<int>(pts.size());
    if (n < 2) return std::numeric_limits<long double>::infinity();

    std::vector<long double> risk(static_cast<std::size_t>(n), 2.0L);
    std::vector<unsigned char> done(static_cast<std::size_t>(n), 0);

    int cur = 0;
    risk[0] = 0.0L;
    while (true) {
        done[static_cast<std::size_t>(cur)] = 1;
        const Point& pc = pts[static_cast<std::size_t>(cur)];
        const long double rc = risk[static_cast<std::size_t>(cur)];

        if (riskmax > 0.0L) {
            const long double dmax = riskmax - 2.0L * rc;
            for (int j = 0; j < n; ++j) {
                if (done[static_cast<std::size_t>(j)]) continue;
                const Point& pj = pts[static_cast<std::size_t>(j)];
                long double f = asin_pi[static_cast<std::size_t>(pc.z)] -
                                asin_pi[static_cast<std::size_t>(pj.z)];
                if (f * f > dmax) {
                    if (j > cur) break;
                    continue;
                }
                f = asin_pi[static_cast<std::size_t>(pc.x)] -
                    asin_pi[static_cast<std::size_t>(pj.x)];
                if (f * f > dmax) continue;
                f = asin_pi[static_cast<std::size_t>(pc.y)] -
                    asin_pi[static_cast<std::size_t>(pj.y)];
                if (f * f > dmax) continue;

                const long double nd = rc + edge_risk(pc, pj, r);
                if (nd < risk[static_cast<std::size_t>(j)]) {
                    risk[static_cast<std::size_t>(j)] = nd;
                }
            }
        } else {
            for (int j = 0; j < n; ++j) {
                if (done[static_cast<std::size_t>(j)] || j == cur) continue;
                const long double nd = rc + edge_risk(pc, pts[static_cast<std::size_t>(j)], r);
                if (nd < risk[static_cast<std::size_t>(j)]) {
                    risk[static_cast<std::size_t>(j)] = nd;
                }
            }
        }

        int next = -1;
        long double best = 2.0L;
        for (int j = 0; j < n; ++j) {
            if (done[static_cast<std::size_t>(j)]) continue;
            if (risk[static_cast<std::size_t>(j)] < best) {
                best = risk[static_cast<std::size_t>(j)];
                next = j;
            }
        }
        if (next < 0) break;
        cur = next;
    }

    long double ans = 2.0L;
    for (int i = 0; i < n; ++i) {
        const long double e = 2.0L * asin_pi[static_cast<std::size_t>(pts[static_cast<std::size_t>(i)].z)];
        const long double cand = 2.0L * risk[static_cast<std::size_t>(i)] + e * e;
        if (cand < ans) ans = cand;
    }
    return ans;
}

std::vector<Point> downsample(const std::vector<Point>& pts, int stride) {
    std::vector<Point> out;
    out.reserve(pts.size() / static_cast<std::size_t>(stride) + 2);
    for (std::size_t i = 0; i < pts.size(); i += static_cast<std::size_t>(stride)) {
        out.push_back(pts[i]);
    }
    return out;
}

long double minimal_risk_fast(int r) {
    std::vector<long double> asin_pi;
    const std::vector<Point> pts = generate_reduced_points(r, asin_pi);
    if (pts.size() < 2) return std::numeric_limits<long double>::infinity();

    long double d = 0.0L;
    for (int stride : {27, 9, 3, 1}) {
        std::vector<Point> sample = (stride == 1) ? pts : downsample(pts, stride);
        if (sample.size() < 2) continue;
        d = dijkstra_pruned(sample, asin_pi, r, d);
    }
    return d;
}

bool approx_equal_10dp(long double a, long double b) {
    return std::fabsl(a - b) < 0.5e-10L;
}

bool run_checkpoints() {
    const long double m7 = minimal_risk_fast(7);
    if (!approx_equal_10dp(m7, 0.1784943998L)) {
        std::cerr << "Checkpoint failed: M(7)\n";
        return false;
    }

    const int r = 31;
    const auto pts = generate_all_points(r);
    const auto [north, south] = pole_indices(pts, r);
    const long double exact = dijkstra_complete(pts, north, south, r);
    const long double fast = minimal_risk_fast(r);
    if (std::fabsl(exact - fast) > 1e-12L) {
        std::cerr << "Checkpoint failed: M(31) exact mismatch\n";
        return false;
    }

    return true;
}

long double solve() {
    struct Task {
        std::atomic<int>* next_n;
        std::array<long double, 16>* values;
    };

    auto worker = [](void* arg) -> void* {
        auto* task = static_cast<Task*>(arg);
        while (true) {
            const int n = task->next_n->fetch_add(1);
            if (n > 15) break;
            const int r = (1 << n) - 1;
            (*task->values)[static_cast<std::size_t>(n)] = minimal_risk_fast(r);
        }
        return nullptr;
    };

    unsigned hw = std::thread::hardware_concurrency();
    if (hw == 0U) hw = 1U;
    const unsigned thread_count = std::min<unsigned>(8U, std::min<unsigned>(15U, hw));

    std::array<long double, 16> values{};
    std::atomic<int> next_n{1};
    Task task{&next_n, &values};
    std::vector<pthread_t> threads(thread_count);
    for (unsigned t = 0; t < thread_count; ++t) {
        pthread_create(&threads[static_cast<std::size_t>(t)], nullptr, worker, &task);
    }
    for (unsigned t = 0; t < thread_count; ++t) {
        pthread_join(threads[static_cast<std::size_t>(t)], nullptr);
    }

    long double sum = 0.0L;
    for (int n = 1; n <= 15; ++n) {
        sum += values[static_cast<std::size_t>(n)];
    }
    return sum;
}

}  // namespace

int main(int argc, char** argv) {
    bool skip_checkpoints = false;
    for (int i = 1; i < argc; ++i) {
        std::string arg(argv[i]);
        if (arg == "--skip-checkpoints") {
            skip_checkpoints = true;
        } else {
            std::cerr << "Unknown argument: " << arg << '\n';
            return 1;
        }
    }

    if (!skip_checkpoints && !run_checkpoints()) {
        return 2;
    }

    const long double answer = solve();
    std::cout << std::fixed << std::setprecision(10) << answer << '\n';
    return 0;
}

Python

import math
import heapq

def solve():
    PI = math.pi

    def edge_risk(a, b, r):
        dot = a[0]*b[0] + a[1]*b[1] + a[2]*b[2]
        c = dot / (r * r)
        c = max(-1.0, min(1.0, c))
        t = math.acos(c) / PI
        return t * t

    def generate_reduced_points(r):
        asin_pi = [math.asin(z / r) / PI for z in range(r + 1)]
        rr = r * r
        pts = []
        z0 = r
        x = 0
        while 3 * x * x < rr:
            y = x
            z = z0
            h = x*x + y*y + z*z - rr
            while h > 0:
                h -= 2*z - 1
                z -= 1
            z0 = z
            while y <= z:
                if h == 0:
                    if y == 0:
                        pts.append((0, 0, r))
                        pts.append((0, r, 0))
                    elif x == 0:
                        pts.append((0, y, z))
                        pts.append((0, z, y))
                        pts.append((y, z, 0))
                    elif y == z:
                        pts.append((x, y, y))
                        pts.append((y, y, x))
                    elif x == y:
                        pts.append((x, x, z))
                        pts.append((x, z, x))
                    else:
                        pts.append((x, y, z))
                        pts.append((x, z, y))
                        pts.append((y, z, x))
                h += 2*y + 1
                y += 1
                if h > 0:
                    h -= 2*z - 1
                    z -= 1
            x += 1
        pts.sort(key=lambda p: (-p[2], -p[1], -p[0]))
        seen = set()
        unique = []
        for p in pts:
            if p not in seen:
                seen.add(p)
                unique.append(p)
        return unique, asin_pi

    def dijkstra_pruned(pts, asin_pi, r, riskmax):
        n = len(pts)
        if n < 2: return float('inf')
        risk = [2.0] * n
        done = [False] * n
        cur = 0
        risk[0] = 0.0
        while True:
            done[cur] = True
            pc = pts[cur]
            rc = risk[cur]
            for j in range(n):
                if done[j]: continue
                pj = pts[j]
                if riskmax > 0:
                    dmax = riskmax - 2.0 * rc
                    f = asin_pi[pc[2]] - asin_pi[pj[2]]
                    if f*f > dmax:
                        if j > cur: break
                        continue
                nd = rc + edge_risk(pc, pj, r)
                if nd < risk[j]:
                    risk[j] = nd
            nxt = -1
            best = 2.0
            for j in range(n):
                if not done[j] and risk[j] < best:
                    best = risk[j]
                    nxt = j
            if nxt < 0: break
            cur = nxt
        ans = 2.0
        for i in range(n):
            e = 2.0 * asin_pi[pts[i][2]]
            cand = 2.0 * risk[i] + e * e
            if cand < ans: ans = cand
        return ans

    def minimal_risk_fast(r):
        pts, asin_pi = generate_reduced_points(r)
        if len(pts) < 2: return float('inf')
        d = 0.0
        for stride in [27, 9, 3, 1]:
            sample = pts[::stride] if stride > 1 else pts
            if len(sample) < 2: continue
            d = dijkstra_pruned(sample, asin_pi, r, d)
        return d

    total = 0.0
    for n in range(1, 16):
        r = (1 << n) - 1
        total += minimal_risk_fast(r)
    return f"{total:.10f}"

if __name__ == '__main__':
    print(solve())

Java

import java.util.*;
import java.util.stream.IntStream;

public class Euler353 {

    static class Point implements Comparable<Point> {
        int x, y, z;

        Point(int x, int y, int z) {
            this.x = x;
            this.y = y;
            this.z = z;
        }

        @Override
        public int compareTo(Point o) {
            if (this.z != o.z)
                return Integer.compare(o.z, this.z);
            if (this.y != o.y)
                return Integer.compare(o.y, this.y);
            return Integer.compare(o.x, this.x);
        }

        @Override
        public boolean equals(Object o) {
            if (this == o)
                return true;
            if (o == null || getClass() != o.getClass())
                return false;
            Point p = (Point) o;
            return x == p.x && y == p.y && z == p.z;
        }

        @Override
        public int hashCode() {
            return Objects.hash(x, y, z);
        }
    }

    static double edgeRisk(Point a, Point b, int r) {
        double dot = (double) a.x * b.x + (double) a.y * b.y + (double) a.z * b.z;
        double c = dot / ((double) r * r);
        if (c > 1.0)
            c = 1.0;
        if (c < -1.0)
            c = -1.0;
        double t = Math.acos(c) / Math.PI;
        return t * t;
    }

    static class ReducedResult {
        List<Point> pts;
        double[] asinPi;
    }

    static ReducedResult generateReducedPoints(int r) {
        ReducedResult res = new ReducedResult();
        res.asinPi = new double[r + 1];
        double invR = 1.0 / r;
        for (int z = 0; z <= r; z++) {
            res.asinPi[z] = Math.asin((double) z * invR) / Math.PI;
        }

        Set<Point> set = new HashSet<>();
        long rr = (long) r * r;
        int z0 = r;

        for (int x = 0; 3L * x * x < rr; x++) {
            int y = x;
            int z = z0;
            long h = (long) x * x + (long) y * y + (long) z * z - rr;
            while (h > 0) {
                h -= 2L * z - 1L;
                z--;
            }
            z0 = z;
            while (y <= z) {
                if (h == 0) {
                    if (y == 0) {
                        set.add(new Point(0, 0, r));
                        set.add(new Point(0, r, 0));
                    } else if (x == 0) {
                        set.add(new Point(0, y, z));
                        set.add(new Point(0, z, y));
                        set.add(new Point(y, z, 0));
                    } else if (y == z) {
                        set.add(new Point(x, y, y));
                        set.add(new Point(y, y, x));
                    } else if (x == y) {
                        set.add(new Point(x, x, z));
                        set.add(new Point(x, z, x));
                    } else {
                        set.add(new Point(x, y, z));
                        set.add(new Point(x, z, y));
                        set.add(new Point(y, z, x));
                    }
                }
                h += 2L * y + 1L;
                y++;
                if (h > 0) {
                    h -= 2L * z - 1L;
                    z--;
                }
            }
        }
        res.pts = new ArrayList<>(set);
        Collections.sort(res.pts);
        return res;
    }

    static double dijkstraPruned(List<Point> pts, double[] asinPi, int r, double riskmax) {
        int n = pts.size();
        if (n < 2)
            return Double.POSITIVE_INFINITY;

        double[] risk = new double[n];
        boolean[] done = new boolean[n];
        Arrays.fill(risk, 2.0);

        int cur = 0;
        risk[0] = 0.0;

        while (true) {
            done[cur] = true;
            Point pc = pts.get(cur);
            double rc = risk[cur];

            if (riskmax > 0.0) {
                double dmax = riskmax - 2.0 * rc;
                for (int j = 0; j < n; j++) {
                    if (done[j])
                        continue;
                    Point pj = pts.get(j);
                    double f = asinPi[pc.z] - asinPi[pj.z];
                    if (f * f > dmax) {
                        if (j > cur)
                            break;
                        continue;
                    }
                    f = asinPi[pc.x] - asinPi[pj.x];
                    if (f * f > dmax)
                        continue;
                    f = asinPi[pc.y] - asinPi[pj.y];
                    if (f * f > dmax)
                        continue;

                    double nd = rc + edgeRisk(pc, pj, r);
                    if (nd < risk[j])
                        risk[j] = nd;
                }
            } else {
                for (int j = 0; j < n; j++) {
                    if (done[j] || j == cur)
                        continue;
                    double nd = rc + edgeRisk(pc, pts.get(j), r);
                    if (nd < risk[j])
                        risk[j] = nd;
                }
            }

            int next = -1;
            double best = 2.0;
            for (int j = 0; j < n; j++) {
                if (!done[j] && risk[j] < best) {
                    best = risk[j];
                    next = j;
                }
            }
            if (next < 0)
                break;
            cur = next;
        }

        double ans = 2.0;
        for (int i = 0; i < n; i++) {
            double e = 2.0 * asinPi[pts.get(i).z];
            double cand = 2.0 * risk[i] + e * e;
            if (cand < ans)
                ans = cand;
        }
        return ans;
    }

    static List<Point> downsample(List<Point> pts, int stride) {
        List<Point> out = new ArrayList<>();
        for (int i = 0; i < pts.size(); i += stride)
            out.add(pts.get(i));
        return out;
    }

    static double minimalRiskFast(int r) {
        ReducedResult res = generateReducedPoints(r);
        if (res.pts.size() < 2)
            return Double.POSITIVE_INFINITY;

        double d = 0.0;
        int[] strides = { 27, 9, 3, 1 };
        for (int stride : strides) {
            List<Point> sample = (stride == 1) ? res.pts : downsample(res.pts, stride);
            if (sample.size() < 2)
                continue;
            d = dijkstraPruned(sample, res.asinPi, r, d);
        }
        return d;
    }

    static double solve() {
        return IntStream.rangeClosed(1, 15).parallel()
                .mapToDouble(n -> minimalRiskFast((1 << n) - 1))
                .sum();
    }

    public static void main(String[] args) {
        System.out.printf(Locale.US, "%.10f\n", solve());
    }
}