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

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

Problemzusammenfassung

Für jedes \(n\) wird \(r=2^n-1\) gesetzt und die Menge der ganzzahligen Kugelpunkte

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

betrachtet. Nord- und Südpol sind \(N=(0,0,r)\) und \(S=(0,0,-r)\). Ein Sprung zwischen zwei Gitterpunkten hat ein Risiko, das als Quadrat des normierten Zentriwinkels definiert ist. Bezeichnet \(M(r)\) das minimale Gesamtrisiko von \(N\) nach \(S\), dann lautet die gesuchte Summe

$$\sum_{n=1}^{15} M(2^n-1).$$

Die lokalen Lösungen modellieren das Problem als exakten dicht gewichteten Graphen und beschleunigen ihn anschließend durch Symmetrieausnutzung.

Mathematischer Ansatz

Schritt 1: Kantengewicht aus der Kugelgeometrie

Für \(a,b\in S_r\) gilt \(\lVert a\rVert=\lVert b\rVert=r\). Der Zentriwinkel \(\theta(a,b)\) erfüllt daher

$$\cos\theta(a,b)=\frac{a\cdot b}{r^2},\qquad \theta(a,b)=\arccos\left(\frac{a\cdot b}{r^2}\right).$$

Der Code normiert mit \(\pi\), klemmt den Kosinus aus Stabilitätsgründen auf \([-1,1]\) und setzt

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

Genau diese Formel steht in edge_risk beziehungsweise edgeRisk.

Schritt 2: Das Problem ist ein Kürzester-Weg-Problem

Da die Implementierungen jeden Übergang zwischen zwei Kugelpunkten erlauben, entsteht ein vollständiger gewichteter Graph auf \(S_r\). Für einen Weg

$$P=(p_0,p_1,\dots,p_k),\qquad p_0=N,\quad p_k=S,$$

ist das Gesamtrisiko

$$R(P)=\sum_{i=0}^{k-1} w(p_i,p_{i+1}).$$

Also

$$M(r)=\min_{P:N\leadsto S} R(P).$$

Die exakte C++-Kontrollrechnung verwendet dafür Dijkstra auf allen durch generate_all_points erzeugten Punkten.

Schritt 3: Kleines Beispiel für die Risikostruktur

Für \(r=1\) gibt es nur die sechs Achsenpunkte. Der direkte Sprung vom Nord- zum Südpol hat \(\theta=\pi\) und damit Risiko \(1\). Über einen Äquatorpunkt erhält man zwei Vierteldrehungen mit normiertem Winkel \(1/2\), also

$$M(1)=2\left(\frac{1}{2}\right)^2=\frac{1}{2}.$$

Das zeigt den zentralen Effekt: Viele kleine Winkel sind oft billiger als ein großer, weil das Risiko quadratisch wächst.

Schritt 4: Spiegelbrücke der schnellen Methode

Die optimierte Suche berechnet zunächst minimale Risiken vom Nordpol zu einer symmetriereduzierten Punktmenge in der oberen Halbkugel. Sei

$$p=(x,y,z),\qquad z\ge 0,$$

ein solcher Repräsentant und

$$p'=(x,y,-z)$$

