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.

Further Reading

  1. Problem page: https://projecteuler.net/problem=332
  2. Solid angle: https://en.wikipedia.org/wiki/Solid_angle
  3. Spherical triangle and spherical excess: https://en.wikipedia.org/wiki/Spherical_triangle
  4. Scalar triple product: https://en.wikipedia.org/wiki/Triple_product

Problemzusammenfassung

Für jeden ganzzahligen Radius \(1 \le r \le 50\) werden die Gitterpunkte auf der Sphäre \(x^2+y^2+z^2=r^2\) betrachtet. Unter allen nicht-degenerierten Punkttripeln sei \(A(r)\) die kleinste sphärische Dreiecksfläche. Gesucht ist

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

Mathematischer Ansatz

Für festes \(r\) liegen alle zulässigen Ecken in

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

Die Lösung besteht aus drei Schritten: \(\mathcal P_r\) erzeugen, für jedes Tripel die sphärische Fläche berechnen und die kleinste positive Fläche behalten.

1) Gitterpunkte auf der Sphäre erzeugen

Der Code durchläuft \(x,y\in[-r,r]\), berechnet

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

und akzeptiert ein Paar \((x,y)\) nur dann, wenn \(z^2\) nichtnegativ und ein perfektes Quadrat ist. Falls \(z\neq 0\), werden sowohl \(z\) als auch \(-z\) eingefügt. So entstehen genau die ganzzahligen Punkte der Sphäre.

Für \(r\le 50\) bleibt diese Punktmenge klein genug, sodass ein direkter Tripel-Scan praktikabel ist.

2) Fläche über eine Raumwinkelformel

Seien \(a,b,c\in\mathcal P_r\). Nach der Normierung

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

liegen \(u,v,w\) auf der Einheitskugel. Die sphärische Dreiecksfläche auf einer Kugel mit Radius \(r\) ist \(r^2\Omega\), wobei \(\Omega\) der vom Ursprung aus gesehene Raumwinkel ist.

Für Einheitsvektoren gilt die Standardformel

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

Zurückskaliert auf die Radius-\(r\)-Vektoren erhält man

$$\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).$$

Damit lautet die exakt im Code verwendete Formel

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

Für \(D>0\) ist das gleich \(2\arctan(\Delta/D)\), aber die `atan2`-Schreibweise ist numerisch robuster, weil Zähler und Nenner getrennt bleiben.

3) Warum die Determinante erscheint

Die Größe \(|\det(a,b,c)|\) ist das Volumen des von den drei Radiusvektoren aufgespannten Parallelepipeds. Sie verschwindet genau dann, wenn die drei Vektoren zusammen mit dem Ursprung koplanar sind, also wenn das sphärische Dreieck degeneriert. Deshalb verwirft der Code alle Tripel mit

$$\Delta=0.$$

4) Warum die Minimierung auf \(\Delta/D\) reduziert

Für die im Code berücksichtigten Kandidaten gilt \(\Delta>0\) und \(D>0\). Auf diesem Bereich ist \(x\mapsto 2\arctan x\) streng monoton wachsend. Also ist die Minimierung der Fläche äquivalent zur Minimierung von

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

Das ist rechnerisch wichtig: Statt für jedes Tripel Fließkommaflächen zu berechnen, vergleicht das Programm zwei Kandidaten \(\Delta_1/D_1\) und \(\Delta_2/D_2\) exakt über Kreuzmultiplikation,

$$\Delta_1D_2<\Delta_2D_1.$$

Die C++-Version verwendet dafür `__int128`, die Java-Version `BigInteger`, und Python nutzt eingebaute Ganzzahlen beliebiger Länge.

Auch \(D\le 0\) kann sicher verworfen werden. Dann gilt \(\Omega\ge\pi\), also ist die Fläche mindestens \(\pi r^2\). Die Achsenpunkte \((r,0,0)\), \((0,r,0)\), \((0,0,r)\) existieren jedoch immer und liefern die Fläche \(\pi r^2/2\). Ein Minimum kann daher nie aus einem Tripel mit \(D\le 0\) stammen.

5) Beispiel: \(r=1\)

Für \(r=1\) bestehen die Gitterpunkte genau aus den sechs Achsenpunkten \(\pm(1,0,0)\), \(\pm(0,1,0)\), \(\pm(0,0,1)\). Wähle

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

Dann ist

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

und alle paarweisen Skalarprodukte sind \(0\), also \(D=1\). Damit folgt

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

