Problem 332: Spherical Triangles

Detailed mathematical approach

Problem Summary

Let \(C(r)\) be the sphere of radius \(r\) centred at the origin and let \(Z(r)\) be the set of its points with integer coordinates,

$$\begin{aligned}Z(r)&=\{(x,y,z)\in\mathbb Z^3:\\&\quad x^2+y^2+z^2=r^2\}.\end{aligned}$$

A spherical triangle is a figure formed on the surface of a sphere by three great circular arcs intersecting pairwise in three vertices. \(T(r)\) is the set of those whose vertices lie in \(Z(r)\), the degenerate ones formed by three points on the same great arc excluded, and \(A(r)\) is the area of the smallest triangle in \(T(r)\). The statement notes that \(A(14)\) is \(3.294040\) rounded to six decimal places, and asks for \(\sum_{r=1}^{50}A(r)\), rounded the same way. ^[1]

The whole solution in one sentence: for each radius list the integer points, build from every triple one determinant \(\Delta\) and one dot-product sum \(D\), both integers, drop the triples with \(\Delta=0\) or \(D\le0\), find the smallest fraction \(\Delta/D\) among the rest by exact integer comparison, and turn it into the area \(2r^2\arctan(\Delta/D)\).

Mathematical Approach

1. The vertices: integer points on a sphere

Finding \(Z(r)\) needs no theory. For each pair \((x,y)\) with \(-r\le x,y\le r\) put \(z^2=r^2-x^2-y^2\); a negative value contributes nothing, and otherwise the pair gives a point exactly when \(z^2\) is a perfect square, with both \(+z\) and \(-z\) (equal only when \(z=0\)). Every point of \(Z(r)\) appears once, at a cost of \((2r+1)^2\) tests.

The six points \((\pm r,0,0)\), \((0,\pm r,0)\), \((0,0,\pm r)\) are always there, so \(Z(r)\) is never empty — as Legendre's three-square theorem also shows, since it excludes only the integers \(4^a(8b+7)\), and \(r=2^am\) with \(m\) odd gives \(r^2=4^am^2\), \(m^2\equiv1\pmod 8\). ^[6]

The size \(n_r=|Z(r)|\) is the number of ways of writing \(r^2\) as a sum of three squares, counting order and signs; it is the coefficient of \(q^{r^2}\) in the theta series of the simple cubic lattice, the generating function of those counts. ^[5] It is very irregular: \(n_1=6\), \(n_3=30\), \(n_7=54\), \(n_{32}=6\), while the largest value for \(r\le 50\) is \(n_{45}=510\).

2. Which triples are triangles

For three points \(u,v,w\in Z(r)\) write

$$\Delta=\bigl|u\cdot(v\times w)\bigr|,$$

the absolute value of the \(3\times3\) determinant with rows \(u,v,w\), the volume of the parallelepiped they span. ^[4] It vanishes precisely when the three vectors are linearly dependent, and since they are nonzero that means they lie in a plane through the origin, whose intersection with the sphere is a great circle. So \(\Delta=0\) says exactly that \(u,v,w\) lie on one great circle, the degenerate case the statement removes; conversely three points on a great circle lie in its plane, which contains the centre.

Because the coordinates are integers, \(\Delta\) is a non-negative integer, so an admissible triple has \(\Delta\ge1\). In particular \(u\) and \(-u\) are never two vertices of a triangle: \(\Delta=0\) for every third point.

Take \(r=1\), where \(Z(1)\) is the six axis points. Of the \(\binom63=20\) triples the \(12\) containing an antipodal pair are degenerate, and the remaining \(8\) are the octant triangles with three right angles.

3. Which of the eight triangles the statement means

The three great circles through the pairs of \(u,v,w\) cut the sphere into eight triangles, one corner of each on each of the three lines through the origin, so "the triangle on \(u,v,w\)" needs a convention. It does not matter which, and the reason also proves the area formula of Section 4.