sein Spiegel an der Äquatorebene. Ist \(d(p)\) das minimale Risiko von \(N\) nach \(p\), dann liefert die Symmetrie denselben Wert von \(p'\) nach \(S\). Ein kompletter Kandidat ist also

$$2d(p)+w(p,p').$$

Wegen

$$p\cdot p'=x^2+y^2-z^2=r^2-2z^2$$

folgt

$$\cos\theta=\frac{r^2-2z^2}{r^2}=1-2\left(\frac{z}{r}\right)^2.$$

Mit \(\alpha=\arcsin(z/r)\) gilt \(\cos\theta=\cos(2\alpha)\), also

$$\theta=2\arcsin\left(\frac{z}{r}\right).$$

Damit wird die Spiegelkante zu

$$w(p,p')=\left(\frac{2\arcsin(z/r)}{\pi}\right)^2,$$

und die schnelle Methode minimiert

$$\boxed{M(r)=\min_{p}\left(2d(p)+\left(\frac{2\arcsin(z_p/r)}{\pi}\right)^2\right).}$$

Schritt 5: Reduzierte Punktmenge und konservatives Pruning

generate_reduced_points speichert nicht alle Vorzeichen- und Permutationskopien. Stattdessen werden geordnete nichtnegative Tripel mit

$$0\le x\le y\le z,\qquad x^2+y^2+z^2=r^2$$

durchlaufen und nur die für das Symmetriemodell nötigen Permutationen ausgegeben. Die Bedingung \(3x^2<r^2\) stammt direkt aus \(x\le y\le z\). Nach Sortierung und Deduplikation bleibt die reduzierte Zustandsmenge für die schnelle Suche.

Zusätzlich wird

$$a_t=\frac{\arcsin(t/r)}{\pi}\qquad (0\le t\le r)$$

vorab berechnet. Hat man aktuell eine obere Schranke \(B\) für den gesamten Nord-Süd-Weg und einen einseitigen Risikowert \(d_c\), dann verwendet der Code das konservative Restbudget \(B-2d_c\). Billige Untergrenzen aus Differenzen der \(a_t\)-Werte werden vor dem teuren acos-Aufruf getestet. C++ und Java prüfen dabei \(z\), \(x\) und \(y\); Python behält dieselbe Struktur, nutzt aber nur die \(z\)-Schranke.

Schritt 6: Grob-zu-fein und Validierung

Die schnelle Lösung arbeitet in vier Durchgängen mit den Strides \(27,9,3,1\). Ein grober Durchgang liefert eine obere Schranke, die den nächsten feineren Durchgang billiger macht. Daher ruft minimal_risk_fast dieselbe geprunte Dijkstra-Routine mehrfach auf immer dichteren Stichproben auf.

Die C++-Datei prüft zwei Kontrollwerte. Erstens wird

$$M(7)=0.1784943998$$

verlangt. Zweitens wird für \(r=31\) der exakte vollständige Graph mit der schnellen symmetriereduzierten Methode verglichen; die Abweichung muss unter \(10^{-12}\) bleiben. So wird abgesichert, dass die Optimierung dieselbe Mathematik löst wie das exakte Modell.

Funktionsweise des Codes

Die C++-Quelle enthält sowohl den exakten Prüfer als auch den optimierten Solver. Python und Java implementieren direkt die reduzierte schnelle Variante, aber mit denselben mathematischen Bausteinen: Kantengewicht, vorab berechnete \(\arcsin(z/r)/\pi\)-Tabelle, Erzeugung der Repräsentantenpunkte und die Spiegelabschlussformel.

In der Ausführung unterscheiden sich die Sprachen leicht. C++ verteilt die 15 Radien auf bis zu acht pthread-Worker, Java verwendet einen parallelen IntStream, und Python rechnet seriell. Die exakten Hilfsfunktionen generate_all_points, pole_indices und dijkstra_complete existieren nur in C++ und dienen dort nur der Validierung.

Komplexitätsanalyse

Sei \(V_r=|S_r|\). Die exakte Dijkstra-Variante auf dem vollständigen Graphen braucht \(O(V_r^2)\) Zeit und \(O(V_r)\) Speicher, weil in jedem Schritt alle verbleibenden Knoten relaxiert werden. Für Kontrollfälle wie \(r=31\) ist das akzeptabel, für die größten Radien der Endsumme aber nicht mehr.

Die schnelle Methode bleibt im Worst Case quadratisch in der Größe der reduzierten Punktmenge. Praktisch ist sie jedoch deutlich schneller, weil die reduzierte Menge viel kleiner als \(S_r\) ist, die Strides \(27,9,3,1\) immer bessere Schranken liefern und das konservative Pruning viele Kandidaten gar nicht mehr vollständig bewertet.

Fußnoten und Quellen

  1. Problemseite: https://projecteuler.net/problem=353
  2. Großkreisabstand und Zentriwinkel: Wikipedia — Großkreis
  3. Dijkstra-Algorithmus: Wikipedia — Dijkstra-Algorithmus
  4. Gitterpunkte auf Kugeln und Drei-Quadrate-Satz: Wikipedia — Drei-Quadrate-Satz
  5. Kugelkosinussatz: Wikipedia — Kugelkosinussatz

Problem Özeti

Her \(n\) için \(r=2^n-1\) alınır ve küre üzerindeki tamsayı noktalar

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

tanımlanır. Kuzey ve güney kutupları \(N=(0,0,r)\) ile \(S=(0,0,-r)\) noktalarıdır. İki kafes noktası arasındaki hamlenin riski, merkez açısının \(\pi\) ile normalize edilip karesinin alınmasıyla ölçülür. Eğer \(M(r)\), \(N\)'den \(S\)'ye giden en düşük toplam riski gösteriyorsa, istenen toplam

$$\sum_{n=1}^{15} M(2^n-1)$$

olur. Yerel çözüm dosyaları bu hesabı önce tam bir yoğun grafik modeliyle kurup sonra simetri indirgemesiyle hızlandırır.

Matematiksel Yaklaşım

Adım 1: Küresel Geometriden Kenar Riski

\(a,b\in S_r\) için \(\lVert a\rVert=\lVert b\rVert=r\) olduğundan merkez açı \(\theta(a,b)\) şu bağıntıyı sağlar:

$$\cos\theta(a,b)=\frac{a\cdot b}{r^2},\qquad \theta(a,b)=\arccos\left(\frac{a\cdot b}{r^2}\right).$$

Kod, sayısal kararlılık için kosinüsü \([-1,1]\) aralığına sıkıştırır ve

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

formülünü kullanır. C++, Python ve Java dosyalarındaki edge_risk/edgeRisk tanımı tam olarak budur.

Adım 2: Problem Yoğun Bir En Kısa Yol Problemidir

Çözüm dosyaları herhangi iki küre noktasını doğrudan birbirine bağlayabildiği için, \(S_r\) üzerinde tam ağırlıklı bir grafik oluşur. Eğer

$$P=(p_0,p_1,\dots,p_k),\qquad p_0=N,\quad p_k=S$$

bir yol ise toplam risk

$$R(P)=\sum_{i=0}^{k-1} w(p_i,p_{i+1})$$

olur. Dolayısıyla

$$M(r)=\min_{P:N\leadsto S} R(P)$$

ve çekirdek problem en kısa yol hesabına dönüşür. C++ içindeki tam doğrulama, generate_all_points ile tüm düğümleri üretip tam graf üzerinde Dijkstra çalıştırır.

Adım 3: Basit Bir Örnek

\(r=1\) olduğunda sadece eksen üzerindeki altı nokta vardır. Kuzeyden güneye doğrudan atlama için \(\theta=\pi\) olur ve risk \(1\)'dir. Bir ekvator noktasından geçilirse iki tane çeyrek dönüş yapılır; her birinin normalize açısı \(1/2\) olduğundan

$$M(1)=2\left(\frac{1}{2}\right)^2=\frac{1}{2}.$$

Bu örnek, neden çok sayıda küçük açının tek büyük açıdan daha ucuz olabildiğini açıkça gösterir: maliyet açıya lineer değil, karesel bağlıdır.

Adım 4: Hızlı Çözümdeki Ayna Köprüsü

Optimizasyonlu algoritma önce üst yarıküredeki simetriyle azaltılmış temsilci noktalara kadar olan en kısa riskleri hesaplar. Böyle bir nokta

$$p=(x,y,z),\qquad z\ge 0$$

ve ekvator düzlemine göre aynası

$$p'=(x,y,-z)$$

olsun. Eğer \(d(p)\), kuzeyden \(p\)'ye minimum risk ise simetriden dolayı \(p'\)'den güneye olan minimum risk de aynıdır. O halde tam bir aday yolun maliyeti

$$2d(p)+w(p,p')$$

şeklindedir. Ayrıca

$$p\cdot p'=x^2+y^2-z^2=r^2-2z^2$$

olduğu için

$$\cos\theta=\frac{r^2-2z^2}{r^2}=1-2\left(\frac{z}{r}\right)^2.$$

\(\alpha=\arcsin(z/r)\) dersek \(\cos\theta=\cos(2\alpha)\) elde edilir ve

$$\theta=2\arcsin\left(\frac{z}{r}\right)$$

olur. Böylece ayna köprüsünün riski

$$w(p,p')=\left(\frac{2\arcsin(z/r)}{\pi}\right)^2$$

ve hızlı çözücünün minimize ettiği ifade

$$\boxed{M(r)=\min_{p}\left(2d(p)+\left(\frac{2\arcsin(z_p/r)}{\pi}\right)^2\right).}$$

Adım 5: Azaltılmış Nokta Kümesi ve Budama

generate_reduced_points bütün işaret ve permütasyon kopyalarını tutmaz. Bunun yerine

$$0\le x\le y\le z,\qquad x^2+y^2+z^2=r^2$$

koşulunu sağlayan sıralı, negatif olmayan üçlüleri tarar; sonra kullanılan simetri modeline uygun gerekli permütasyonları üretir. \(3x^2<r^2\) dış sınırı, doğrudan \(x\le y\le z\) koşulundan gelir. Sıralama ve tekrar silmeden sonra hızlı Dijkstra'nın düğüm kümesi elde edilir.

Buna ek olarak

$$a_t=\frac{\arcsin(t/r)}{\pi}\qquad (0\le t\le r)$$

önceden hesaplanır. Eğer eldeki en iyi kuzey-güney üst sınırı \(B\), mevcut tek taraflı risk de \(d_c\) ise kod, muhafazakar kalan bütçe olarak \(B-2d_c\) kullanır. \(a_t\) farklarından türetilen ucuz alt sınırlar pahalı acos hesabından önce denenir. C++ ve Java sürümleri bunu \(z\), \(x\) ve \(y\) için uygular; Python aynı ana yapıyı korur fakat yalnızca \(z\)-tabanlı testi kullanır.

Adım 6: Kademeli Yoğunlaştırma ve Kontrol Noktaları

Hızlı çözücü, azaltılmış nokta kümesi üzerinde \(27,9,3,1\) adım aralıklarıyla dört geçiş yapar. Kaba geçişler üst sınır üretir; bu sınırlar daha yoğun sonraki geçişleri ucuzlatır. Bu yüzden minimal_risk_fast, aynı budanmış Dijkstra yordamını giderek daha sık örneklenmiş listeler üzerinde tekrar çağırır.

C++ kaynağı iki önemli doğrulama yapar. Birincisi,

$$M(7)=0.1784943998$$

değerini kontrol eder. İkincisi, \(r=31\) için tam yoğun grafik Dijkstra sonucu ile hızlı simetri-indirgenmiş sonucun \(10^{-12}\) hassasiyet içinde aynı olmasını ister. Böylece hızlı yöntem, farklı bir yaklaşım değil, tam modelin güvenli bir optimizasyonu olarak doğrulanmış olur.

Kodun Çalışma Mantığı

C++ dosyası hem tam doğrulayıcıyı hem de optimize edilmiş çözücüyü içerir. Python ve Java ise doğrudan azaltılmış hızlı yöntemi uygular. Buna rağmen üç dilde de aynı matematiksel parçalar bulunur: kenar riski, \(\arcsin(z/r)/\pi\) tablosu, temsilci nokta üretimi ve ayna ile kapanış formülü.

Uygulama ayrıntılarında küçük farklar vardır. C++, 15 yarıçapı en fazla sekiz pthread iş parçacığına dağıtır; Java paralel IntStream kullanır; Python seri döngü çalıştırır. Tam doğrulama için kullanılan generate_all_points, pole_indices ve dijkstra_complete yordamları yalnızca C++ kaynağında vardır.

Karmaşıklık Analizi

\(V_r=|S_r|\) olsun. Tam yoğun grafik Dijkstra yaklaşımı \(O(V_r^2)\) zaman ve \(O(V_r)\) bellek kullanır; çünkü her adımda geride kalan bütün düğümler gevşetilir. Bu, \(r=31\) gibi doğrulama örnekleri için uygundur ama son toplamın en büyük yarıçapları için pahalıdır.

Hızlı yöntem de en kötü durumda azaltılmış düğüm sayısına göre karesel kalır. Fakat azaltılmış küme tam kafes nokta kümesinden çok daha küçüktür, \(27,9,3,1\) kademeleri gittikçe daha iyi üst sınırlar sağlar ve muhafazakar budama çok sayıda kenarı hiç tam değerlendirmeden eler. Pratik hız kazancı bu üç etkiden gelir.

Dipnotlar ve Kaynakça

  1. Problem sayfası: https://projecteuler.net/problem=353
  2. Büyük çember uzaklığı ve merkez açı: Wikipedia — Great-circle distance
  3. Dijkstra algoritması: Wikipedia — Dijkstra algoritması
  4. Küre üzerindeki tamsayı noktalar ve üç kare teoremi: Wikipedia — Sum of three squares theorem
  5. Küresel kosinüs bağıntısı: Wikipedia — Spherical law of cosines

Resumen del Problema

Para cada \(n\), el problema fija \(r=2^n-1\) y considera los puntos enteros de la esfera

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

Los polos norte y sur son \(N=(0,0,r)\) y \(S=(0,0,-r)\). Un movimiento entre dos puntos de red tiene un riesgo igual al cuadrado del ángulo central normalizado por \(\pi\). Si \(M(r)\) es el riesgo total mínimo para ir de \(N\) a \(S\), entonces la cantidad pedida es

$$\sum_{n=1}^{15} M(2^n-1).$$

Los archivos locales resuelven esto combinando un modelo exacto de grafo denso con una versión rápida basada en simetrías.

Enfoque Matemático

Paso 1: Riesgo de una arista a partir del ángulo central

Si \(a,b\in S_r\), entonces \(\lVert a\rVert=\lVert b\rVert=r\). El ángulo central \(\theta(a,b)\) cumple

$$\cos\theta(a,b)=\frac{a\cdot b}{r^2},\qquad \theta(a,b)=\arccos\left(\frac{a\cdot b}{r^2}\right).$$

Las implementaciones normalizan por \(\pi\), recortan el coseno al intervalo \([-1,1]\) por estabilidad numérica y definen

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

Esa es exactamente la fórmula usada en edge_risk y edgeRisk.

Paso 2: Convertir la luna en un grafo completo ponderado

Como el código permite pasar directamente entre cualesquiera dos puntos de \(S_r\), obtenemos un grafo completo. Para un camino

$$P=(p_0,p_1,\dots,p_k),\qquad p_0=N,\quad p_k=S,$$

el riesgo total es

$$R(P)=\sum_{i=0}^{k-1} w(p_i,p_{i+1}).$$

Por tanto,

$$M(r)=\min_{P:N\leadsto S} R(P).$$

La rutina exacta de C++ resuelve este problema con Dijkstra sobre todos los puntos generados por generate_all_points.

Paso 3: Un ejemplo pequeño

Cuando \(r=1\), solo existen los seis puntos de los ejes coordenados. El salto directo del norte al sur tiene \(\theta=\pi\), así que su riesgo es \(1\). Si se pasa por un punto ecuatorial, aparecen dos giros de un cuarto de vuelta; cada uno tiene ángulo normalizado \(1/2\), luego

$$M(1)=2\left(\frac{1}{2}\right)^2=\frac{1}{2}.$$

Este caso muestra por qué dividir un gran salto en varios saltos cortos puede ser mejor: la función cuadrática castiga mucho los ángulos grandes.

Paso 4: Fórmula del espejo en el solucionador rápido

La versión optimizada calcula primero riesgos mínimos desde el polo norte hasta un conjunto reducido de puntos representativos en el hemisferio superior. Sea

$$p=(x,y,z),\qquad z\ge 0,$$

y su reflejo respecto del ecuador

$$p'=(x,y,-z).$$

Si \(d(p)\) es el riesgo mínimo desde \(N\) hasta \(p\), entonces por simetría el riesgo mínimo desde \(p'\) hasta \(S\) es el mismo. Por ello, un candidato completo norte-sur vale

$$2d(p)+w(p,p').$$

Además,

$$p\cdot p'=x^2+y^2-z^2=r^2-2z^2,$$

de modo que

$$\cos\theta=\frac{r^2-2z^2}{r^2}=1-2\left(\frac{z}{r}\right)^2.$$

Si \(\alpha=\arcsin(z/r)\), entonces \(\cos\theta=\cos(2\alpha)\), y por tanto

$$\theta=2\arcsin\left(\frac{z}{r}\right).$$

Así, el riesgo del puente espejo es

$$w(p,p')=\left(\frac{2\arcsin(z/r)}{\pi}\right)^2,$$

y el algoritmo rápido minimiza

$$\boxed{M(r)=\min_{p}\left(2d(p)+\left(\frac{2\arcsin(z_p/r)}{\pi}\right)^2\right).}$$

Paso 5: Conjunto reducido y poda conservadora

generate_reduced_points no conserva todas las copias por signos y permutaciones. En su lugar recorre ternas ordenadas no negativas con

$$0\le x\le y\le z,\qquad x^2+y^2+z^2=r^2,$$

y luego emite solo las permutaciones necesarias para el modelo de simetría usado por el programa. La cota \(3x^2<r^2\) viene directamente de \(x\le y\le z\). Después de ordenar y eliminar duplicados, esa lista reducida es el espacio de estados del buscador rápido.

También se precalcula

$$a_t=\frac{\arcsin(t/r)}{\pi}\qquad (0\le t\le r).$$

Si la mejor cota superior actual para el camino completo es \(B\) y el riesgo unilateral actual es \(d_c\), el código usa el presupuesto conservador restante \(B-2d_c\). Antes de llamar a acos, prueba cotas inferiores baratas obtenidas de diferencias de los valores \(a_t\). En C++ y Java esa prueba se aplica a \(z\), \(x\) e \(y\); en Python se mantiene la misma idea general, pero solo con la cota basada en \(z\).

Paso 6: Búsqueda coarse-to-fine y comprobaciones

El solucionador rápido hace cuatro pasadas con saltos \(27,9,3,1\). Las pasadas gruesas generan cotas superiores que abaratan las pasadas más finas. Por eso minimal_risk_fast llama varias veces al Dijkstra podado sobre muestras cada vez más densas del mismo conjunto reducido.

El archivo C++ valida el método de dos maneras. Primero comprueba el benchmark

$$M(7)=0.1784943998.$$

Después, para \(r=31\), ejecuta tanto el Dijkstra exacto del grafo completo como el solver rápido reducido, y exige coincidencia a \(10^{-12}\). Eso confirma que la optimización acelera el mismo modelo matemático exacto.

Cómo Funciona el Código

La fuente en C++ contiene tanto el validador exacto como el solver optimizado. Los archivos Python y Java implementan directamente la versión reducida rápida. En los tres casos aparecen las mismas piezas matemáticas: la fórmula del riesgo, la tabla precalculada \(\arcsin(z/r)/\pi\), la construcción de puntos representativos y la fórmula de cierre por reflexión.

Hay pequeñas diferencias de implementación. C++ reparte los 15 radios entre hasta ocho hilos pthread; Java usa un IntStream paralelo; Python recorre los radios en serie. Las funciones exactas generate_all_points, pole_indices y dijkstra_complete solo existen en C++ y se usan como verificación, no en la ejecución principal.

Complejidad

Si \(V_r=|S_r|\), el Dijkstra exacto sobre el grafo completo cuesta \(O(V_r^2)\) tiempo y \(O(V_r)\) memoria, porque cada extracción relaja todos los vértices restantes. Eso es razonable para comprobaciones como \(r=31\), pero demasiado caro para los radios mayores de la suma final.

El solucionador rápido sigue siendo cuadrático en el peor caso respecto del tamaño reducido, ya que todavía relaja un grafo denso. Sin embargo, el conjunto reducido es mucho menor que \(S_r\), el esquema \(27,9,3,1\) aporta cotas cada vez mejores y la poda conservadora evita evaluar muchas aristas completas. Esa combinación es la clave del rendimiento práctico.

Notas y Referencias

  1. Página del problema: https://projecteuler.net/problem=353
  2. Distancia de gran círculo y ángulo central: Wikipedia — Distancia ortodrómica
  3. Algoritmo de Dijkstra: Wikipedia — Algoritmo de Dijkstra
  4. Puntos enteros en esferas y suma de tres cuadrados: Wikipedia — Sum of three squares theorem
  5. Ley esférica de los cosenos: Wikipedia — Trigonometría esférica

问题概述

对每个 \(n\),题目令 \(r=2^n-1\),并考虑球面上的整点集合

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

北极和南极分别是 \(N=(0,0,r)\) 与 \(S=(0,0,-r)\)。从一个整点跳到另一个整点的风险,定义为对应球心夹角除以 \(\pi\) 后再平方。如果 \(M(r)\) 表示从 \(N\) 到 \(S\) 的最小总风险,那么题目要求的是

$$\sum_{n=1}^{15} M(2^n-1).$$

本仓库里的解法先把问题写成精确的稠密加权图最短路,再利用对称性把状态空间大幅压缩。

数学方法

步骤 1:由球面几何得到边权

若 \(a,b\in S_r\),则 \(\lVert a\rVert=\lVert b\rVert=r\)。设它们的球心夹角为 \(\theta(a,b)\),由点积公式可得

$$\cos\theta(a,b)=\frac{a\cdot b}{r^2},\qquad \theta(a,b)=\arccos\left(\frac{a\cdot b}{r^2}\right).$$

程序把角度除以 \(\pi\) 进行归一化,并把余弦值裁剪到 \([-1,1]\) 以避免浮点误差,然后定义

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

这正是三份源码里 edge_riskedgeRisk 的核心公式。

步骤 2:把问题变成完全图最短路

由于代码允许任意两个球面整点之间直接连边,所以 \(S_r\) 上得到的是一个完全加权图。对一条路径

$$P=(p_0,p_1,\dots,p_k),\qquad p_0=N,\quad p_k=S,$$

其总风险为

$$R(P)=\sum_{i=0}^{k-1} w(p_i,p_{i+1}).$$

于是

$$M(r)=\min_{P:N\leadsto S} R(P).$$

因此数学本质就是一个最短路问题。C++ 中的精确校验路径使用 generate_all_points 生成全部整点,再在完整图上运行 Dijkstra。

步骤 3:一个小半径示例

当 \(r=1\) 时,球面整点只有 6 个坐标轴点。北极直接跳到南极时 \(\theta=\pi\),风险等于 \(1\)。如果经过任意一个赤道轴点,则会走两次四分之一圆,每次的归一化角度是 \(1/2\),因此

$$M(1)=2\left(\frac{1}{2}\right)^2=\frac{1}{2}.$$

这个例子说明了为什么把大跳跃拆成若干小跳跃会更优:代价与角度成平方关系,而不是线性关系。

步骤 4:快速算法的镜像闭合公式

优化后的算法先只在上半球的一个对称代表集合上计算从北极出发的最小风险。设某个代表点为

$$p=(x,y,z),\qquad z\ge 0,$$

它关于赤道平面的镜像点为

$$p'=(x,y,-z).$$

若 \(d(p)\) 是从 \(N\) 到 \(p\) 的最小风险,则由对称性可知,从 \(p'\) 到 \(S\) 的最小风险相同。因此一个完整的北到南候选值是

$$2d(p)+w(p,p').$$

又因为

$$p\cdot p'=x^2+y^2-z^2=r^2-2z^2,$$

所以

$$\cos\theta=\frac{r^2-2z^2}{r^2}=1-2\left(\frac{z}{r}\right)^2.$$

令 \(\alpha=\arcsin(z/r)\),则 \(\cos\theta=\cos(2\alpha)\),从而

$$\theta=2\arcsin\left(\frac{z}{r}\right).$$

因此镜像桥边的风险为

$$w(p,p')=\left(\frac{2\arcsin(z/r)}{\pi}\right)^2,$$

快速求解器最终最小化的是

$$\boxed{M(r)=\min_{p}\left(2d(p)+\left(\frac{2\arcsin(z_p/r)}{\pi}\right)^2\right).}$$

步骤 5:约化点集与保守剪枝

generate_reduced_points 并不会保留所有符号变化和坐标置换后的副本,而是先枚举满足

$$0\le x\le y\le z,\qquad x^2+y^2+z^2=r^2$$

的非负有序三元组,再按源码中定义的对称规则输出所需的若干排列。外层条件 \(3x^2<r^2\) 就来自 \(x\le y\le z\)。排序并去重后,这个列表就是快速搜索真正使用的顶点集合。

程序还会预处理

$$a_t=\frac{\arcsin(t/r)}{\pi}\qquad (0\le t\le r).$$

若当前完整北南路径的最好上界为 \(B\),当前单侧风险为 \(d_c\),则代码使用保守剩余预算 \(B-2d_c\)。在计算昂贵的 acos 之前,先用这些 \(a_t\) 的差值构造廉价下界。C++ 与 Java 会同时检查 \(z\)、\(x\)、\(y\) 三个坐标;Python 保留相同的大框架,但只使用基于 \(z\) 的那条下界。

步骤 6:分层细化与校验

快速算法按步长 \(27,9,3,1\) 进行四次搜索。粗采样先给出一个上界,后续更密的采样就可以用这个上界做更多剪枝。因此 minimal_risk_fast 会在同一组约化点上多次调用带剪枝的 Dijkstra。

C++ 源码用两种方式验证正确性。第一,检查基准值

$$M(7)=0.1784943998.$$

第二,对 \(r=31\) 同时运行完整图上的精确 Dijkstra 和约化后的快速算法,并要求两者误差不超过 \(10^{-12}\)。这说明快速算法并不是更换了数学模型,而是对精确模型做了经过验证的优化。

代码如何对应这些公式

C++ 文件同时包含精确校验器和优化求解器。Python 与 Java 则直接实现快速约化版本。三种语言共享同一套数学骨架:边风险公式、\(\arcsin(z/r)/\pi\) 预处理表、代表点生成,以及镜像闭合公式。

实现上的差别主要在调度层。C++ 使用最多 8 个 pthread 工作线程并行计算 15 个半径;Java 用并行 IntStream;Python 按顺序循环。只有 C++ 还保留了 generate_all_pointspole_indicesdijkstra_complete 这组完整图校验工具,用于检查快速算法而不是用于最终大规模计算。

复杂度分析

记 \(V_r=|S_r|\)。完整图上的精确 Dijkstra 需要 \(O(V_r^2)\) 时间和 \(O(V_r)\) 空间,因为每次都要扫描所有尚未确定的顶点。这对 \(r=31\) 这类检查点仍可接受,但对最终求和中的最大半径就太贵了。

快速求解器在最坏情况下依然对约化顶点数呈二次复杂度,因为它仍然在做稠密 relax。真正的加速来自三件事:约化点集远小于完整点集;\(27,9,3,1\) 四层搜索逐步收紧上界;保守剪枝会跳过大量不可能改进答案的边。这三者结合后,\(n=1\) 到 \(15\) 的总和才变得可计算。

参考资料

  1. 题目页面: https://projecteuler.net/problem=353
  2. 大圆距离与球心夹角: Wikipedia — 大圆
  3. Dijkstra 最短路算法: Wikipedia — 迪杰斯特拉算法
  4. 球面整点与三平方定理: Wikipedia — 三平方和定理
  5. 球面余弦定理: Wikipedia — 球面三角学

Краткое описание задачи

Для каждого \(n\) берется \(r=2^n-1\) и рассматривается множество целочисленных точек сферы

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

Северный и южный полюса равны \(N=(0,0,r)\) и \(S=(0,0,-r)\). Переход между двумя узлами имеет риск, равный квадрату центрального угла, нормированного на \(\pi\). Если \(M(r)\) обозначает минимальный суммарный риск пути из \(N\) в \(S\), то требуется вычислить

$$\sum_{n=1}^{15} M(2^n-1).$$

Локальные реализации строят точную модель как плотный взвешенный граф, а затем ускоряют ее за счет симметрий и проверки верхних границ.

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

Шаг 1: Вес ребра через центральный угол

Если \(a,b\in S_r\), то \(\lVert a\rVert=\lVert b\rVert=r\). Центральный угол \(\theta(a,b)\) определяется по формуле

$$\cos\theta(a,b)=\frac{a\cdot b}{r^2},\qquad \theta(a,b)=\arccos\left(\frac{a\cdot b}{r^2}\right).$$

Код делит этот угол на \(\pi\), ограничивает косинус интервалом \([-1,1]\) для устойчивости и задает

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

Именно это и реализуют функции edge_risk и edgeRisk в трех языках.

Шаг 2: Задача превращается в кратчайший путь на полном графе

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

$$P=(p_0,p_1,\dots,p_k),\qquad p_0=N,\quad p_k=S,$$

суммарный риск равен

$$R(P)=\sum_{i=0}^{k-1} w(p_i,p_{i+1}).$$

Следовательно,

$$M(r)=\min_{P:N\leadsto S} R(P).$$

Точная контрольная ветка в C++ решает именно эту задачу алгоритмом Дейкстры на полном множестве точек, построенном функцией generate_all_points.

Шаг 3: Малый пример

При \(r=1\) на сфере есть только 6 осевых точек. Прямой прыжок из \(N\) в \(S\) имеет \(\theta=\pi\), а значит риск \(1\). Если идти через любую экваториальную осевую точку, получаются два поворота на четверть окружности, то есть два риска \((1/2)^2\). Поэтому

$$M(1)=2\left(\frac{1}{2}\right)^2=\frac{1}{2}.$$

Этот пример хорошо показывает, почему несколько коротких переходов могут быть выгоднее одного длинного: квадратичная функция сильно штрафует большие углы.

Шаг 4: Формула зеркального моста в быстрой версии

Оптимизированный алгоритм сначала считает минимальные риски от северного полюса до уменьшенного множества представителей в верхнем полушарии. Пусть

$$p=(x,y,z),\qquad z\ge 0,$$

а его отражение относительно экваториальной плоскости равно

$$p'=(x,y,-z).$$

Если \(d(p)\) есть минимальный риск пути из \(N\) в \(p\), то по симметрии такой же риск нужен для пути из \(p'\) в \(S\). Поэтому полный кандидат север-юг имеет вид

$$2d(p)+w(p,p').$$

Заметим, что

$$p\cdot p'=x^2+y^2-z^2=r^2-2z^2,$$

откуда

$$\cos\theta=\frac{r^2-2z^2}{r^2}=1-2\left(\frac{z}{r}\right)^2.$$

Если \(\alpha=\arcsin(z/r)\), то \(\cos\theta=\cos(2\alpha)\), а значит

$$\theta=2\arcsin\left(\frac{z}{r}\right).$$

Следовательно, риск зеркального мостика равен

$$w(p,p')=\left(\frac{2\arcsin(z/r)}{\pi}\right)^2,$$

и быстрая схема минимизирует

$$\boxed{M(r)=\min_{p}\left(2d(p)+\left(\frac{2\arcsin(z_p/r)}{\pi}\right)^2\right).}$$

Шаг 5: Редуцированное множество и консервативная отсечка

Функция generate_reduced_points не хранит все копии, получаемые перестановками координат и сменой знаков. Вместо этого она перебирает упорядоченные неотрицательные тройки

$$0\le x\le y\le z,\qquad x^2+y^2+z^2=r^2,$$

после чего добавляет только те перестановки, которые нужны выбранной симметрийной модели. Ограничение \(3x^2<r^2\) непосредственно следует из \(x\le y\le z\). После сортировки и удаления дублей именно это множество подается на быструю версию Дейкстры.

Также заранее вычисляются значения

$$a_t=\frac{\arcsin(t/r)}{\pi}\qquad (0\le t\le r).$$

Если текущая лучшая верхняя граница полного пути равна \(B\), а односторонний риск в текущей вершине равен \(d_c\), то код использует консервативный остаточный бюджет \(B-2d_c\). До вычисления дорогого acos проверяются дешевые нижние оценки, построенные из разностей \(a_t\). В C++ и Java такие проверки делаются для координат \(z\), \(x\) и \(y\); в Python сохраняется та же общая схема, но используется только оценка по \(z\).

Шаг 6: Переход от грубого к точному и контрольные точки

Быстрый решатель выполняет четыре прохода со stride \(27,9,3,1\). Грубые проходы дают верхнюю границу, которая ускоряет последующие более плотные проходы. Поэтому minimal_risk_fast несколько раз вызывает одну и ту же процедуру Дейкстры на все более насыщенных подвыборках редуцированного множества.

Файл C++ проверяет метод двумя контрольными условиями. Во-первых, требуется

$$M(7)=0.1784943998.$$

Во-вторых, при \(r=31\) сравниваются точный Дейкстра на полном графе и быстрая редуцированная версия, причем расхождение должно быть меньше \(10^{-12}\). Это важно: быстрый метод не меняет задачу, а лишь ускоряет точную постановку.

Как это отражено в коде

Источник на C++ содержит и точную проверку, и оптимизированный решатель. Python и Java реализуют сразу быструю редуцированную версию. Во всех трех языках совпадают основные математические компоненты: формула риска, таблица \(\arcsin(z/r)/\pi\), генерация представителей и зеркальная формула завершения пути.

Небольшие различия касаются организации вычислений. C++ распределяет 15 радиусов между максимум восемью pthread-потоками, Java использует параллельный IntStream, а Python считает последовательно. Полные функции generate_all_points, pole_indices и dijkstra_complete присутствуют только в C++ и нужны для проверки корректности.

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

Пусть \(V_r=|S_r|\). Точный алгоритм Дейкстры на полном графе требует \(O(V_r^2)\) времени и \(O(V_r)\) памяти, поскольку на каждом шаге просматриваются все еще не обработанные вершины. Для проверок вроде \(r=31\) этого достаточно, но для наибольших радиусов в финальной сумме такой подход слишком дорог.

Быстрая схема в худшем случае все равно квадратична по размеру редуцированного множества, так как релаксации остаются плотными. Практическое ускорение достигается за счет трех факторов: редуцированное множество существенно меньше полного; проходы \(27,9,3,1\) последовательно улучшают верхнюю границу; консервативная отсечка не дает тратить время на ребра, которые уже не могут улучшить ответ.

Ссылки и литература

  1. Страница задачи: https://projecteuler.net/problem=353
  2. Большая окружность и центральный угол: Wikipedia — Дуга большого круга
  3. Алгоритм Дейкстры: Wikipedia — Алгоритм Дейкстры
  4. Теорема о сумме трех квадратов: Wikipedia — Теорема о трёх квадратах
  5. Сферический закон косинусов: Wikipedia — Сферическая тригонометрия

ملخص المسألة

لكل \(n\)، نضع \(r=2^n-1\) وننظر إلى مجموعة النقاط الصحيحة على الكرة

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

يمثل القطبان الشمالي والجنوبي النقطتين \(N=(0,0,r)\) و\(S=(0,0,-r)\). الانتقال بين نقطتين من هذه النقاط له مخاطرة تساوي مربع الزاوية المركزية بعد تطبيعها بالقسمة على \(\pi\). وإذا رمزنا إلى أقل مخاطرة كلية من \(N\) إلى \(S\) بالرمز \(M(r)\)، فالمطلوب هو

$$\sum_{n=1}^{15} M(2^n-1).$$

الملفات المحلية تحل هذه المسألة أولًا كنموذج دقيق على رسم بياني كثيف الأوزان، ثم تسرعه باستغلال التماثل وتقليل عدد الحالات.

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

الخطوة 1: اشتقاق وزن الحافة من هندسة الكرة

إذا كانت \(a,b\in S_r\)، فلدينا \(\lVert a\rVert=\lVert b\rVert=r\). وعليه فإن الزاوية المركزية \(\theta(a,b)\) تحقق

$$\cos\theta(a,b)=\frac{a\cdot b}{r^2},\qquad \theta(a,b)=\arccos\left(\frac{a\cdot b}{r^2}\right).$$

يقوم الكود بتطبيع هذه الزاوية على \(\pi\)، ثم يقيد قيمة جيب التمام في المجال \([-1,1]\) لتفادي أخطاء الفاصلة العائمة، ويعرّف

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

وهذه هي بالضبط الصيغة الموجودة في edge_risk وedgeRisk في ملفات C++ وPython وJava.

الخطوة 2: تحويل المسألة إلى أقصر مسار على رسم كامل

بما أن التنفيذ يسمح بالانتقال المباشر بين أي نقطتين من \(S_r\)، فإننا نحصل على رسم بياني كامل موزون. إذا كان

$$P=(p_0,p_1,\dots,p_k),\qquad p_0=N,\quad p_k=S,$$

مسارًا من الشمال إلى الجنوب، فإن مخاطرة هذا المسار هي

$$R(P)=\sum_{i=0}^{k-1} w(p_i,p_{i+1}).$$

ومن ثم

$$M(r)=\min_{P:N\leadsto S} R(P).$$

إذًا جوهر المسألة هو مسألة أقصر مسار. ولهذا يستخدم مسار التحقق الدقيق في C++ خوارزمية Dijkstra على جميع النقاط المولدة بواسطة generate_all_points.

الخطوة 3: مثال صغير يوضح الفكرة

عندما \(r=1\)، لا توجد على الكرة إلا ست نقاط محورية. القفزة المباشرة من الشمال إلى الجنوب تعطي \(\theta=\pi\)، ومن ثم مخاطرة مقدارها \(1\). أما إذا مررنا عبر نقطة على خط الاستواء، فسننفذ ربع دورة مرتين، أي مخاطرتين كل منهما \((1/2)^2\)، ولذلك

$$M(1)=2\left(\frac{1}{2}\right)^2=\frac{1}{2}.$$

هذا المثال يبين لماذا قد يكون تفكيك القفزة الكبيرة إلى قفزات أصغر أفضل: لأن دالة الخطر تربيعية في الزاوية وليست خطية.

الخطوة 4: صيغة الجسر المرآتي في الحل السريع

الخوارزمية المحسنة تحسب أولًا أقل المخاطر من القطب الشمالي إلى مجموعة مختزلة من النقاط الممثلة في نصف الكرة العلوي. لتكن

$$p=(x,y,z),\qquad z\ge 0,$$

ولنأخذ انعكاسها حول مستوى خط الاستواء

$$p'=(x,y,-z).$$

إذا كانت \(d(p)\) هي أقل مخاطرة من \(N\) إلى \(p\)، فإن التناظر يعطينا القيمة نفسها من \(p'\) إلى \(S\). لذلك فإن أي مرشح كامل من الشمال إلى الجنوب يساوي

$$2d(p)+w(p,p').$$

ولدينا أيضًا

$$p\cdot p'=x^2+y^2-z^2=r^2-2z^2,$$

ومنها

$$\cos\theta=\frac{r^2-2z^2}{r^2}=1-2\left(\frac{z}{r}\right)^2.$$

إذا وضعنا \(\alpha=\arcsin(z/r)\)، فإن \(\cos\theta=\cos(2\alpha)\)، وبالتالي

$$\theta=2\arcsin\left(\frac{z}{r}\right).$$

وعليه فإن مخاطرة الحافة المرآتية تصبح

$$w(p,p')=\left(\frac{2\arcsin(z/r)}{\pi}\right)^2,$$

والكمية التي تصغرها الخوارزمية السريعة هي

$$\boxed{M(r)=\min_{p}\left(2d(p)+\left(\frac{2\arcsin(z_p/r)}{\pi}\right)^2\right).}$$

الخطوة 5: مجموعة النقاط المختزلة والتقليم المحافظ

الدالة generate_reduced_points لا تحتفظ بكل النسخ الناتجة عن تبديل الإشارات وترتيب الإحداثيات. بدلًا من ذلك، تمر على الثلاثيات غير السالبة المرتبة

$$0\le x\le y\le z,\qquad x^2+y^2+z^2=r^2,$$

ثم تصدر فقط التبديلات اللازمة لنموذج التماثل المستخدم في الحل. أما القيد \(3x^2<r^2\) فيأتي مباشرة من الشرط \(x\le y\le z\). وبعد الفرز وحذف التكرارات نحصل على فضاء الحالات الذي تعمل عليه النسخة السريعة.

كذلك يسبق الحساب الجدول

$$a_t=\frac{\arcsin(t/r)}{\pi}\qquad (0\le t\le r).$$

إذا كانت أفضل قيمة كلية معروفة حاليًا هي \(B\)، وكانت المخاطرة الأحادية الحالية هي \(d_c\)، فإن الكود يستخدم الميزانية المحافظة المتبقية \(B-2d_c\). قبل تنفيذ استدعاء acos المكلف، يختبر حدودًا دنيا رخيصة مبنية على فروق القيم \(a_t\). في C++ وJava يُطبَّق هذا الاختبار على الإحداثيات \(z\) و\(x\) و\(y\)، أما Python فيحتفظ بالبنية نفسها لكنه يستخدم فقط الحد المعتمد على \(z\).

الخطوة 6: البحث من الخشن إلى الدقيق ونقاط التحقق

الحل السريع يعمل على أربع مراحل بخطوات \(27,9,3,1\). المرحلة الخشنة تعطي حدًا علويًا، وهذا الحد يجعل المرحلة الأدق أقل تكلفة. لهذا السبب تستدعي الدالة minimal_risk_fast خوارزمية Dijkstra المقلمة عدة مرات على عينات تزداد كثافتها تدريجيًا.

ملف C++ يتحقق من صحة الفكرة بطريقتين. أولًا، يفرض القيمة المرجعية

$$M(7)=0.1784943998.$$

ثانيًا، عند \(r=31\) يشغل كلًا من الحل الدقيق على الرسم الكامل والحل السريع المختزل، ويطلب أن يتطابقا حتى \(10^{-12}\). وهذا مهم لأنه يثبت أن النسخة السريعة ليست نموذجًا مختلفًا، بل تسريعًا محسوبًا للنموذج الدقيق نفسه.

كيف يعمل الكود

المصدر في C++ يحتوي على أداة التحقق الدقيقة وعلى الحل المحسن معًا. أما ملفا Python وJava فيطبقان مباشرة النسخة السريعة المختزلة. ومع ذلك فالعناصر الرياضية الأساسية واحدة في اللغات الثلاث: صيغة المخاطرة، وجدول \(\arcsin(z/r)/\pi\)، وتوليد نقاط التمثيل، وصيغة الإغلاق بالمرآة.

الاختلافات التنفيذية طفيفة. C++ يوزع أنصاف الأقطار الخمسة عشر على ما يصل إلى ثمانية خيوط pthread، وJava يستخدم IntStream متوازيًا، وPython ينفذ حلقة تسلسلية. أما الدوال الدقيقة generate_all_points وpole_indices وdijkstra_complete فهي موجودة فقط في C++ وتستخدم للتحقق من الصحة، لا للحساب النهائي الكبير.

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

إذا رمزنا إلى عدد نقاط الكرة الصحيحة بـ \(V_r=|S_r|\)، فإن Dijkstra الدقيق على الرسم الكامل يستهلك \(O(V_r^2)\) زمنًا و\(O(V_r)\) ذاكرة، لأن كل خطوة تسترخي جميع الرؤوس المتبقية. هذا مقبول للحالات المرجعية مثل \(r=31\)، لكنه غير عملي لأكبر أنصاف الأقطار في المجموع النهائي.

أما الحل السريع فيبقى تربيعيًا في أسوأ حالة بالنسبة إلى عدد النقاط المختزلة، لأنه ما زال يتعامل مع استرخاءات كثيفة. لكن الأداء العملي يتحسن بشدة لثلاثة أسباب: المجموعة المختزلة أصغر بكثير من المجموعة الكاملة، والمراحل \(27,9,3,1\) تعطي حدودًا علوية أفضل فأفضل، والتقليم المحافظ يستبعد عددًا كبيرًا من الحواف قبل تقييمها كاملًا.

مراجع

  1. صفحة المسألة: https://projecteuler.net/problem=353
  2. المسافة على الدائرة العظمى والزاوية المركزية: Wikipedia — دائرة عظمى
  3. خوارزمية Dijkstra: Wikipedia — خوارزمية ديكسترا
  4. تمثيل الأعداد كمجموع ثلاثة مربعات: Wikipedia — Sum of three squares theorem
  5. قانون جيب التمام الكروي: Wikipedia — مثلثات كروية