Das ist genau das sphärische Oktanten-Dreieck mit drei rechten Winkeln und ein sehr guter Plausibilitätstest für die Formel.

Funktionsweise des Codes

Für jedes \(r\) baut die Implementierung zunächst die Liste der Gitterpunkte auf der Sphäre auf. Danach werden alle paarweisen Skalarprodukte in einer \(n_r\times n_r\)-Tabelle vorkalkuliert, wobei \(n_r=|\mathcal P_r|\). Damit ist der Nenner

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

für jedes Tripel in \(O(1)\) verfügbar.

Anschließend durchläuft das Programm alle \(i<j<k\), überspringt degenerierte Tripel, ignoriert \(D\le 0\) und speichert den besten Bruch \(\Delta/D\) per exaktem Ganzzahlvergleich. Erst ganz am Ende wird einmal `atan2` ausgewertet und mit \(r^2\) multipliziert. Der C++-Code prüft zusätzlich den veröffentlichten Kontrollwert \(A(14)\approx 3.294040\), bevor über alle Radien summiert wird.

Komplexitätsanalyse

Sei \(n_r=|\mathcal P_r|\). Die Punkterzeugung benötigt \(O(r^2)\) Testpaare \((x,y)\), die Vorberechnung der Skalarprodukte benötigt \(O(n_r^2)\) Zeit und Speicher, und der dominante Tripel-Scan benötigt \(O(n_r^3)\) Zeit.

Da das Problem bei \(r=50\) endet, bleibt jedes \(n_r\) klein genug, sodass diese brute-forceartige geometrische Suche völlig ausreicht. Die exakten Ganzzahlvergleiche halten das Ergebnis außerdem in allen drei Programmiersprachen stabil.

Weiterführende Hinweise

  1. Problemseite: https://projecteuler.net/problem=332
  2. Raumwinkel: https://en.wikipedia.org/wiki/Solid_angle
  3. Sphärisches Dreieck und sphärischer Exzess: https://en.wikipedia.org/wiki/Spherical_triangle
  4. Spatprodukt / gemischtes Produkt: https://en.wikipedia.org/wiki/Triple_product

Problem Özeti

Her \(1 \le r \le 50\) tam sayı yarıçapı için \(x^2+y^2+z^2=r^2\) küresi üzerindeki tam sayı kafes noktaları inceleniyor. Bu noktalardan oluşan tüm dejenere olmayan üçlüler arasında en küçük küresel üçgen alanına \(A(r)\) diyelim. İstenen değer

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

toplamıdır.

Matematiksel Yaklaşım

Sabit bir \(r\) için tüm uygun köşeler şu kümede yer alır:

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

Çözüm üç parçadan oluşur: \(\mathcal P_r\) kümesini üretmek, her üçlü için küresel alanı hesaplamak ve en küçük pozitif alanı saklamak.

1) Küre Üzerindeki Kafes Noktalarını Üretmek

Kod \(x,y\in[-r,r]\) aralığını dolaşır ve

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

değerini hesaplar. \(z^2\) negatif değilse ve tam kare ise nokta geçerlidir. Ayrıca \(z\neq 0\) durumunda hem \(z\) hem de \(-z\) eklenir. Böylece küre üzerindeki tüm tam sayı noktaları eksiksiz üretilmiş olur.

\(r\le 50\) için oluşan nokta kümesi yeterince küçük kaldığından tüm üçlüleri doğrudan taramak mümkündür.

2) Katı Açıdan Küresel Alan Formülü

\(a,b,c\in\mathcal P_r\) olsun. Normalize edersek

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

vektörleri birim küre üzerindedir. Yarıçapı \(r\) olan küredeki küresel üçgen alanı \(r^2\Omega\)'dır; burada \(\Omega\) orijinden görülen katı açıdır.

Birim vektörler için standart özdeşlik

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

şeklindedir. Bunu yarıçapı \(r\) olan vektörlere geri ölçekleyince

$$\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)$$

elde edilir. Dolayısıyla kodun kullandığı tam formül

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

olur. \(D>0\) iken bu ifade \(2\arctan(\Delta/D)\) ile aynıdır; ancak `atan2` biçimi pay ve paydayı ayrı tuttuğu için sayısal olarak daha güvenlidir.

3) Determinant Neden Ortaya Çıkıyor?