Since \(\Delta\ne0\) the vectors \(u,v,w\) form a basis, so every point of space is \(\lambda u+\mu v+\nu w\) for unique reals. The plane of \(v\) and \(w\) is \(\{\lambda=0\}\), and likewise for the other two, so the eight regions on the sphere are the eight sign patterns of \((\lambda,\mu,\nu)\), that of \((\varepsilon_1,\varepsilon_2,\varepsilon_3)\) being the part inside the cone spanned by \(\varepsilon_1u,\varepsilon_2v,\varepsilon_3w\): the triangle on those vertices bounded by the three shorter arcs, unique because no two vertices are antipodal (Section 2).

Now \(Z(r)=-Z(r)\) and a sign flip only flips the sign of the determinant, so all eight triples \((\pm u,\pm v,\pm w)\) lie in \(T(r)\) too. The eight triangles of all triples are therefore the same collection as the shorter-arc triangles of all triples, so \(A(r)\) does not depend on the reading. From here on the triangle on \(u,v,w\) is the shorter-arc one,

$$\{\lambda u+\mu v+\nu w:\lambda,\mu,\nu\ge0\}\cap C(r).$$

4. The area of a spherical triangle

Let \(\alpha,\beta,\gamma\) be the angles at \(u,v,w\), between the tangents of the two arcs meeting there. Girard's theorem says the area is \(r^2\) times the spherical excess \(\alpha+\beta+\gamma-\pi\). ^[3]

Here is the classical proof in the notation of Section 3. A lune, the region between two great half-circles from a point to its antipode, of angle \(\alpha\), is the fraction \(\alpha/(2\pi)\) of the sphere, so its area is \(2\alpha r^2\). The two planes through \(u\) cut out a lune of angle \(\alpha\) containing \(u\), which by the description of the regions is the union of our triangle \(T\) and the triangle \(T_u\) on \((-u,v,w)\), so

$$|T|+|T_u|=2\alpha r^2,$$

and likewise for \(v\) and \(w\). The four regions \(T,T_u,T_v,T_w\) contain one of each antipodal pair of the eight, and antipodal regions have equal areas, so \(|T|+|T_u|+|T_v|+|T_w|=2\pi r^2\), half the sphere. Adding the three lune equations gives \(3|T|+(2\pi r^2-|T|)=2(\alpha+\beta+\gamma)r^2\), that is

$$|T|=r^2(\alpha+\beta+\gamma-\pi).$$

The excess is thus the solid angle \(\Omega\) subtended at the centre, the area it cuts out on the unit sphere, ^[2] so the area is \(r^2\Omega\) and the radius enters only through that factor.

5. The solid angle straight from the coordinates

Measuring \(\alpha,\beta,\gamma\) would mean several inverse cosines per triple. The formula of Van Oosterom and Strackee, also proved by Eriksson, avoids them: for three linearly independent vectors with \(\rho\ne0\), ^[2]

$$\tan\frac{\Omega}{2}=\frac{|u\cdot(v\times w)|}{\rho},$$

$$\begin{aligned}\rho&=|u|\,|v|\,|w|+(u\cdot v)\,|w|\\&\quad+(u\cdot w)\,|v|+(v\cdot w)\,|u|.\end{aligned}$$

Every term of \(\rho\) has degree three in the lengths, which makes the formula convenient here: all three are \(r\), so no normalisation is needed and

$$\begin{aligned}\rho=D&:=r^3+r\sigma,\\\sigma&=p+q+t,\end{aligned}$$

where \(p=u\cdot v\), \(q=u\cdot w\) and \(t=v\cdot w\). Both \(\Delta\) and \(D\) are integers, built from the coordinates with multiplications and additions only.

Here is why it holds. Let \(a,b,c\) be the angular side lengths, the central angles opposite \(u,v,w\), so that \(t=r^2\cos a\), \(q=r^2\cos b\), \(p=r^2\cos c\) and \(D=r^3(1+\cos a+\cos b+\cos c)\). The sides \(b\) and \(c\) meet at \(u\) in the angle \(\alpha\), and the excess in terms of two sides and their included angle satisfies, with \(h=\tan\frac b2\tan\frac c2\), ^[3]

