Problem 532: Nanobots on Geodesics

Project Euler Problem 532 Solution

EulerSolve provides an optimized solution for Project Euler Problem 532, Nanobots on Geodesics, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary In the symmetric nanobot configuration, every bot follows the same geodesic-type path on the sphere, so the total travelled distance is $T(n)=nL(n)$, where $L(n)$ is the length assigned to one bot when there are $n$ bots. The C++, Python, and Java implementations do not simulate the full multi-bot motion directly. They reduce the geometry to one numerical integral depending only on $n$, then search for the smallest integer $n\ge 3$ such that $$L(n) > 1000.$$ After locating that first threshold crossing, the reported value is the total distance $$T(n)=nL(n),$$ printed to two decimal places. Mathematical Approach Set $$s=\sin\left(\frac{\pi}{n}\right),\qquad u_0=0.999.$$ After the substitution $u=\sin\theta$, where $\theta$ is the colatitude from the geometric setup, the implementations evaluate $$L(n)=\frac{1}{s}\int_0^{u_0}\frac{\sqrt{1-(su)^2}}{1-u^2}\,du.$$ The rest of the method consists of evaluating this integral accurately and then searching the first $n$ for which its value exceeds the target. Step 1: Reduce the $n$-bot motion to one representative path The rotational symmetry means every bot contributes the same path length....

Detailed mathematical approach

Problem Summary

In the symmetric nanobot configuration, every bot follows the same geodesic-type path on the sphere, so the total travelled distance is $T(n)=nL(n)$, where $L(n)$ is the length assigned to one bot when there are $n$ bots.

The C++, Python, and Java implementations do not simulate the full multi-bot motion directly. They reduce the geometry to one numerical integral depending only on $n$, then search for the smallest integer $n\ge 3$ such that

$$L(n) > 1000.$$

After locating that first threshold crossing, the reported value is the total distance

$$T(n)=nL(n),$$

printed to two decimal places.

Mathematical Approach

Set

$$s=\sin\left(\frac{\pi}{n}\right),\qquad u_0=0.999.$$

After the substitution $u=\sin\theta$, where $\theta$ is the colatitude from the geometric setup, the implementations evaluate

$$L(n)=\frac{1}{s}\int_0^{u_0}\frac{\sqrt{1-(su)^2}}{1-u^2}\,du.$$

The rest of the method consists of evaluating this integral accurately and then searching the first $n$ for which its value exceeds the target.

Step 1: Reduce the $n$-bot motion to one representative path

The rotational symmetry means every bot contributes the same path length. Therefore the global quantity of interest is obtained from a single scalar function $L(n)$, and the full answer is simply

$$T(n)=nL(n).$$

This is why the implementations can ignore the detailed positions of the other bots once the single representative geodesic has been derived.

Step 2: Express the path length as a one-dimensional integral

The dependence on $n$ enters through the spacing parameter $s=\sin(\pi/n)$. With the change of variable $u=\sin\theta$, the arc-length formula becomes

$$L(n)=\int_0^{u_0}\frac{\sqrt{\frac{1}{s^2}-u^2}}{1-u^2}\,du=\frac{1}{s}\int_0^{u_0}\frac{\sqrt{1-(su)^2}}{1-u^2}\,du.$$

This form isolates all $n$-dependence in a single parameter. The denominator $1-u^2$ becomes small near $u=1$, so the implementations evaluate the truncated model with the fixed cutoff $u_0=0.999$; that is the exact numerical quantity used for the threshold test.

Step 3: Evaluate the integral with adaptive Simpson refinement

For $f(u)=\sqrt{1-(su)^2}/(1-u^2)$, Simpson's rule on an interval $[a,b]$ is

$$S(a,b)=\frac{b-a}{6}\left(f(a)+4f\left(\frac{a+b}{2}\right)+f(b)\right).$$

The interval is bisected at $m=(a+b)/2$. The implementations compare the coarse estimate $S(a,b)$ with the refined estimate $S(a,m)+S(m,b)$, using

$$\Delta=S(a,m)+S(m,b)-S(a,b)$$

as the local error signal. When $|\Delta|\le 15\varepsilon$, the accepted corrected value is

$$S(a,m)+S(m,b)+\frac{\Delta}{15}.$$

Otherwise both halves are refined recursively. The tolerance is $10^{-13}$, and the recursion depth cap is $30$.

Step 4: Show that the computed length is monotone in $n$

For fixed $u\in[0,u_0]$, the integrand can be rewritten as

$$\frac{1}{s}\frac{\sqrt{1-(su)^2}}{1-u^2}=\frac{\sqrt{\frac{1}{s^2}-u^2}}{1-u^2}.$$

As $n$ increases, $s=\sin(\pi/n)$ decreases, so $1/s^2$ increases. Therefore the integrand increases pointwise on the whole integration interval, and the computed $L(n)$ is increasing as well. For large $n$, $s\sim \pi/n$, so $L(n)$ grows on the order of $1/s$, hence roughly linearly in $n$. That guarantees the existence of a first threshold index

