Problem 332: Spherical Triangles

Project Euler Problem 332 Solution

EulerSolve provides an optimized solution for Project Euler Problem 332, Spherical Triangles, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For each integer radius $1 \le r \le 50$, consider the lattice points on the sphere $x^2+y^2+z^2=r^2$. Among all non-degenerate triples of such points, let $A(r)$ be the minimum spherical triangle area. The task is to compute $$\sum_{r=1}^{50} A(r).$$ Mathematical Approach Fix a radius $r$. Every admissible vertex lies in $$\mathcal P_r=\{(x,y,z)\in\mathbb Z^3:x^2+y^2+z^2=r^2\}.$$ The solution has three parts: generate $\mathcal P_r$, evaluate the spherical area of each triple, and keep the minimum positive area. 1) Enumerating Lattice Points on the Sphere The code loops over $x,y\in[-r,r]$, computes $$z^2=r^2-x^2-y^2,$$ and keeps the pair $(x,y)$ only when $z^2$ is a non-negative perfect square. Then both $z$ and $-z$ are inserted when $z\neq 0$. This generates exactly the integer points on the sphere. For $r\le 50$, the resulting set is small enough that a direct triple search is still practical. 2) Area from a Solid-Angle Formula Take three lattice points $a,b,c\in\mathcal P_r$. After normalization, $$u=\frac{a}{r},\qquad v=\frac{b}{r},\qquad w=\frac{c}{r},$$ lie on the unit sphere. The spherical triangle area on the sphere of radius $r$ is $r^2\Omega$, where $\Omega$ is the solid angle seen from the origin....

Detailed mathematical approach

Problem Summary

For each integer radius $1 \le r \le 50$, consider the lattice points on the sphere $x^2+y^2+z^2=r^2$. Among all non-degenerate triples of such points, let $A(r)$ be the minimum spherical triangle area. The task is to compute

$$\sum_{r=1}^{50} A(r).$$

Mathematical Approach

Fix a radius $r$. Every admissible vertex lies in

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

The solution has three parts: generate $\mathcal P_r$, evaluate the spherical area of each triple, and keep the minimum positive area.

1) Enumerating Lattice Points on the Sphere

The code loops over $x,y\in[-r,r]$, computes

$$z^2=r^2-x^2-y^2,$$

and keeps the pair $(x,y)$ only when $z^2$ is a non-negative perfect square. Then both $z$ and $-z$ are inserted when $z\neq 0$. This generates exactly the integer points on the sphere.

For $r\le 50$, the resulting set is small enough that a direct triple search is still practical.

2) Area from a Solid-Angle Formula

Take three lattice points $a,b,c\in\mathcal P_r$. After normalization,

$$u=\frac{a}{r},\qquad v=\frac{b}{r},\qquad w=\frac{c}{r},$$

lie on the unit sphere. The spherical triangle area on the sphere of radius $r$ is $r^2\Omega$, where $\Omega$ is the solid angle seen from the origin.

For unit vectors, a standard identity is

$$\Omega=2\operatorname{atan2}\!\left(|u\cdot(v\times w)|,\ 1+u\cdot v+u\cdot w+v\cdot w\right).$$

Scaling back to the radius-$r$ vectors gives

$$\Delta=|\det(a,b,c)|=r^3|u\cdot(v\times w)|,$$

$$D=r^3+r(a\cdot b+a\cdot c+b\cdot c)=r^3(1+u\cdot v+u\cdot w+v\cdot w).$$

Therefore the exact formula implemented in the code is

$$\Omega=2\operatorname{atan2}(\Delta,D),\qquad \operatorname{Area}=r^2\Omega.$$

When $D>0$, this is the same as $2\arctan(\Delta/D)$, but the `atan2` form is numerically safer because it keeps numerator and denominator separate.

3) Why the Determinant Appears

The quantity $|\det(a,b,c)|$ is the volume of the parallelepiped spanned by the three radius vectors. It is zero exactly when the vectors are coplanar with the origin, which means the spherical triangle is degenerate. So the code simply skips all triples with

$$\Delta=0.$$

4) Why Minimizing Area Reduces to Minimizing $\Delta/D$

For the candidates kept by the code we have $\Delta>0$ and $D>0$. On that domain, $x\mapsto 2\arctan x$ is strictly increasing, so minimizing the area is equivalent to minimizing

$$\frac{\Delta}{D}.$$