$$\tan\frac{\Omega}{2}=\frac{h\sin\alpha}{1+h\cos\alpha}.$$

Write \(\tan\frac b2=\sin b/(1+\cos b)\) and the same for \(c\), then multiply numerator and denominator by \((1+\cos b)(1+\cos c)\). The spherical cosine rule \(\cos a=\cos b\cos c+\sin b\sin c\cos\alpha\) ^[3] collapses the denominator to \(1+\cos a+\cos b+\cos c\), leaving

$$\tan\frac{\Omega}{2}=\frac{\sin b\sin c\sin\alpha}{1+\cos a+\cos b+\cos c}.$$

The numerator is \(\Delta/r^3\): with \(u\) as the vertical axis, the horizontal parts of \(v/r\) and \(w/r\) have lengths \(\sin c\) and \(\sin b\) and meet in the angle \(\alpha\), and \(|u\cdot(v\times w)|/r^3\) is the area of the parallelogram they span. The denominator is \(D/r^3\), so the quotient is \(\Delta/D\).

One detail matters for the implementation. Since \(\Omega/2\) runs over \([0,\pi)\) the tangent alone does not fix the quadrant, and \(D=0\) does occur: at \(r=3\) the triple \((2,2,1)\), \((2,-2,-1)\), \((-2,2,-1)\) has \(\Delta=16\), \(D=0\). The two-argument arctangent settles both, returning a value in \([0,\pi]\) when its first argument is non-negative, ^[7] so

$$\Omega=2\operatorname{atan2}(\Delta,D)$$

is valid for every non-degenerate triple and gives \(\Omega\in[0,2\pi]\).

6. One pair of integers per triple

Both integers are small: by Cauchy–Schwarz and \(|v\times w|\le|v|\,|w|\),

$$1\le\Delta\le|u|\,|v|\,|w|=r^3 .$$

For \(D\), \(p=u\cdot v\le r^2\) with equality only when \(u=v\), since both have length \(r\); the points are distinct, so \(\sigma\lt 3r^2\) and

$$D\lt r^3+r\cdot3r^2=4r^3 .$$

Over all \(85\,584\,072\) triples with \(r\le50\) the largest determinant is \(125\,000=50^3\) and the largest denominator \(495\,400\), just under \(4\cdot50^3\).

7. Comparing candidates exactly

On \(D\gt0\) we have \(\Omega=2\arctan(\Delta/D)\) and \(\arctan\) is strictly increasing, so at a fixed radius the smallest area belongs to the smallest fraction \(\Delta/D\). Two candidates can be ordered without forming a quotient, both denominators being positive:

$$\frac{\Delta_1}{D_1}\lt\frac{\Delta_2}{D_2}\iff\Delta_1D_2\lt\Delta_2D_1 .$$

Both sides are integers below \(125\,000\cdot495\,400\approx6.2\cdot10^{10}\) by Section 6, so the comparison is exact. The smallest fractions are as tiny as \(1/255\,389\), and only the winner of each radius is ever turned into a real number, by a single arctangent at the end.

8. Triples with a non-positive denominator

If \(D\le0\) then \(\operatorname{atan2}(\Delta,D)\ge\pi/2\), so \(\Omega\ge\pi\) and the area is at least \(\pi r^2\). They can be dropped unseen, because a much smaller triangle is always available: the octant on \((r,0,0)\), \((0,r,0)\), \((0,0,r)\) has \(\Delta=r^3\), \(\sigma=0\) and \(D=r^3\), hence

$$\begin{aligned}\Omega&=2\arctan1=\frac{\pi}{2},\\|T|&=\frac{\pi r^2}{2}.\end{aligned}$$

Since \(\pi r^2/2\lt\pi r^2\), the minimum is attained with \(D\gt0\) and the search may be restricted to those. For \(r\le50\) the rule discards \(14\,110\,504\) of the \(84\,614\,064\) triples that survive the determinant test. The same octant gives \(A(r)\le\pi r^2/2\) for every radius.

9. Worked example: the smallest triangle for \(r=3\)