$$n_*=\min\{n\ge 3:L(n)>1000\}.$$

Step 5: Locate the first crossing with doubling and binary search

Because $L(n)$ is monotone, the search is standard. Start from $n=3$, repeatedly double an upper bound until the computed length exceeds $1000$, and then binary-search inside that bracket. Once $n_*$ is found, the final reported quantity is

$$T(n_*)=n_*L(n_*).$$

Worked Example: $n=3$

For three bots,

$$s=\sin\left(\frac{\pi}{3}\right)=\frac{\sqrt{3}}{2},$$

so the implemented integral becomes

$$L(3)=\frac{2}{\sqrt{3}}\int_0^{0.999}\frac{\sqrt{1-\frac{3u^2}{4}}}{1-u^2}\,du.$$

Numerically, the implementation obtains

$$L(3)\approx 2.84.$$

This is the checkpoint used by the C++ version. The corresponding total distance is

$$T(3)\approx 3\times 2.84=8.52,$$

which is still far below $1000$, so the threshold search must continue to much larger $n$.

How the Code Works

The C++, Python, and Java implementations all follow the same numerical pipeline. For a chosen $n$, they compute $s=\sin(\pi/n)$, evaluate the integrand at the two endpoints and the midpoint of $[0,0.999]$, build the initial Simpson estimate, and recursively refine only the subintervals whose local correction is too large.

Each successful integral evaluation returns one value $L(n)$. Around that numerical kernel, the implementation first doubles an upper bound until $L(n)>1000$, then uses binary search to isolate the smallest valid integer $n$. The printed answer is the product $nL(n)$, rounded to two decimal places.

The C++ implementation also performs explicit sanity checks: it verifies the $n=3$ checkpoint and confirms that the final answer is truly the first threshold crossing by checking both $L(n)$ and $L(n-1)$.

Complexity Analysis

Let $Q(n)$ be the number of integrand evaluations used by the adaptive Simpson routine for one input $n$. The threshold search requires $O(\log n_*)$ length evaluations, so the total running time is

$$O\!\left(Q_{\max}\log n_*\right),\qquad Q_{\max}=\max_{3\le n\le n_*}Q(n).$$

The memory usage is the recursion stack of the integrator, namely $O(d)$, where $d$ is the recursion depth. In the actual implementations, $d$ is capped at $30$.

Footnotes and References

Problem 532 source code

C++

#include <cmath>
#include <cstdint>
#include <iomanip>
#include <iostream>

namespace {

long double simpson(long double fa, long double fm, long double fb, long double a, long double b) {
    return (b - a) * (fa + 4 * fm + fb) / 6;
}

long double adaptive_simpson(long double (*f)(long double, long double), long double p, long double a, long double b,
                             long double eps, long double whole, long double fa, long double fm, long double fb,
                             int depth) {
    const long double m = (a + b) / 2;
    const long double l = (a + m) / 2;
    const long double r = (m + b) / 2;
    const long double fl = f(l, p);
    const long double fr = f(r, p);

    const long double left = simpson(fa, fl, fm, a, m);
    const long double right = simpson(fm, fr, fb, m, b);
    const long double delta = left + right - whole;
    if (depth <= 0 || std::fabsl(delta) <= 15 * eps) {
        return left + right + delta / 15;
    }
    return adaptive_simpson(f, p, a, m, eps / 2, left, fa, fl, fm, depth - 1) +
           adaptive_simpson(f, p, m, b, eps / 2, right, fm, fr, fb, depth - 1);
}

// Integrand after substituting u = sin(theta), where theta is colatitude:
// L(n) = (1/s) * ∫_0^{u0} sqrt(1 - (s*u)^2) / (1 - u^2) du,  s = sin(pi/n).
long double integrand_u(long double u, long double s) {
    const long double uu = u * u;
    const long double su = s * u;
    return std::sqrt(1.0L - su * su) / (1.0L - uu);
}

long double length_per_bot(int n) {
    const long double pi = acosl(-1.0L);
    const long double s = std::sin(pi / static_cast<long double>(n));
    const long double u0 = 0.999L;

    const long double a = 0.0L;
    const long double b = u0;
    const long double m = (a + b) / 2;
    const long double fa = integrand_u(a, s);
    const long double fb = integrand_u(b, s);
    const long double fm = integrand_u(m, s);
    const long double whole = simpson(fa, fm, fb, a, b);
    const long double eps = 1e-13L;
    const long double I = adaptive_simpson(&integrand_u, s, a, b, eps, whole, fa, fm, fb, 30);
    return I / s;
}

bool run_checkpoints() {
    const long double L3 = length_per_bot(3);
    // Given: per-bot length rounds to 2.84 for n=3.
    if (std::fabsl(L3 - 2.84L) > 0.01L) {
        std::cerr << std::setprecision(15) << "Checkpoint failed: L(3) got " << L3 << '\n';
        return false;
    }
    return true;
}

}  // namespace