\(|\det(a,b,c)|\), üç yarıçap vektörünün oluşturduğu paralelyüzlünün hacmidir. Bu değer yalnızca vektörler orijinle birlikte aynı düzleme düştüğünde sıfır olur; bu da küresel üçgenin dejenere olduğu durumdur. Bu nedenle kod

$$\Delta=0$$

olan tüm üçlüleri doğrudan eler.

4) Neden Minimizasyon \(\Delta/D\)'ye İniyor?

Kodun tuttuğu adaylarda \(\Delta>0\) ve \(D>0\) vardır. Bu bölgede \(x\mapsto 2\arctan x\) fonksiyonu sıkı biçimde artandır. O halde alanı minimize etmek,

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

oranını minimize etmekle eşdeğerdir.

Bu, algoritma açısından çok önemlidir: Program her üçlü için kayan noktalı alan hesaplamak yerine iki adayı \(\Delta_1/D_1\) ve \(\Delta_2/D_2\) tam sayı çapraz çarpımıyla karşılaştırır:

$$\Delta_1D_2<\Delta_2D_1.$$

C++ sürümü bunun için `__int128`, Java sürümü `BigInteger`, Python ise doğal büyük tamsayı aritmetiği kullanır.

\(D\le 0\) durumu da güvenle atılabilir. Böyle adaylarda \(\Omega\ge\pi\) olur; dolayısıyla alan en az \(\pi r^2\)'dir. Oysa \((r,0,0)\), \((0,r,0)\), \((0,0,r)\) eksen noktaları her zaman vardır ve \(\pi r^2/2\) alanlı bir üçgen oluşturur. Bu yüzden minimum alan hiçbir zaman \(D\le 0\) olan bir adaydan gelmez.

5) Örnek: \(r=1\)

\(r=1\) için kafes noktaları tam olarak altı eksen noktasıdır: \(\pm(1,0,0)\), \(\pm(0,1,0)\), \(\pm(0,0,1)\). Şu üçlüyü seçelim:

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

Bu durumda

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

ve tüm ikili skaler çarpımlar \(0\) olduğu için \(D=1\) olur. Sonuç olarak

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

Bu üçgen, üç dik açıya sahip küresel oktant üçgenidir ve formül için temiz bir doğrulama örneği verir.

Kodun Çalışma Mantığı

Her \(r\) için uygulama önce küre üzerindeki tüm kafes noktalarını çıkarır. Sonra \(n_r=|\mathcal P_r|\) olmak üzere tüm ikili skaler çarpımları \(n_r\times n_r\) boyutlu bir tabloda önceden hesaplar. Böylece

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

ifadesi her üçlü için \(O(1)\) sürede elde edilir.

Daha sonra tüm \(i<j<k\) üçlüleri taranır; dejenere olanlar atılır, \(D\le 0\) olanlar göz ardı edilir ve en iyi \(\Delta/D\) kesri tam sayı karşılaştırmasıyla saklanır. Sadece kazanan aday bulunduktan sonra bir kez `atan2` çağrılır ve sonuç \(r^2\) ile çarpılır. C++ sürümü ayrıca tüm yarıçapları toplamadan önce yayımlanmış kontrol değeri \(A(14)\approx 3.294040\)'ı sınar.

Karmaşıklık Analizi

\(n_r=|\mathcal P_r|\) olsun. Nokta üretimi \(O(r^2)\) adet \((x,y)\) denemesi gerektirir. Skaler çarpım önhesabı \(O(n_r^2)\) zaman ve bellek kullanır. Baskın adım olan üçlü taraması ise \(O(n_r^3)\) zamandır.

Problem yalnızca \(r=50\)'ye kadar gittiği için her \(n_r\) yeterince küçük kalır ve bu doğrudan geometrik arama rahatlıkla çalışır. Tam sayı tabanlı karşılaştırmalar da üç dilde aynı sonucu güvenilir biçimde korur.

İleri Okuma

  1. Problem sayfası: https://projecteuler.net/problem=332
  2. Katı açı: https://en.wikipedia.org/wiki/Solid_angle
  3. Küresel üçgen ve küresel fazlalık: https://en.wikipedia.org/wiki/Spherical_triangle
  4. Üçlü skaler çarpım: https://en.wikipedia.org/wiki/Triple_product

Resumen del Problema

Para cada radio entero \(1 \le r \le 50\), se consideran los puntos de red sobre la esfera \(x^2+y^2+z^2=r^2\). Entre todas las ternas no degeneradas de esos puntos, \(A(r)\) es el área esférica mínima. Hay que calcular

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