Here \(9=9+0+0=4+4+1\), so \(Z(3)\) is the six axis points together with the \(3\cdot2^3=24\) points with coordinates \(\pm2,\pm2,\pm1\) in some order: \(n_3=30\) and \(\binom{30}3=4060\) triples, \(564\) degenerate and \(592\) of the rest with \(D\le0\). Take

$$\begin{aligned}u&=(-2,-2,1),\\v&=(-2,-2,-1),\\w&=(-2,-1,2).\end{aligned}$$

Then \(p=7\), \(q=8\), \(t=4\), so \(\sigma=19\) and \(D=27+3\cdot19=84\). Also \(v\times w=(-5,6,-2)\) and \(u\cdot(v\times w)=-4\), so \(\Delta=4\) and

$$\begin{aligned}\frac{\Delta}{D}&=\frac{4}{84}=\frac1{21},\\A(3)&=18\arctan\frac1{21},\end{aligned}$$

which is \(0.856496\) to six places. No triple does better, as the exhaustive search confirms; the octant here has area \(9\pi/2\), or \(14.137167\) to six places, larger by a factor of about \(16.51\).

10. The checkpoint \(r=14\)

Here \(n_{14}=54\) and \(\binom{54}3=24\,804\) triples. The minimum is at

$$\begin{aligned}u&=(-12,-6,4),\\v&=(-12,-4,6),\\w&=(-6,-12,-4),\end{aligned}$$

all three of squared norm \(196=14^2\). Their dot products are \(p=192\), \(q=128\), \(t=96\), so \(\sigma=416\) and \(D=2744+14\cdot416=8568\), while the determinant is \(-72\), so \(\Delta=72\). Since \(8568=119\cdot72\),

$$\begin{aligned}A(14)&=392\arctan\frac1{119}\\&=3.294040110677619\ldots,\end{aligned}$$

which rounds to the \(3.294040\) quoted in the statement. This single test exercises the enumeration, the degeneracy test, the area formula and the comparison at once.

11. Doubling the radius, and how small \(A(r)\) can be

Squares are \(0\) or \(1\) modulo \(4\), so three of them add up to a multiple of \(4\) only as \(0+0+0\): every solution of \(x^2+y^2+z^2=(2r)^2\) has all three coordinates even, and therefore

$$Z(2r)=2\,Z(r),$$

so the triples of radius \(2r\) are the doubled triples of radius \(r\). Doubling all three vectors multiplies \(\Delta\) and \(D\) both by \(8\), leaving \(\Delta/D\) and hence \(\Omega\) unchanged, so

$$A(2r)=4A(r).$$

Only the \(25\) odd radii need a search. Iterating from \(Z(1)\), where every triangle is an octant, gives \(A(2^k)=4^k\pi/2\) exactly. The converse is an observation about this range only: for \(r\le50\) enumeration finds exactly six points at \(r=1,2,4,8,16,32\) and nowhere else, which are also the radii attaining the bound \(A(r)=\pi r^2/2\) of Section 8. The largest of the fifty is \(A(32)=512\pi\).

In the other direction the minimising fraction satisfies \(\Delta/D\gt1/(4r^3)\) by Section 6, and it is at most \(1\) because the octant reaches \(1\). On \([0,1]\) the second derivative of \(\arctan\) is \(-2x/(1+x^2)^2\le0\), so its graph lies above the chord from \((0,0)\) to \((1,\pi/4)\): \(\arctan x\ge\pi x/4\) there. Hence

$$\frac{\pi}{8r}\lt A(r)\le\frac{\pi r^2}{2}.$$

The smallest of the fifty values is \(A(41)\), which is \(0.013164\) to six places, against the lower bound \(\pi/328=0.009578\ldots\), so the estimate is close to sharp. A small area needs a small ratio: at \(r=31\) and \(r=41\) the best triple has \(\Delta=1\), the least positive volume an integer parallelepiped can have. That means near-coplanarity, but not a minimal distance: the distance from \(u\) to the plane of \(v,w\) is \(\Delta/|v\times w|\), and the area depends on \(D\) too.

How the Code Works