int main() {
    if (!run_checkpoints()) {
        return 1;
    }

    int lo = 3;
    int hi = 3;
    while (length_per_bot(hi) <= 1000.0L) {
        hi *= 2;
    }
    while (lo + 1 < hi) {
        const int mid = lo + (hi - lo) / 2;
        if (length_per_bot(mid) > 1000.0L) {
            hi = mid;
        } else {
            lo = mid;
        }
    }

    const int n = hi;
    const long double per = length_per_bot(n);
    const long double per_prev = length_per_bot(n - 1);
    if (!(per > 1000.0L && per_prev <= 1000.0L)) {
        std::cerr << std::setprecision(15) << "Threshold check failed: n=" << n << " L(n)=" << per
                  << " L(n-1)=" << per_prev << '\n';
        return 2;
    }
    const long double total = static_cast<long double>(n) * per;
    std::cout << std::fixed << std::setprecision(2) << total << '\n';
    return 0;
}

Python

import math
import sys

sys.setrecursionlimit(20000)

def simpson(fa, fm, fb, a, b):
    return (b - a) * (fa + 4 * fm + fb) / 6.0

def integrand_u(u, s):
    uu = u * u
    su = s * u
    return math.sqrt(1.0 - su * su) / (1.0 - uu)

def adaptive_simpson(s, a, b, eps, whole, fa, fm, fb, depth):
    m = (a + b) / 2.0
    l = (a + m) / 2.0
    r = (m + b) / 2.0
    fl = integrand_u(l, s)
    fr = integrand_u(r, s)

    left = simpson(fa, fl, fm, a, m)
    right = simpson(fm, fr, fb, m, b)
    delta = left + right - whole
    
    if depth <= 0 or abs(delta) <= 15.0 * eps:
        return left + right + delta / 15.0
        
    return (adaptive_simpson(s, a, m, eps / 2.0, left, fa, fl, fm, depth - 1) +
            adaptive_simpson(s, m, b, eps / 2.0, right, fm, fr, fb, depth - 1))

def length_per_bot(n):
    s = math.sin(math.pi / n)
    u0 = 0.999
    
    a = 0.0
    b = u0
    m = (a + b) / 2.0
    fa = integrand_u(a, s)
    fb = integrand_u(b, s)
    fm = integrand_u(m, s)
    whole = simpson(fa, fm, fb, a, b)
    eps = 1e-13
    
    I = adaptive_simpson(s, a, b, eps, whole, fa, fm, fb, 30)
    return I / s

def solve():
    lo = 3
    hi = 3
    while length_per_bot(hi) <= 1000.0:
        hi *= 2
        
    while lo + 1 < hi:
        mid = lo + (hi - lo) // 2
        if length_per_bot(mid) > 1000.0:
            hi = mid
        else:
            lo = mid
            
    n = hi
    per = length_per_bot(n)
    total = n * per
    return f"{total:.2f}"

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

Java

public class Euler532 {
    private static double simpson(double fa, double fm, double fb, double a, double b) {
        return (b - a) * (fa + 4 * fm + fb) / 6.0;
    }

    private static double integrandU(double u, double s) {
        double uu = u * u;
        double su = s * u;
        return Math.sqrt(1.0 - su * su) / (1.0 - uu);
    }

    private static double adaptiveSimpson(double s, double a, double b, double eps, double whole,
            double fa, double fm, double fb, int depth) {
        double m = (a + b) / 2.0;
        double l = (a + m) / 2.0;
        double r = (m + b) / 2.0;
        double fl = integrandU(l, s);
        double fr = integrandU(r, s);

        double left = simpson(fa, fl, fm, a, m);
        double right = simpson(fm, fr, fb, m, b);
        double delta = left + right - whole;

        if (depth <= 0 || Math.abs(delta) <= 15.0 * eps) {
            return left + right + delta / 15.0;
        }

        return adaptiveSimpson(s, a, m, eps / 2.0, left, fa, fl, fm, depth - 1) +
                adaptiveSimpson(s, m, b, eps / 2.0, right, fm, fr, fb, depth - 1);
    }

    private static double lengthPerBot(int n) {
        double s = Math.sin(Math.PI / n);
        double u0 = 0.999;

        double a = 0.0;
        double b = u0;
        double m = (a + b) / 2.0;
        double fa = integrandU(a, s);
        double fb = integrandU(b, s);
        double fm = integrandU(m, s);
        double whole = simpson(fa, fm, fb, a, b);
        double eps = 1e-13;

        double I = adaptiveSimpson(s, a, b, eps, whole, fa, fm, fb, 30);
        return I / s;
    }

    public static void main(String[] args) {
        int lo = 3;
        int hi = 3;
        while (lengthPerBot(hi) <= 1000.0) {
            hi *= 2;
        }

        while (lo + 1 < hi) {
            int mid = lo + (hi - lo) / 2;
            if (lengthPerBot(mid) > 1000.0) {
                hi = mid;
            } else {
                lo = mid;
            }
        }

        int n = hi;
        double per = lengthPerBot(n);
        double total = n * per;
        System.out.printf(java.util.Locale.US, "%.2f\n", total);
    }
}