Enfoque Matemático

Fijado un radio \(r\), todos los vértices permitidos pertenecen a

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

La solución tiene tres pasos: generar \(\mathcal P_r\), calcular el área esférica de cada terna y conservar la menor área positiva.

1) Enumeración de Puntos de Red sobre la Esfera

El código recorre \(x,y\in[-r,r]\), calcula

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

y solo acepta \((x,y)\) cuando \(z^2\) es no negativo y además un cuadrado perfecto. Si \(z\neq 0\), se añaden tanto \(z\) como \(-z\). Así se generan exactamente todos los puntos enteros de la esfera.

Para \(r\le 50\), el conjunto obtenido es lo bastante pequeño como para permitir una búsqueda directa sobre todas las ternas.

2) Área a partir del Ángulo Sólido

Sean \(a,b,c\in\mathcal P_r\). Tras normalizar,

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

los tres vectores quedan sobre la esfera unitaria. El área del triángulo esférico en la esfera de radio \(r\) es \(r^2\Omega\), donde \(\Omega\) es el ángulo sólido visto desde el origen.

Para vectores unitarios se usa la identidad

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

Al volver a escalar a radio \(r\), se obtiene

$$\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).$$

Por tanto, la fórmula exacta implementada por el código es

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

Cuando \(D>0\), esto coincide con \(2\arctan(\Delta/D)\), pero la forma con `atan2` es más segura numéricamente porque mantiene separados numerador y denominador.

3) Por qué aparece el determinante

La cantidad \(|\det(a,b,c)|\) es el volumen del paralelepípedo generado por los tres vectores radiales. Vale cero exactamente cuando los tres vectores son coplanarios con el origen, es decir, cuando el triángulo esférico es degenerado. Por eso el código descarta todas las ternas con

$$\Delta=0.$$

4) Por qué basta minimizar \(\Delta/D\)

En los candidatos que conserva el código se cumple \(\Delta>0\) y \(D>0\). En ese dominio, la función \(x\mapsto 2\arctan x\) es estrictamente creciente. Entonces minimizar el área equivale a minimizar

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

Eso es crucial para la implementación: en lugar de calcular áreas en coma flotante para cada terna, el programa compara dos candidatos \(\Delta_1/D_1\) y \(\Delta_2/D_2\) mediante multiplicación cruzada exacta,

$$\Delta_1D_2<\Delta_2D_1.$$

La versión en C++ usa `__int128`, la versión en Java usa `BigInteger` y Python utiliza enteros de precisión arbitraria.

También se puede descartar con seguridad el caso \(D\le 0\). En esas ternas se tiene \(\Omega\ge\pi\), así que el área es al menos \(\pi r^2\). Pero los puntos de los ejes \((r,0,0)\), \((0,r,0)\), \((0,0,r)\) siempre existen y producen un triángulo de área \(\pi r^2/2\). Por lo tanto, el mínimo nunca puede provenir de \(D\le 0\).

5) Ejemplo: \(r=1\)

Para \(r=1\), los puntos de red son exactamente los seis puntos de los ejes \(\pm(1,0,0)\), \(\pm(0,1,0)\), \(\pm(0,0,1)\). Tomemos

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

Entonces

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

y todos los productos escalares por pares son \(0\), así que \(D=1\). Por consiguiente,

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

Es el triángulo esférico de un octante, con tres ángulos rectos, y sirve como comprobación directa de la fórmula.

Cómo Funciona el Código

Para cada radio \(r\), la implementación construye primero la lista de puntos enteros sobre la esfera. Luego precalcula todos los productos escalares por pares en una tabla \(n_r\times n_r\), donde \(n_r=|\mathcal P_r|\). Así, el denominador

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

queda disponible en \(O(1)\) para cada terna.

Después recorre todas las ternas \(i<j<k\), omite las degeneradas, ignora \(D\le 0\) y conserva la mejor fracción \(\Delta/D\) usando comparación entera exacta. Solo al final evalúa una única llamada a `atan2` y multiplica por \(r^2\). El código C++ además verifica el punto de control publicado \(A(14)\approx 3.294040\) antes de sumar todos los radios.

Análisis de Complejidad

Sea \(n_r=|\mathcal P_r|\). La generación de puntos cuesta \(O(r^2)\) pares de prueba \((x,y)\), la precalculación de productos escalares cuesta \(O(n_r^2)\) tiempo y memoria, y el barrido dominante de ternas cuesta \(O(n_r^3)\) tiempo.