This matters algorithmically: instead of computing floating-point areas for every triple, the program compares two candidates $\Delta_1/D_1$ and $\Delta_2/D_2$ using exact cross multiplication,

$$\Delta_1D_2\lt \Delta_2D_1.$$

The C++ version uses `__int128`, the Java version uses `BigInteger`, and Python relies on arbitrary-precision integers.

The condition $D\le 0$ can also be discarded safely. Such triples satisfy $\Omega\ge\pi$, hence their area is at least $\pi r^2$. But the axis points $(r,0,0)$, $(0,r,0)$, $(0,0,r)$ are always present and give area $\pi r^2/2$, so a minimum can never come from $D\le 0$.

5) Worked Example: $r=1$

For $r=1$, the lattice points are the six axis points $\pm(1,0,0)$, $\pm(0,1,0)$, $\pm(0,0,1)$. Choose

$$a=(1,0,0),\qquad b=(0,1,0),\qquad c=(0,0,1).$$

Then

$$\Delta=|\det(a,b,c)|=1,$$

and every pairwise dot product is $0$, so $D=1$. Therefore

$$\Omega=2\operatorname{atan2}(1,1)=\frac{\pi}{2},\qquad A(1)=1^2\cdot\frac{\pi}{2}=\frac{\pi}{2}.$$

This is exactly the spherical octant triangle with three right angles, and it is a clean sanity check for the formula.

How the Code Works

For each radius $r$, the implementation first builds the list of lattice points on the sphere. It then precomputes all pairwise dot products in an $n_r\times n_r$ table, where $n_r=|\mathcal P_r|$. This makes the denominator

$$D=r^3+r(a\cdot b+a\cdot c+b\cdot c)$$

available in $O(1)$ time for each triple.

Next it scans all $i<j<k$, skips degenerate triples, ignores $D\le 0$, and keeps the best fraction $\Delta/D$ using exact integer comparison. Only after the best triple is known does it evaluate one final `atan2` and multiply by $r^2$. The C++ code also checks the published checkpoint $A(14)\approx 3.294040$ before summing all radii.

Complexity Analysis

Let $n_r=|\mathcal P_r|$. Point generation takes $O(r^2)$ trial pairs $(x,y)$, dot-product precomputation takes $O(n_r^2)$ time and memory, and the dominant triple scan takes $O(n_r^3)$ time.

Because the problem stops at $r=50$, every $n_r$ stays small enough that this brute-force geometric search is completely feasible. The exact-integer comparison also keeps the result stable across all three language implementations.

Problem 332 source code

C++

#include <algorithm>
#include <cmath>
#include <cstdint>
#include <iomanip>
#include <iostream>
#include <limits>
#include <string>
#include <vector>