The vertices of Section 1 come from lattice_points_on_sphere(r), a free function in C++ and one nested inside solve in Python; Java inlines the same loop in minSphericalTriangleArea. It loops over \(x,y\) in \([-r,r]\), forms z2\({}=r^2-x^2-y^2\), skips negative values and rounds its square root to the nearest integer z; the exact test z * z != z2 then decides, so the floating-point root is only a guess. The point \((x,y,z)\) is stored, and \((x,y,-z)\) too when \(z\ne0\), which produces each point once. The C++ version reserves \(600\) entries, safely above the largest set, \(n_{45}=510\).

The search itself is min_spherical_triangle_area in C++, min_area in Python and minSphericalTriangleArea in Java. It first builds a symmetric \(n\times n\) table dots of the pairwise dot products, so each of the three products of a triple later costs one look-up. The triple loop runs over \(i\lt j\lt k\); C++ and Python call abs_det, which expands the determinant of Section 2 along its first row, while Java writes that expansion out in the loop, and a zero value is skipped as degenerate. Otherwise sum_dots is the \(\sigma\) of Section 5, denv\({}=r^3+r\sigma\) is \(D\), and a non-positive denv is skipped — the discard proved in Section 8. Survivors are compared with the running best pair best_det, best_den by the cross-multiplication rule of Section 7; the first candidate of each radius is taken without any multiplication, through a sentinel — the largest signed 64-bit value in C++, Long.MAX_VALUE in Java, None in Python. When the loops end, one call to atan2 on the winning pair and a multiplication by \(r^2\) give the area. The guard returning \(0\) for \(n\lt3\) is defensive only: every sphere carries its six axis points.

solve_sum adds the fifty values in order of radius, printed with six decimals. The C++ program also has run_checkpoints, which recomputes \(A(14)\) and exits with a nonzero status unless it agrees with \(3.294040\) to within \(5\cdot10^{-7}\); --skip-checkpoints turns it off. The ports carry no such check; their solve functions format the sum directly.

The three versions run the same search and differ in numeric types. C++ accumulates areas in long double and widens the comparison to __int128; Java uses double and compares in isLess with BigInteger; ^[8] Python needs neither, its integers being unbounded. By Section 7 the widening is a safety margin, not a necessity: the products never exceed \(6.2\cdot10^{10}\), inside the signed 64-bit range. The dot table is a long[n][n] in Java and a vector of vectors in C++, about \(2\) MB at \(n=510\); Python's lists of integer objects use several times more.

Complexity and Verification

Per radius the enumeration costs \((2r+1)^2\le10\,201\) tests and the dot table \(O(n_r^2)\) time and memory; the triple loop dominates with \(\binom{n_r}3\) iterations of constant work, so the method is \(O(r^2+n_r^3)\) per radius. The summation pass visits \(\sum_{r=1}^{50}\binom{n_r}3=85\,584\,072\) triples; the determinant test rejects \(970\,008\) and the denominator test a further \(14\,110\,504\), leaving \(70\,503\,560\) candidates and, one sentinel per radius aside, \(70\,503\,510\) comparisons. Radius \(45\) alone contributes \(21\,978\,620\) visits, and the C++ checkpoint adds a pass over the \(24\,804\) triples of radius \(14\). Measured on one machine, C++ built with -O2 takes \(0.17\) s, the Java port \(1.8\) s and CPython 3.11 \(23\) s, all three printing the same six decimals.

Growth in \(n_r\), not in \(r\), is the real limit: the loop is cubic in \(n_r\), already \(510\) at \(r=45\). The doubling rule of Section 11 helps far less than the halved number of radii suggests, because \(n_{2r}=n_r\) makes the even spheres copies of smaller ones: they are only \(6\,044\,996\) of the \(85\,584\,072\) visits, about \(7\) percent. A real saving comes from \(D\lt4r^3\): once the best fraction is \(b/c\), a candidate with \(\Delta c\ge4r^3b\) has \(\Delta/D\gt\Delta/(4r^3)\ge b/c\) and can go as soon as its determinant is known. Neither shortcut is used by the published programs.

Beyond the \(A(14)\) checkpoint built into the C++ program, these checks were run for this article.