Como el problema solo llega hasta \(r=50\), todos los valores de \(n_r\) se mantienen lo bastante pequeños y esta búsqueda geométrica por fuerza bruta es completamente viable. Las comparaciones exactas en enteros también evitan inestabilidades entre lenguajes.

Lecturas Adicionales

  1. Página del problema: https://projecteuler.net/problem=332
  2. Ángulo sólido: https://en.wikipedia.org/wiki/Solid_angle
  3. Triángulo esférico y exceso esférico: https://en.wikipedia.org/wiki/Spherical_triangle
  4. Producto triple escalar: https://en.wikipedia.org/wiki/Triple_product

问题概述

对每个整数半径 \(1 \le r \le 50\),考虑球面 \(x^2+y^2+z^2=r^2\) 上的整数格点。在所有非退化三点组中,记最小球面三角形面积为 \(A(r)\)。目标是计算

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

数学方法

固定半径 \(r\) 后,所有允许的顶点都落在

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

之中。整个解法分为三步:生成 \(\mathcal P_r\),计算每个三元组的球面面积,并保留最小的正面积。

1)枚举球面上的整数点

代码枚举 \(x,y\in[-r,r]\),计算

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

只有当 \(z^2\) 非负且是完全平方数时,\((x,y)\) 才能对应球面上的整数点。若 \(z\neq 0\),则同时加入 \(z\) 和 \(-z\)。这样就能恰好生成球面上的全部整数格点。

由于这里只需要处理 \(r\le 50\),所以得到的点集规模足够小,可以直接枚举所有三元组。

2)由立体角得到球面面积

取 \(a,b,c\in\mathcal P_r\)。归一化后

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

都位于单位球面上。半径为 \(r\) 的球面三角形面积等于 \(r^2\Omega\),其中 \(\Omega\) 是从原点看到的立体角。

对于单位向量,有标准公式

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

把它缩放回半径为 \(r\) 的向量后,得到

$$\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).$$

因此代码实际使用的精确公式是

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

当 \(D>0\) 时,这也等价于 \(2\arctan(\Delta/D)\)。不过 `atan2` 形式在数值上更稳健,因为它不会过早把分子分母做浮点除法。

3)为什么会出现行列式

\(|\det(a,b,c)|\) 正是三个半径向量张成的平行六面体体积。它为零当且仅当这三个向量与原点共面,也就是对应的球面三角形退化。因此代码会直接跳过所有

$$\Delta=0$$

的三元组。

4)为什么只需最小化 \(\Delta/D\)

在代码保留的候选中,\(\Delta>0\) 且 \(D>0\)。在这个范围内,函数 \(x\mapsto 2\arctan x\) 严格单调递增,所以最小化面积等价于最小化

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

这对实现很关键:程序不必为每个三元组都计算浮点面积,而是用精确的交叉乘法比较两个候选 \(\Delta_1/D_1\) 和 \(\Delta_2/D_2\):

$$\Delta_1D_2<\Delta_2D_1.$$

C++ 版本用 `__int128`,Java 版本用 `BigInteger`,Python 则直接依赖内建高精度整数。

\(D\le 0\) 的情况也可以安全跳过。因为此时 \(\Omega\ge\pi\),对应面积至少为 \(\pi r^2\)。而轴点 \((r,0,0)\)、\((0,r,0)\)、\((0,0,r)\) 总是存在,并且它们构成的三角形面积是 \(\pi r^2/2\)。所以最小值绝不可能来自 \(D\le 0\) 的候选。

5)示例:\(r=1\)

当 \(r=1\) 时,整数格点恰好就是六个坐标轴点 \(\pm(1,0,0)\)、\(\pm(0,1,0)\)、\(\pm(0,0,1)\)。取

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

此时

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

且两两点积都为 \(0\),所以 \(D=1\)。于是

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

这正是单位球面的一个八分体三角形,三个角都是直角,是检验公式是否正确的理想例子。

代码如何工作

对每个半径 \(r\),实现先构造球面上的整数点列表,然后预计算所有两两点积并存入 \(n_r\times n_r\) 表,其中 \(n_r=|\mathcal P_r|\)。这样每个三元组的分母

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

都可以在 \(O(1)\) 时间得到。

随后程序遍历所有 \(i<j<k\),跳过退化三元组,忽略 \(D\le 0\) 的情况,并通过精确整数比较维护当前最优的 \(\Delta/D\)。只有在找到最优候选后,才做一次最终的 `atan2` 计算并乘上 \(r^2\)。C++ 代码在累加所有半径之前,还会先验证公开检查点 \(A(14)\approx 3.294040\)。