namespace {

using i64 = std::int64_t;
using i128 = __int128_t;

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

std::vector<Point> lattice_points_on_sphere(const int r) {
    const int rr = r * r;
    std::vector<Point> points;
    points.reserve(600);

    for (int x = -r; x <= r; ++x) {
        const int x2 = x * x;
        for (int y = -r; y <= r; ++y) {
            const int y2 = y * y;
            const int z2 = rr - x2 - y2;
            if (z2 < 0) {
                continue;
            }
            const int z = static_cast<int>(std::llround(std::sqrt(static_cast<long double>(z2))));
            if (z * z != z2) {
                continue;
            }
            points.push_back({x, y, z});
            if (z != 0) {
                points.push_back({x, y, -z});
            }
        }
    }

    return points;
}

i64 dot(const Point& a, const Point& b) {
    return static_cast<i64>(a.x) * b.x + static_cast<i64>(a.y) * b.y + static_cast<i64>(a.z) * b.z;
}

i64 abs_det(const Point& a, const Point& b, const Point& c) {
    const i64 det = static_cast<i64>(a.x) * (static_cast<i64>(b.y) * c.z - static_cast<i64>(b.z) * c.y) -
                    static_cast<i64>(a.y) * (static_cast<i64>(b.x) * c.z - static_cast<i64>(b.z) * c.x) +
                    static_cast<i64>(a.z) * (static_cast<i64>(b.x) * c.y - static_cast<i64>(b.y) * c.x);
    return det >= 0 ? det : -det;
}

long double min_spherical_triangle_area(const int r) {
    const std::vector<Point> points = lattice_points_on_sphere(r);
    const int n = static_cast<int>(points.size());
    if (n < 3) {
        return 0.0L;
    }

    std::vector<std::vector<i64>> dots(static_cast<std::size_t>(n), std::vector<i64>(static_cast<std::size_t>(n), 0));
    for (int i = 0; i < n; ++i) {
        for (int j = i + 1; j < n; ++j) {
            const i64 v = dot(points[static_cast<std::size_t>(i)], points[static_cast<std::size_t>(j)]);
            dots[static_cast<std::size_t>(i)][static_cast<std::size_t>(j)] = v;
            dots[static_cast<std::size_t>(j)][static_cast<std::size_t>(i)] = v;
        }
    }

    i64 best_det = std::numeric_limits<i64>::max();
    i64 best_den = 1;

    const i64 r2 = static_cast<i64>(r) * r;
    const i64 r3 = r2 * static_cast<i64>(r);

    for (int i = 0; i < n; ++i) {
        for (int j = i + 1; j < n; ++j) {
            for (int k = j + 1; k < n; ++k) {
                const i64 detv = abs_det(points[static_cast<std::size_t>(i)], points[static_cast<std::size_t>(j)],
                                         points[static_cast<std::size_t>(k)]);
                if (detv == 0) {
                    continue;
                }

                const i64 sum_dots = dots[static_cast<std::size_t>(i)][static_cast<std::size_t>(j)] +
                                     dots[static_cast<std::size_t>(i)][static_cast<std::size_t>(k)] +
                                     dots[static_cast<std::size_t>(j)][static_cast<std::size_t>(k)];
                const i64 denv = r3 + static_cast<i64>(r) * sum_dots;

                if (denv <= 0) {
                    continue;
                }

                // Minimize det/den for positive numerator/denominator.
                if (best_det == std::numeric_limits<i64>::max() ||
                    static_cast<i128>(detv) * static_cast<i128>(best_den) <
                        static_cast<i128>(best_det) * static_cast<i128>(denv)) {
                    best_det = detv;
                    best_den = denv;
                }
            }
        }
    }

    if (best_det == std::numeric_limits<i64>::max()) {
        return 0.0L;
    }

    const long double omega = 2.0L * std::atan2(static_cast<long double>(best_det), static_cast<long double>(best_den));
    return static_cast<long double>(r2) * omega;
}

long double solve_sum() {
    long double total = 0.0L;
    for (int r = 1; r <= 50; ++r) {
        total += min_spherical_triangle_area(r);
    }
    return total;
}

bool run_checkpoints() {
    const long double a14 = min_spherical_triangle_area(14);
    const long double expected = 3.294040L;
    if (std::fabsl(a14 - expected) > 0.0000005L) {
        std::cerr << "Checkpoint failed: A(14) mismatch, got=" << std::fixed << std::setprecision(9) << a14 << '\n';
        return false;
    }
    return true;
}

}  // namespace