• An independent re-implementation reproduced \(A(r)\) for every \(r\le50\) — also when minimising the floating-point area over all non-degenerate triples, with no \(D\le0\) filter — and its total agrees with the printed output of all three published programs.
• Each winning area, recomputed from the three angles by Girard's theorem at \(60\)-digit precision, agrees with \(2r^2\arctan(\Delta/D)\) to over \(50\) digits.
• Every step of the chain in Section 5, up to \(\Omega=2\operatorname{atan2}(\Delta,D)\), holds to \(10^{-57}\) on \(300\) random triangles, \(51\) of them with \(D\lt0\).
• The eight sign variants of a triple form four antipodal pairs of equal area; minimising over one representative of each pair, not just the shorter-arc triangle, returns the same \(A(r)\) for every \(r\le15\), as Section 3 predicts.
• \(Z(2r)=2Z(r)\) was confirmed for \(r\le25\), the reduced fractions \(\Delta/D\) agree for \(r\) and \(2r\), and the values satisfy \(A(2r)=4A(r)\). The counts \(n_r\) match the representation numbers of \(r^2\) for every \(r\le50\).

Footnotes and References

The references below supply the background facts used above; the equivalence of the eight readings, the integer form of the area, the discard of the non-positive denominators, the doubling rule and the bounds on \(A(r)\) are derived here.

1. Project Euler 332 — Spherical Triangles. The official statement defines the sphere C(r), its integer points Z(r), the set T(r) of spherical triangles with vertices in Z(r) — those on one great arc excluded — and the smallest such area A(r). It records A(14) as 3.294040 to six places and asks for the sum over r = 1 to 50, rounded the same way.
2. Wikipedia — Solid angle. A solid angle is measured in steradians as the area it covers on the unit sphere about the apex. The article states the vector formula tan(Ω/2) = |a·(b×c)| / (abc + (a·b)c + (a·c)b + (b·c)a) for a triangular cone, citing F. Eriksson, Mathematics Magazine 63 (1990), and A. Van Oosterom and J. Strackee, IEEE Transactions on Biomedical Engineering (1983), which Section 5 specialises to vectors of equal length.
3. Wikipedia — Spherical trigonometry. The spherical excess is E = A + B + C − π; the area is E on the unit sphere and R²E on a sphere of radius R. The cosine rule reads cos a = cos b cos c + sin b sin c cos A, and the section on area and spherical excess gives tan(E/2) = tan(a/2) tan(b/2) sin C / (1 + tan(a/2) tan(b/2) cos C) for two sides and their included angle; Section 5 combines these with the scalar triple product.
4. Wikipedia — Triple product. The scalar triple product a·(b×c) is the determinant of the 3×3 matrix with the three vectors as rows and the signed volume of the parallelepiped they define; a zero value means they are coplanar, the parallelepiped being flat. Section 2 uses this in both directions.
5. OEIS A005875. The theta series of the simple cubic lattice, whose coefficients count the ways of writing n as a sum of three squares with order and signs: 1, 6, 12, 8, 6, 24, 24, 0, 12, 30, … The entry writes the count r₃(n); the size of Z(r) in Section 1 is its value at r².
6. Wikipedia — Legendre's three-square theorem. A natural number is a sum of three integer squares exactly when it is not of the form 4ᵃ(8b+7), the necessity following because every square is 0, 1 or 4 modulo 8. Section 1 uses it to confirm that r² is always representable.
7. cppreference — std::atan2. The function computes the arc tangent of y/x using the signs of both arguments to pick the quadrant and returns a value in [−π, π]; its example shows atan2(1, −1) = 3π/4. This lets Section 5 write 2·atan2(Δ, D) for solid angles up to 2π.
8. GCC manual — 128-bit integers documents __int128 as an extension supported where an integer mode is wide enough for 128 bits, the type the C++ comparison widens to. Java — BigInteger provides immutable arbitrary-precision integers with multiply and compareTo, used by the Java port for the same comparison. Section 7 bounds those products below 6.2·10¹⁰, so both are precautions.