复杂度分析

记 \(n_r=|\mathcal P_r|\)。生成格点需要 \(O(r^2)\) 个 \((x,y)\) 试探,点积预计算需要 \(O(n_r^2)\) 的时间和空间,而主导步骤是三元组枚举,其时间复杂度为 \(O(n_r^3)\)。

由于本题只要求到 \(r=50\),所以每个 \(n_r\) 都不大,直接进行这种几何暴力搜索完全可行。整个搜索阶段都使用精确整数比较,也让不同语言版本更稳定一致。

延伸阅读

  1. 题目页面: https://projecteuler.net/problem=332
  2. 立体角: https://en.wikipedia.org/wiki/Solid_angle
  3. 球面三角形与球面超额: https://en.wikipedia.org/wiki/Spherical_triangle
  4. 标量三重积: https://en.wikipedia.org/wiki/Triple_product

Краткое описание

Для каждого целого радиуса \(1 \le r \le 50\) рассматриваются решетчатые точки на сфере \(x^2+y^2+z^2=r^2\). Среди всех невырожденных троек таких точек обозначим через \(A(r)\) минимальную площадь сферического треугольника. Требуется вычислить

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

Математический подход

При фиксированном \(r\) все допустимые вершины принадлежат множеству

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

Решение состоит из трех частей: построить \(\mathcal P_r\), вычислить площадь для каждой тройки и оставить минимальную положительную площадь.

1) Перебор решетчатых точек на сфере

Код перебирает \(x,y\in[-r,r]\), вычисляет

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

и принимает пару \((x,y)\) только тогда, когда \(z^2\) неотрицательно и является точным квадратом. Если \(z\neq 0\), добавляются обе точки с координатами \(z\) и \(-z\). Так получаются ровно все целочисленные точки сферы.

Для \(r\le 50\) множество точек остается достаточно маленьким, поэтому полный перебор троек вполне реалистичен.

2) Площадь через телесный угол

Пусть \(a,b,c\in\mathcal P_r\). После нормировки

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

они лежат на единичной сфере. Площадь сферического треугольника на сфере радиуса \(r\) равна \(r^2\Omega\), где \(\Omega\) — телесный угол, видимый из начала координат.

Для единичных векторов используется стандартная формула

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

После обратного масштабирования к векторам длины \(r\) получаем

$$\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).$$

Следовательно, точная формула, используемая в коде, такова:

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

При \(D>0\) это совпадает с \(2\arctan(\Delta/D)\), но форма через `atan2` численно надежнее, потому что числитель и знаменатель не делятся заранее.

3) Почему появляется определитель

Величина \(|\det(a,b,c)|\) равна объему параллелепипеда, натянутого на три радиус-вектора. Она обращается в ноль тогда и только тогда, когда эти векторы компланарны с началом координат, то есть когда сферический треугольник вырожден. Поэтому код сразу пропускает все тройки с

$$\Delta=0.$$

4) Почему достаточно минимизировать \(\Delta/D\)

Для кандидатов, которые реально рассматривает код, выполнено \(\Delta>0\) и \(D>0\). На этой области функция \(x\mapsto 2\arctan x\) строго возрастает. Значит, минимизация площади эквивалентна минимизации отношения

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

Это важно вычислительно: вместо того чтобы считать вещественную площадь для каждой тройки, программа сравнивает две дроби \(\Delta_1/D_1\) и \(\Delta_2/D_2\) точным перекрестным умножением,

$$\Delta_1D_2<\Delta_2D_1.$$

Версия на C++ использует `__int128`, версия на Java — `BigInteger`, а Python опирается на встроенные длинные целые.

Случай \(D\le 0\) тоже можно безопасно отбросить. Тогда \(\Omega\ge\pi\), а значит площадь не меньше \(\pi r^2\). Но точки \((r,0,0)\), \((0,r,0)\), \((0,0,r)\) всегда присутствуют и дают треугольник площади \(\pi r^2/2\). Следовательно, минимум никогда не достигается при \(D\le 0\).

5) Пример: \(r=1\)

При \(r=1\) решетчатые точки — это ровно шесть осевых точек \(\pm(1,0,0)\), \(\pm(0,1,0)\), \(\pm(0,0,1)\). Возьмем

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

Тогда

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