int main(int argc, char** argv) {
    bool skip_checkpoints = false;
    for (int i = 1; i < argc; ++i) {
        const 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_sum();
    std::cout << std::fixed << std::setprecision(6) << answer << '\n';
    return 0;
}

Python

import math

def solve():
    def lattice_points_on_sphere(r):
        rr = r * r
        points = []
        for x in range(-r, r+1):
            x2 = x * x
            for y in range(-r, r+1):
                y2 = y * y
                z2 = rr - x2 - y2
                if z2 < 0:
                    continue
                z = round(math.sqrt(z2))
                if z * z != z2:
                    continue
                points.append((x, y, z))
                if z != 0:
                    points.append((x, y, -z))
        return points

    def dot(a, b):
        return a[0]*b[0] + a[1]*b[1] + a[2]*b[2]

    def abs_det(a, b, c):
        det = (a[0]*(b[1]*c[2]-b[2]*c[1]) -
               a[1]*(b[0]*c[2]-b[2]*c[0]) +
               a[2]*(b[0]*c[1]-b[1]*c[0]))
        return abs(det)

    def min_area(r):
        pts = lattice_points_on_sphere(r)
        n = len(pts)
        if n < 3:
            return 0.0
        dots = [[0]*n for _ in range(n)]
        for i in range(n):
            for j in range(i+1, n):
                v = dot(pts[i], pts[j])
                dots[i][j] = v
                dots[j][i] = v

        best_det = None
        best_den = 1
        r2 = r * r

        for i in range(n):
            for j in range(i+1, n):
                for k in range(j+1, n):
                    detv = abs_det(pts[i], pts[j], pts[k])
                    if detv == 0:
                        continue
                    sd = dots[i][j] + dots[i][k] + dots[j][k]
                    denv = r * r2 + r * sd
                    if denv <= 0:
                        continue
                    if best_det is None or detv * best_den < best_det * denv:
                        best_det = detv
                        best_den = denv

        if best_det is None:
            return 0.0
        omega = 2.0 * math.atan2(best_det, best_den)
        return r2 * omega

    total = 0.0
    for r in range(1, 51):
        total += min_area(r)
    return f"{total:.6f}"

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

Java

import java.util.*;

public class Euler332 {

    static class Point {
        int x, y, z;

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

    static double minSphericalTriangleArea(int r) {
        int rr = r * r;
        List<Point> points = new ArrayList<>();

        for (int x = -r; x <= r; x++) {
            int x2 = x * x;
            for (int y = -r; y <= r; y++) {
                int y2 = y * y;
                int z2 = rr - x2 - y2;
                if (z2 < 0)
                    continue;
                int z = (int) Math.round(Math.sqrt(z2));
                if (z * z != z2)
                    continue;
                points.add(new Point(x, y, z));
                if (z != 0) {
                    points.add(new Point(x, y, -z));
                }
            }
        }

        int n = points.size();
        if (n < 3)
            return 0.0;

        long[][] dots = new long[n][n];
        for (int i = 0; i < n; i++) {
            for (int j = i + 1; j < n; j++) {
                long v = (long) points.get(i).x * points.get(j).x +
                        (long) points.get(i).y * points.get(j).y +
                        (long) points.get(i).z * points.get(j).z;
                dots[i][j] = v;
                dots[j][i] = v;
            }
        }

        long bestDet = Long.MAX_VALUE;
        long bestDen = 1;

        long r2 = (long) r * r;
        long r3 = r2 * r;

        for (int i = 0; i < n; i++) {
            Point pi = points.get(i);
            for (int j = i + 1; j < n; j++) {
                Point pj = points.get(j);
                for (int k = j + 1; k < n; k++) {
                    Point pk = points.get(k);

                    long det = (long) pi.x * ((long) pj.y * pk.z - (long) pj.z * pk.y) -
                            (long) pi.y * ((long) pj.x * pk.z - (long) pj.z * pk.x) +
                            (long) pi.z * ((long) pj.x * pk.y - (long) pj.y * pk.x);

                    long detv = Math.abs(det);
                    if (detv == 0)
                        continue;

                    long sumDots = dots[i][j] + dots[i][k] + dots[j][k];
                    long denv = r3 + r * sumDots;

                    if (denv <= 0)
                        continue;

                    // detv/denv < bestDet/bestDen => detv*bestDen < bestDet*denv
                    // need using BigInteger or float point approximations, or simply rely on
                    // Math.multiplyHigh, but since values can be up to ~50^6 (1.5e10), detv *
                    // bestDen can be (~1.5e10)^2 = 2.25e20 which slightly exceeds Long.MAX_VALUE
                    // (9e18).
                    // We can use double division for quick check, then BigInteger if close, or just
                    // BigInteger.
                    if (bestDet == Long.MAX_VALUE || isLess(detv, denv, bestDet, bestDen)) {
                        bestDet = detv;
                        bestDen = denv;
                    }
                }
            }
        }

        if (bestDet == Long.MAX_VALUE)
            return 0.0;

        double omega = 2.0 * Math.atan2((double) bestDet, (double) bestDen);
        return r2 * omega;
    }

    static boolean isLess(long aNum, long aDen, long bNum, long bDen) {
        // aNum/aDen < bNum/bDen -> aNum*bDen < bNum*aDen
        java.math.BigInteger left = java.math.BigInteger.valueOf(aNum).multiply(java.math.BigInteger.valueOf(bDen));
        java.math.BigInteger right = java.math.BigInteger.valueOf(bNum).multiply(java.math.BigInteger.valueOf(aDen));
        return left.compareTo(right) < 0;
    }

    public static String solve() {
        double total = 0.0;
        for (int r = 1; r <= 50; r++) {
            total += minSphericalTriangleArea(r);
        }
        return String.format(Locale.US, "%.6f", total);
    }

    public static void main(String[] args) {
        System.out.println(solve());
    }
}