а все попарные скалярные произведения равны \(0\), так что \(D=1\). Получаем

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

Это сферический треугольник одного октанта с тремя прямыми углами, и он служит удобной проверкой корректности формулы.

Как работает код

Для каждого \(r\) программа сначала строит список решетчатых точек на сфере. Затем она заранее вычисляет все попарные скалярные произведения и сохраняет их в таблице размера \(n_r\times n_r\), где \(n_r=|\mathcal P_r|\). Благодаря этому знаменатель

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

можно получить за \(O(1)\) для каждой тройки.

После этого код перебирает все \(i<j<k\), пропускает вырожденные тройки, игнорирует случаи \(D\le 0\) и поддерживает лучшую дробь \(\Delta/D\) с помощью точного целочисленного сравнения. Только после выбора победителя вызывается один `atan2`, а результат умножается на \(r^2\). C++-версия также проверяет опубликованный контрольный ориентир \(A(14)\approx 3.294040\) перед суммированием по всем радиусам.

Анализ сложности

Пусть \(n_r=|\mathcal P_r|\). Генерация точек требует \(O(r^2)\) пробных пар \((x,y)\), предвычисление скалярных произведений требует \(O(n_r^2)\) времени и памяти, а доминирующий перебор троек занимает \(O(n_r^3)\) времени.

Так как в задаче нужно рассматривать только радиусы до \(50\), каждое \(n_r\) остается небольшим, и такой полный геометрический перебор легко укладывается в разумное время. Точные целочисленные сравнения также обеспечивают стабильность результата во всех языках реализации.

Дополнительные материалы

  1. Страница задачи: https://projecteuler.net/problem=332
  2. Телесный угол: https://en.wikipedia.org/wiki/Solid_angle
  3. Сферический треугольник и сферический избыток: https://en.wikipedia.org/wiki/Spherical_triangle
  4. Смешанное произведение: https://en.wikipedia.org/wiki/Triple_product

ملخص المسألة

لكل نصف قطر صحيح \(1 \le r \le 50\)، ننظر إلى نقاط الشبكة الصحيحة على الكرة \(x^2+y^2+z^2=r^2\). بين جميع الثلاثيات غير المنحلة من هذه النقاط، نرمز إلى أصغر مساحة مثلث كروي بالرمز \(A(r)\). المطلوب هو حساب

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

المنهج الرياضي

عند تثبيت \(r\)، فإن جميع الرؤوس المسموح بها تنتمي إلى المجموعة

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

يتكون الحل من ثلاث مراحل: توليد \(\mathcal P_r\)، حساب المساحة الكروية لكل ثلاثية، ثم الاحتفاظ بأصغر مساحة موجبة.

1) توليد نقاط الشبكة على الكرة

تقوم الشيفرة بالمرور على \(x,y\in[-r,r]\)، ثم تحسب

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

ولا تقبل الزوج \((x,y)\) إلا إذا كان \(z^2\) غير سالب وكان مربعًا كاملاً. وإذا كان \(z\neq 0\)، تُضاف النقطتان المقابلتان \(z\) و\(-z\). بهذه الطريقة نحصل بدقة على جميع النقاط الصحيحة الواقعة على الكرة.

وبما أن \(r\le 50\)، فإن عدد النقاط الناتج يظل صغيرًا بما يكفي للمرور المباشر على جميع الثلاثيات.

2) اشتقاق المساحة من الزاوية المجسمة

لتكن \(a,b,c\in\mathcal P_r\). بعد التطبيع تصبح

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

نقاطًا على الكرة الواحدة. مساحة المثلث الكروي على كرة نصف قطرها \(r\) تساوي \(r^2\Omega\)، حيث \(\Omega\) هي الزاوية المجسمة المرئية من الأصل.

ولمتجهات الوحدة توجد الصيغة القياسية

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

وبالعودة إلى المتجهات ذات الطول \(r\) نحصل على

$$\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).$$

إذن الصيغة الدقيقة المستخدمة في الشيفرة هي

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

وعندما يكون \(D>0\)، فهذا يساوي أيضًا \(2\arctan(\Delta/D)\)، لكن صيغة `atan2` أكثر أمانًا عدديًا لأنها لا تحسب القسمة العائمة مبكرًا.

3) لماذا يظهر المحدد

الكمية \(|\det(a,b,c)|\) هي حجم متوازي السطوح المولَّد من المتجهات الشعاعية الثلاثة. وهي تساوي صفرًا فقط عندما تكون هذه المتجهات واقعة في مستوى واحد مع الأصل، أي عندما يكون المثلث الكروي منحلًا. لذلك تتجاوز الشيفرة كل ثلاثية تحقق

$$\Delta=0.$$

4) لماذا يكفي تصغير \(\Delta/D\)

في جميع المرشحين الذين تحتفظ بهم الشيفرة لدينا \(\Delta>0\) و\(D>0\). وعلى هذا المجال تكون الدالة \(x\mapsto 2\arctan x\) متزايدة تمامًا، لذا فإن تصغير المساحة يكافئ تصغير النسبة

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

وهذا مهم من الناحية الحسابية: بدلًا من حساب مساحة عائمة لكل ثلاثية، تقارن الشيفرة بين مرشحين \(\Delta_1/D_1\) و\(\Delta_2/D_2\) بالضرب التبادلي الدقيق:

$$\Delta_1D_2<\Delta_2D_1.$$

نسخة C++ تستخدم `__int128`، ونسخة Java تستخدم `BigInteger`، أما Python فتعتمد على الأعداد الصحيحة كبيرة الدقة بشكل طبيعي.

كما يمكن استبعاد الحالة \(D\le 0\) بأمان. ففي هذه الحالة يكون \(\Omega\ge\pi\)، وبالتالي تكون المساحة على الأقل \(\pi r^2\). لكن النقاط المحورية \((r,0,0)\) و\((0,r,0)\) و\((0,0,r)\) موجودة دائمًا، وهي تعطي مثلثًا مساحته \(\pi r^2/2\). لذا لا يمكن أن يأتي الحد الأدنى من مرشح يحقق \(D\le 0\).

5) مثال: \(r=1\)

عندما \(r=1\)، تكون نقاط الشبكة هي بالضبط النقاط الست على المحاور \(\pm(1,0,0)\) و\(\pm(0,1,0)\) و\(\pm(0,0,1)\). خذ

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

عندئذٍ

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

وجميع الضربات النقطية الثنائية تساوي \(0\)، وبالتالي \(D=1\). ومن ثم

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

وهذا هو مثلث الأوكتانت الكروي المعروف ذو الزوايا الثلاث القائمة، وهو اختبار ممتاز لصحة الصيغة.

كيف تعمل الشيفرة

لكل \(r\)، تبني الشيفرة أولًا قائمة نقاط الشبكة على الكرة. ثم تحسب مسبقًا جميع الضربات النقطية الثنائية وتخزنها في جدول بحجم \(n_r\times n_r\)، حيث \(n_r=|\mathcal P_r|\). وبذلك يصبح المقام

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

متاحًا في زمن \(O(1)\) لكل ثلاثية.

بعد ذلك تمر الشيفرة على جميع \(i<j<k\)، وتتجاوز الثلاثيات المنحلة، وتهمل الحالات \(D\le 0\)، وتحافظ على أفضل كسر \(\Delta/D\) باستخدام مقارنة صحيحة دقيقة. ولا تُجرى عملية `atan2` إلا مرة واحدة بعد معرفة المرشح الفائز، ثم يُضرب الناتج في \(r^2\). كما تتحقق نسخة C++ من نقطة الفحص المنشورة \(A(14)\approx 3.294040\) قبل جمع النتائج لجميع أنصاف الأقطار.

تحليل التعقيد

لنكتب \(n_r=|\mathcal P_r|\). توليد النقاط يحتاج إلى \(O(r^2)\) من أزواج الاختبار \((x,y)\)، وحساب الضربات النقطية مسبقًا يحتاج إلى \(O(n_r^2)\) زمنًا وذاكرة، أما المسح الثلاثي المسيطر فيحتاج إلى \(O(n_r^3)\) زمنًا.

وبما أن المسألة تتوقف عند \(r=50\)، فإن كل \(n_r\) يبقى صغيرًا بما يكفي لجعل هذا البحث الهندسي المباشر عمليًا تمامًا. كما أن المقارنات الصحيحة الدقيقة تجعل النتيجة مستقرة بين اللغات المختلفة.

مراجع إضافية

  1. صفحة المسألة: https://projecteuler.net/problem=332
  2. الزاوية المجسمة: https://en.wikipedia.org/wiki/Solid_angle
  3. المثلث الكروي والفائض الكروي: https://en.wikipedia.org/wiki/Spherical_triangle
  4. الجداء الثلاثي: https://en.wikipedia.org/wiki/Triple_product