Problem Summary

Define

$$S(r)=\sum_{x^2+y^2+z^2=r^2}\left(|x|+|y|+|z|\right),$$

where the sum runs over all integer lattice points on the sphere of radius \(r\). For Project Euler 360 the target radius is \(r=10^{10}=2^{10}5^{10}\). A direct enumeration of all triples \((x,y,z)\) is far too expensive, so the solution rewrites the three-dimensional problem as a one-dimensional scan combined with the classical sum-of-two-squares function.

Mathematical Approach

The repository solutions all use the same two ideas:

$$\text{(i) symmetry reduces the sum to }r_2(r^2-x^2), \qquad \text{(ii) }10^{10}=2^{10}5^{10}\text{ lets us separate the powers of }2\text{ and }5.$$

Step 1: Reduce the Three-Variable Sum to One Coordinate

The three coordinates play identical roles, so the total contribution of \(|x|\), \(|y|\), and \(|z|\) is the same. It is therefore enough to compute the \(|x|\)-part and multiply by \(3\).

For a fixed positive \(x\), every pair \((y,z)\in\mathbb Z^2\) satisfying

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

appears twice in the sphere sum, once with coordinate \(x\) and once with coordinate \(-x\). Hence the total contribution coming from the \(x\)-coordinate is

$$2x\,r_2(r^2-x^2),$$

where \(r_2(n)\) denotes the number of ordered integer pairs \((u,v)\) with \(u^2+v^2=n\). Multiplying by \(3\) for the three symmetric coordinates gives

$$\boxed{S(r)=6\sum_{x=1}^{r} x\,r_2(r^2-x^2).}$$

The term \(x=0\) does not matter because it would be multiplied by \(x\). The endpoint \(x=r\) is special: then \(r^2-x^2=0\), so \(r_2(0)=1\) from the single pair \((0,0)\).

Step 2: Use the Sum-of-Two-Squares Function

For each \(x\), write

$$r^2-x^2=(r-x)(r+x)=ab,$$

with \(a=r-x\) and \(b=r+x\). The problem is now to evaluate \(r_2(ab)\).

If

$$n=2^\alpha \prod_{p\equiv 1 \pmod{4}} p^{e_p}\prod_{q\equiv 3 \pmod{4}} q^{f_q},$$

then the classical formula for the ordered signed representation count is

$$r_2(n)= \begin{cases} 0, & \text{if some } f_q \text{ is odd},\\ 4\displaystyle\prod_{p\equiv 1 \pmod{4}}(e_p+1), & \text{otherwise}. \end{cases}$$

So primes congruent to \(3 \pmod{4}\) are only a parity test: any odd exponent kills the term. Primes congruent to \(1 \pmod{4}\) contribute multiplicative factors \(e_p+1\). The exponent of \(2\) does not affect the value of \(r_2(n)\), so the implementation does not need a separate table for powers of \(2\).

Step 3: Why the Factor Tables for \(r=5^b\) Are Enough

The code first handles radii of the form \(r=5^b\). For \(x\in\{1,\dots,r-1\}\), set

$$a=r-x,\qquad b=r+x.$$

Any common divisor of \(a\) and \(b\) divides both \(a+b=2r\) and \(b-a=2x\), hence

$$\gcd(a,b)\mid 2r=2\cdot 5^b.$$

This simple fact is the key simplification. It implies that any odd prime \(p\neq 5\) can divide at most one of \(a\) and \(b\). Therefore:

\(1.\) A prime \(q\equiv 3 \pmod{4}\) has odd exponent in \(ab\) exactly when it has odd exponent in \(a\) or in \(b\). It is enough to know whether each factor is “bad”.

\(2.\) A prime \(p\equiv 1 \pmod{4}\), \(p\neq 5\), contributes independently from \(a\) and \(b\), because it cannot appear in both factors.

\(3.\) The prime \(5\) must be treated separately, because it is the only prime congruent to \(1 \pmod{4}\) that may appear in both \(a\) and \(b\).

That is why the solutions precompute, for every \(m\le 2r\):

$$\texttt{bad}(m)= \begin{cases} 1,& \text{if some } q\equiv 3 \pmod{4} \text{ has odd exponent in }m,\\ 0,& \text{otherwise}, \end{cases}$$

$$\nu_5(m)=v_5(m),$$

$$g(m)=\prod_{\substack{p\equiv 1 \pmod{4}\\ p\neq 5}}(e_p(m)+1).$$

If either \(a\) or \(b\) is bad, then \(r_2(ab)=0\). Otherwise the exponent of \(5\) in \(ab\) is \(\nu_5(a)+\nu_5(b)\), and the remaining \(1 \pmod{4}\) primes factor cleanly, so

$$\boxed{r_2(ab)=4\,g(a)g(b)\bigl(\nu_5(a)+\nu_5(b)+1\bigr).}$$

This is exactly the expression used in the C++, Python, and Java implementations.

Step 4: Recover the Radius \(2^a5^b\)

Now let the radius be \(R=2^a5^b\). Because \(R^2\) is divisible by \(4\), any solution of

$$x^2+y^2+z^2=R^2$$

must have all three coordinates even: modulo \(4\), each square is \(0\) or \(1\), and the only way three such residues can sum to \(0\) is \(0+0+0\). Dividing by \(2\) and repeating this argument \(a\) times shows that every solution has the form

$$x=2^a x',\qquad y=2^a y',\qquad z=2^a z',$$

with

$$x'^2+y'^2+z'^2=5^{2b}.$$

This is a bijection between lattice points on the sphere of radius \(2^a5^b\) and those on the sphere of radius \(5^b\). The absolute-value sum scales by the same factor \(2^a\), so

$$\boxed{S(2^a5^b)=2^aS(5^b).}$$

That explains the final left shift in the code after the \(5^b\) case has been computed.

Step 5: Worked Example with \(r=5\)

The small checkpoint \(r=5\) is useful because the whole derivation can be checked by hand:

$$S(5)=6\sum_{x=1}^{5}x\,r_2(25-x^2).$$

Evaluate each term:

\(x=1\): \(25-1=24=2^3\cdot 3\), so a prime \(3 \pmod{4}\) has odd exponent and \(r_2(24)=0\).

\(x=2\): \(25-4=21=3\cdot 7\), again impossible, so \(r_2(21)=0\).

\(x=3\): \(25-9=16\), so \(r_2(16)=4\) from \((\pm4,0)\) and \((0,\pm4)\).

\(x=4\): \(25-16=9\), so \(r_2(9)=4\) from \((\pm3,0)\) and \((0,\pm3)\).

\(x=5\): \(25-25=0\), so \(r_2(0)=1\).

Therefore

$$S(5)=6\cdot 3\cdot 4 + 6\cdot 4\cdot 4 + 6\cdot 5 = 72+96+30=198.$$

The scaling law then gives

$$S(10)=2S(5)=396,$$

which matches the checkpoint used by the implementation.

How the Code Works

The C++ solution builds an SPF table up to \(2r\), where \(r=5^{10}\), with build_spf. Then build_factor_tables factors every integer \(1\le m\le 2r\) and stores the three quantities described above: the bad flag, the exponent of \(5\), and the multiplicative factor from primes congruent to \(1 \pmod{4}\) other than \(5\).

The main summation loops over \(x=1,\dots,r-1\), computes \(a=r-x\) and \(b=r+x\), rejects the pair immediately if either factor is bad, and otherwise evaluates

$$\texttt{ways}=4\,g(a)g(b)\bigl(\nu_5(a)+\nu_5(b)+1\bigr).$$

It then adds

$$6x\cdot \texttt{ways}$$

to the total. After the loop it adds the endpoint contribution \(6r\) for \(x=r\), and finally multiplies by \(2^{10}\) via a left shift.

The Python and Java files implement the same arithmetic formula. The Java version stores the final answer in BigInteger. The C++ version uses unsigned __int128 and parallelizes the main \(x\)-range across worker threads. It also contains explicit checkpoints: \(S(5)=198\), \(S(10)=396\), \(S(25)=5766\), and a brute-force verification \(S(45)=34518\).

Complexity Analysis

Let \(r=5^{10}\) and \(\texttt{limit}=2r\). Building the SPF sieve costs \(O(\texttt{limit}\log\log \texttt{limit})\) time and \(O(\texttt{limit})\) memory. Filling the factor tables by repeatedly stripping SPF factors is near-linear on average, again about \(O(\texttt{limit}\log\log \texttt{limit})\) time. The final scan over \(x=1,\dots,r-1\) is \(O(r)\) with only constant-time table lookups per iteration. The overall memory usage is \(O(r)\), dominated by the precomputed arrays. The C++ threading changes wall-clock time, but not the asymptotic complexity.

Footnotes and References

  1. Problem page: https://projecteuler.net/problem=360
  2. Sum of two squares theorem: Wikipedia — Sum of two squares theorem
  3. Sum of two squares function: Wikipedia — Sum of two squares function
  4. Sieve of Eratosthenes and SPF tables: Wikipedia — Sieve of Eratosthenes

Problemzusammenfassung

Definiert ist

$$S(r)=\sum_{x^2+y^2+z^2=r^2}\left(|x|+|y|+|z|\right),$$

wobei über alle ganzzahligen Gitterpunkte auf der Kugeloberfläche mit Radius \(r\) summiert wird. Für Problem 360 ist \(r=10^{10}=2^{10}5^{10}\). Eine direkte Enumeration aller Tripel \((x,y,z)\) ist völlig unpraktisch, daher reduziert die Lösung das dreidimensionale Problem auf eine eindimensionale Summe mit der klassischen Darstellungsfunktion als Summe zweier Quadrate.

Mathematischer Ansatz

Alle Lösungsdateien im Repository verwenden dieselbe Struktur:

$$\text{(i) Symmetrie reduziert auf }r_2(r^2-x^2), \qquad \text{(ii) }10^{10}=2^{10}5^{10}\text{ trennt Zweier- und Fünferpotenzen.}$$

Schritt 1: Die Dreivariablensumme auf eine Koordinate reduzieren

Die drei Koordinaten spielen exakt dieselbe Rolle. Daher ist der Gesamtbeitrag von \(|x|\), \(|y|\) und \(|z|\) jeweils gleich, und es genügt, den \(|x|\)-Anteil zu berechnen und anschließend mit \(3\) zu multiplizieren.

Für ein festes positives \(x\) liefert jede Lösung \((y,z)\in\mathbb Z^2\) von

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

zwei Punkte auf der Kugel, nämlich mit \(x\) und mit \(-x\). Der Beitrag der \(x\)-Koordinate ist also

$$2x\,r_2(r^2-x^2),$$

wobei \(r_2(n)\) die Anzahl geordneter ganzzahliger Paare \((u,v)\) mit \(u^2+v^2=n\) bezeichnet. Durch Multiplikation mit \(3\) erhält man

$$\boxed{S(r)=6\sum_{x=1}^{r} x\,r_2(r^2-x^2).}$$

Der Fall \(x=0\) entfällt wegen des Faktors \(x\). Der Randfall \(x=r\) ist besonders: Dann ist \(r^2-x^2=0\), also \(r_2(0)=1\) wegen des einzigen Paares \((0,0)\).

Schritt 2: Die Funktion für Summen zweier Quadrate verwenden

Für jedes \(x\) schreibe

$$r^2-x^2=(r-x)(r+x)=ab,$$

mit \(a=r-x\) und \(b=r+x\). Damit muss nur noch \(r_2(ab)\) effizient bestimmt werden.

Ist

$$n=2^\alpha \prod_{p\equiv 1 \pmod{4}} p^{e_p}\prod_{q\equiv 3 \pmod{4}} q^{f_q},$$

dann gilt für die Anzahl geordneter Darstellungen mit Vorzeichen

$$r_2(n)= \begin{cases} 0, & \text{falls ein } f_q \text{ ungerade ist},\\ 4\displaystyle\prod_{p\equiv 1 \pmod{4}}(e_p+1), & \text{sonst}. \end{cases}$$

Primzahlen \(3 \pmod{4}\) sind also nur ein Paritätstest: Ein ungerader Exponent macht den Term unmöglich. Primzahlen \(1 \pmod{4}\) liefern die Faktoren \(e_p+1\). Die Zweierpotenz beeinflusst den Wert von \(r_2(n)\) nicht; deshalb braucht der Code keine eigene Tabelle für den Exponenten von \(2\).

Schritt 3: Warum bei \(r=5^b\) genau diese Tabellen genügen

Der Code behandelt zunächst Radien der Form \(r=5^b\). Für \(x\in\{1,\dots,r-1\}\) setze

$$a=r-x,\qquad b=r+x.$$

Jeder gemeinsame Teiler von \(a\) und \(b\) teilt sowohl \(a+b=2r\) als auch \(b-a=2x\). Also gilt

$$\gcd(a,b)\mid 2r=2\cdot 5^b.$$

Genau das macht die Faktorisierung handhabbar. Daraus folgt nämlich, dass jede ungerade Primzahl \(p\neq 5\) höchstens einen der beiden Faktoren \(a\) oder \(b\) teilen kann. Insbesondere:

\(1.\) Für \(q\equiv 3 \pmod{4}\) ist der Exponent in \(ab\) genau dann ungerade, wenn er in \(a\) oder in \(b\) ungerade ist. Daher genügt eine boolesche „bad“-Information pro Faktor.

\(2.\) Für \(p\equiv 1 \pmod{4}\), \(p\neq 5\), faktorisiert der Beitrag unabhängig über \(a\) und \(b\), weil ein solches \(p\) nicht beide Faktoren zugleich teilen kann.

\(3.\) Die Primzahl \(5\) muss separat behandelt werden, denn sie ist die einzige Primzahl \(1 \pmod{4}\), die sowohl in \(a\) als auch in \(b\) auftreten kann.

Deshalb berechnen die Lösungen für jedes \(m\le 2r\) vorab:

$$\texttt{bad}(m)= \begin{cases} 1,& \text{falls ein } q\equiv 3 \pmod{4} \text{ in }m\text{ ungeraden Exponenten hat},\\ 0,& \text{sonst}, \end{cases}$$

$$\nu_5(m)=v_5(m),$$

$$g(m)=\prod_{\substack{p\equiv 1 \pmod{4}\\ p\neq 5}}(e_p(m)+1).$$

Ist einer der beiden Faktoren bad, dann ist \(r_2(ab)=0\). Andernfalls hat \(5\) in \(ab\) den Exponenten \(\nu_5(a)+\nu_5(b)\), und die übrigen Primzahlen \(1 \pmod{4}\) zerfallen sauber, also

$$\boxed{r_2(ab)=4\,g(a)g(b)\bigl(\nu_5(a)+\nu_5(b)+1\bigr).}$$

Genau diese Formel steht in den C++-, Python- und Java-Lösungen.

Schritt 4: Den Radius \(2^a5^b\) zurückgewinnen

Sei nun der Radius \(R=2^a5^b\). Weil \(R^2\) durch \(4\) teilbar ist, müssen in jeder Lösung von

$$x^2+y^2+z^2=R^2$$

alle drei Koordinaten gerade sein: Modulo \(4\) ist jedes Quadrat entweder \(0\) oder \(1\), und \(0\) kann nur als \(0+0+0\) entstehen. Teilt man durch \(2\) und wiederholt das Argument \(a\)-mal, so hat jede Lösung die Form

$$x=2^a x',\qquad y=2^a y',\qquad z=2^a z',$$

mit

$$x'^2+y'^2+z'^2=5^{2b}.$$

Das ist eine Bijektion zwischen den Gitterpunkten auf der Kugel mit Radius \(2^a5^b\) und denen mit Radius \(5^b\). Der Ausdruck \(|x|+|y|+|z|\) skaliert dabei ebenfalls mit \(2^a\), also

$$\boxed{S(2^a5^b)=2^aS(5^b).}$$

Damit ist auch klar, warum der Code nach der Berechnung für \(5^b\) am Ende nur noch nach links schiebt.

Schritt 5: Durchgerechnetes Beispiel \(r=5\)

Der kleine Kontrollwert \(r=5\) lässt sich vollständig von Hand ausrechnen:

$$S(5)=6\sum_{x=1}^{5}x\,r_2(25-x^2).$$

Die einzelnen Terme sind:

\(x=1\): \(25-1=24=2^3\cdot 3\), also ein ungerader Exponent einer Primzahl \(3 \pmod{4}\), daher \(r_2(24)=0\).

\(x=2\): \(25-4=21=3\cdot 7\), ebenfalls unmöglich, also \(r_2(21)=0\).

\(x=3\): \(25-9=16\), daher \(r_2(16)=4\) mit \((\pm4,0)\) und \((0,\pm4)\).

\(x=4\): \(25-16=9\), daher \(r_2(9)=4\) mit \((\pm3,0)\) und \((0,\pm3)\).

\(x=5\): \(25-25=0\), also \(r_2(0)=1\).

Damit folgt

$$S(5)=6\cdot 3\cdot 4 + 6\cdot 4\cdot 4 + 6\cdot 5 = 72+96+30=198.$$

Aus dem Skalierungsgesetz folgt sofort

$$S(10)=2S(5)=396,$$

genau wie im Checkpoint der Implementierung.

Wie der Code arbeitet

Die C++-Lösung baut mit build_spf zunächst eine Tabelle der kleinsten Primfaktoren bis \(2r\) auf, wobei \(r=5^{10}\) ist. Danach faktorisiert build_factor_tables alle Zahlen \(1\le m\le 2r\) und speichert die drei benötigten Größen: das bad-Flag, den Exponenten von \(5\) und den multiplikativen Faktor der Primzahlen \(1 \pmod{4}\) außer \(5\).

Die Hauptsumme läuft über \(x=1,\dots,r-1\), berechnet \(a=r-x\) und \(b=r+x\), verwirft bad-Fälle sofort und setzt sonst

$$\texttt{ways}=4\,g(a)g(b)\bigl(\nu_5(a)+\nu_5(b)+1\bigr).$$

Zum Ergebnis addiert der Code dann

$$6x\cdot \texttt{ways}.$$

Nach der Schleife kommt noch der Randbeitrag \(6r\) für \(x=r\) hinzu, und ganz am Ende erfolgt die Multiplikation mit \(2^{10}\) durch eine Linksschiebung.

Die Python- und Java-Dateien verwenden dieselbe Formel. Die Java-Version speichert die Endsumme in BigInteger. Die C++-Version nutzt unsigned __int128 und parallelisiert den \(x\)-Bereich über mehrere Threads. Zusätzlich enthält sie Kontrollwerte: \(S(5)=198\), \(S(10)=396\), \(S(25)=5766\) sowie eine brute-force-Prüfung \(S(45)=34518\).

Komplexitätsanalyse

Sei \(r=5^{10}\) und \(\texttt{limit}=2r\). Das SPF-Sieb benötigt \(O(\texttt{limit}\log\log \texttt{limit})\) Zeit und \(O(\texttt{limit})\) Speicher. Das Füllen der Faktortabellen durch wiederholtes Entfernen der SPF-Faktoren ist im Mittel ebenfalls nahezu linear, also etwa \(O(\texttt{limit}\log\log \texttt{limit})\). Der abschließende Durchlauf über \(x=1,\dots,r-1\) ist \(O(r)\) mit konstanten Tabellenzugriffen pro Schritt. Insgesamt dominiert \(O(r)\)-Speicher für die vorberechneten Arrays. Die Parallelisierung in C++ verbessert nur die Laufzeit in der Praxis, nicht die asymptotische Ordnung.

Fußnoten und Quellen

  1. Problemseite: https://projecteuler.net/problem=360
  2. Satz über Summen zweier Quadrate: Wikipedia — Fermatscher Zwei-Quadrate-Satz
  3. Darstellungsfunktion \(r_2(n)\): Wikipedia — Sum of two squares function
  4. Sieb des Eratosthenes und SPF-Tabellen: Wikipedia — Sieb des Eratosthenes

Problem Özeti

Problemde

$$S(r)=\sum_{x^2+y^2+z^2=r^2}\left(|x|+|y|+|z|\right)$$

tanımlanır; toplam, yarıçapı \(r\) olan küre üzerindeki tüm tam sayı kafes noktaları üzerinde alınır. Euler 360 için hedef \(r=10^{10}=2^{10}5^{10}\) değeridir. Tüm \((x,y,z)\) üçlülerini doğrudan taramak mümkün olmadığından, çözüm üç değişkenli toplamı tek değişkenli bir taramaya ve iki kare toplamı sayım fonksiyonuna indirger.

Matematiksel Yaklaşım

Depodaki tüm çözüm dosyaları aynı iki fikri kullanır:

$$\text{(i) simetri toplamı }r_2(r^2-x^2)\text{ biçimine indirger,} \qquad \text{(ii) }10^{10}=2^{10}5^{10}\text{ ayrışımı }2\text{ ve }5\text{ kuvvetlerini ayırır.}$$

Adım 1: Üç değişkenli toplamı tek koordinata indirgeme

\(|x|\), \(|y|\) ve \(|z|\) terimleri tamamen simetriktir. Bu yüzden önce yalnızca \(|x|\) katkısını hesaplayıp sonucu \(3\) ile çarpmak yeterlidir.

Sabit bir pozitif \(x\) için

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

denklemini sağlayan her \((y,z)\in\mathbb Z^2\) çifti, küre üzerinde biri \(x\), diğeri \(-x\) olan iki nokta üretir. Dolayısıyla \(x\)-koordinatının toplam katkısı

$$2x\,r_2(r^2-x^2)$$

olur. Burada \(r_2(n)\), \(u^2+v^2=n\) denklemini sağlayan sıralı tam sayı çiftlerinin sayısıdır. Üç koordinat için çarpınca

$$\boxed{S(r)=6\sum_{x=1}^{r} x\,r_2(r^2-x^2)}$$

elde edilir. \(x=0\) terimi zaten \(x\) ile çarpıldığı için katkı vermez. \(x=r\) sınır durumunda ise \(r^2-x^2=0\) olur ve tek çift \((0,0)\) olduğundan \(r_2(0)=1\) gelir.

Adım 2: İki kare toplamı fonksiyonunu kullanma

Her \(x\) için

$$r^2-x^2=(r-x)(r+x)=ab$$

yazarız; burada \(a=r-x\), \(b=r+x\). Böylece mesele \(r_2(ab)\) değerini bulmaya indirgenir.

Eğer

$$n=2^\alpha \prod_{p\equiv 1 \pmod{4}} p^{e_p}\prod_{q\equiv 3 \pmod{4}} q^{f_q}$$

şeklinde asal çarpanlara ayrılıyorsa, sıralı ve işaretli gösterim sayısı için klasik formül

$$r_2(n)= \begin{cases} 0, & \text{eğer bir } f_q \text{ tek ise},\\ 4\displaystyle\prod_{p\equiv 1 \pmod{4}}(e_p+1), & \text{aksi halde} \end{cases}$$

şeklindedir. Yani \(3 \pmod{4}\) sınıfındaki asallar yalnızca bir parite testi olarak davranır: tek üs gelirse terim yok olur. \(1 \pmod{4}\) sınıfındaki asallar ise \((e_p+1)\) çarpanları getirir. \(2\)'nin üssü \(r_2(n)\) değerini değiştirmediği için kodun ayrıca \(2\)-üssü tablosu tutmasına gerek yoktur.

Adım 3: \(r=5^b\) için neden bu tablolar yeterli?

Kod önce \(r=5^b\) biçimindeki yarıçapları çözer. \(x\in\{1,\dots,r-1\}\) için

$$a=r-x,\qquad b=r+x$$

tanımlanır. \(a\) ve \(b\)'nin her ortak böleni hem \(a+b=2r\)'yi hem de \(b-a=2x\)'i böler; dolayısıyla

$$\gcd(a,b)\mid 2r=2\cdot 5^b.$$

Asıl sadeleşme burada gelir. Bu sonuç, tek her asal \(p\neq 5\)'in en fazla \(a\) veya \(b\)'den birini bölebileceğini söyler. Buradan şu üç sonuç çıkar:

\(1.\) \(q\equiv 3 \pmod{4}\) tipindeki bir asalın \(ab\) içindeki üssü, ancak \(a\) veya \(b\) içindeki üssü tekse tek olur. Bu yüzden her faktör için yalnızca “bad” bilgisi yeterlidir.

\(2.\) \(p\equiv 1 \pmod{4}\), \(p\neq 5\), türündeki asallar için katkı \(a\) ve \(b\) üzerinden bağımsızca çarpanlara ayrılır; çünkü böyle bir asal iki faktörde birden bulunamaz.

\(3.\) \(5\) ayrıca ele alınmalıdır; çünkü \(1 \pmod{4}\) sınıfında olup hem \(a\) hem \(b\)'yi bölebilen tek asal odur.

Bu yüzden çözümler, \(m\le 2r\) için şu üç tabloyu hazırlar:

$$\texttt{bad}(m)= \begin{cases} 1,& \text{eğer }m\text{ içinde } q\equiv 3 \pmod{4}\text{ türünden bir asalın üssü tekse},\\ 0,& \text{aksi halde}, \end{cases}$$

$$\nu_5(m)=v_5(m),$$

$$g(m)=\prod_{\substack{p\equiv 1 \pmod{4}\\ p\neq 5}}(e_p(m)+1).$$

Eğer \(a\) ya da \(b\) bad ise \(r_2(ab)=0\) olur. Aksi halde \(ab\) içindeki \(5\) üssü \(\nu_5(a)+\nu_5(b)\) olur ve diğer \(1 \pmod{4}\) asalları temizce ayrışır:

$$\boxed{r_2(ab)=4\,g(a)g(b)\bigl(\nu_5(a)+\nu_5(b)+1\bigr).}$$

C++, Python ve Java kodlarının kullandığı ifade tam olarak budur.

Adım 4: \(2^a5^b\) yarıçapına geri dönme

Şimdi yarıçap \(R=2^a5^b\) olsun. \(R^2\) sayısı \(4\)'e bölündüğü için

$$x^2+y^2+z^2=R^2$$

denklemini sağlayan her çözümde tüm koordinatlar çift olmalıdır. Çünkü mod \(4\)'te her kare yalnızca \(0\) veya \(1\) olabilir ve üç karenin toplamı ancak \(0+0+0\) ile \(0\) verebilir. Bu argümanı \(a\) kez tekrarlayınca her çözümün

$$x=2^a x',\qquad y=2^a y',\qquad z=2^a z'$$

şeklinde olduğu ve

$$x'^2+y'^2+z'^2=5^{2b}$$

eşitliğini sağladığı görülür. Bu, \(2^a5^b\) yarıçaplı küre üzerindeki noktalardan \(5^b\) yarıçaplı küre üzerindekilere bir birebir eşlemedir. Ayrıca mutlak değerlerin toplamı da \(2^a\) ile ölçeklenir. Sonuç:

$$\boxed{S(2^a5^b)=2^aS(5^b).}$$

Kodun son adımda sola kaydırma yapmasının nedeni budur.

Adım 5: \(r=5\) için çözümlü örnek

Küçük kontrol değeri \(r=5\), tüm türetimi elde doğrulamaya uygundur:

$$S(5)=6\sum_{x=1}^{5}x\,r_2(25-x^2).$$

Terimleri tek tek inceleyelim:

\(x=1\): \(25-1=24=2^3\cdot 3\), yani \(3 \pmod{4}\) sınıfından bir asal tek üsle geliyor; bu yüzden \(r_2(24)=0\).

\(x=2\): \(25-4=21=3\cdot 7\), yine temsil yok; \(r_2(21)=0\).

\(x=3\): \(25-9=16\), dolayısıyla \(r_2(16)=4\); çözümler \((\pm4,0)\) ve \((0,\pm4)\).

\(x=4\): \(25-16=9\), dolayısıyla \(r_2(9)=4\); çözümler \((\pm3,0)\) ve \((0,\pm3)\).

\(x=5\): \(25-25=0\), yani \(r_2(0)=1\).

Böylece

$$S(5)=6\cdot 3\cdot 4 + 6\cdot 4\cdot 4 + 6\cdot 5 = 72+96+30=198$$

elde edilir. Ölçekleme kuralı da hemen

$$S(10)=2S(5)=396$$

sonucunu verir; bu da implementasyondaki kontrol değeriyle aynıdır.

Kodun Çalışma Mantığı

C++ çözümü önce build_spf ile \(2r\)'ye kadar en küçük asal bölen tablosunu kurar; burada \(r=5^{10}\). Sonra build_factor_tables, \(1\le m\le 2r\) aralığındaki her sayıyı çarpanlarına ayırarak üç tabloyu doldurur: bad bayrağı, \(5\) üssü ve \(5\) dışındaki \(1 \pmod{4}\) asallardan gelen çarpımsal katsayı.

Ana toplam döngüsü \(x=1,\dots,r-1\) için \(a=r-x\) ve \(b=r+x\) değerlerini hesaplar, bad durumlarını anında atlar ve aksi halde

$$\texttt{ways}=4\,g(a)g(b)\bigl(\nu_5(a)+\nu_5(b)+1\bigr)$$

ifadesini kullanır. Toplama eklenen terim

$$6x\cdot \texttt{ways}$$

olur. Döngü bittikten sonra \(x=r\) için \(6r\) eklenir ve son olarak \(2^{10}\) çarpanı sola kaydırma ile uygulanır.

Python ve Java sürümleri aynı aritmetik formülü kullanır. Java tarafı nihai sonucu BigInteger içinde tutar. C++ tarafı unsigned __int128 kullanır ve ana \(x\)-aralığını iş parçacıklarına bölerek paralelleştirir. Ayrıca şu kontrol değerleri bulunur: \(S(5)=198\), \(S(10)=396\), \(S(25)=5766\) ve brute-force doğrulaması olarak \(S(45)=34518\).

Karmaşıklık Analizi

\(r=5^{10}\) ve \(\texttt{limit}=2r\) olsun. SPF süzgecinin kurulması \(O(\texttt{limit}\log\log \texttt{limit})\) zaman ve \(O(\texttt{limit})\) bellek gerektirir. SPF tablosunu kullanarak faktör tablolarını doldurmak ortalama olarak yine yaklaşık \(O(\texttt{limit}\log\log \texttt{limit})\) maliyettedir. Son tarama \(x=1,\dots,r-1\) üzerinde \(O(r)\) zamanda çalışır ve her adımda yalnızca sabit sayıda tablo erişimi yapar. Toplam bellek tüketimi \(O(r)\)'dir; baskın kısım önceden hesaplanan dizilerdir. C++ sürümündeki paralelleştirme yalnızca pratik çalışma süresini düşürür, asimptotik karmaşıklığı değiştirmez.

Dipnotlar ve Kaynakça

  1. Problem sayfası: https://projecteuler.net/problem=360
  2. İki kare toplamı teoremi: Wikipedia — Sum of two squares theorem
  3. \(r_2(n)\) fonksiyonu: Wikipedia — Sum of two squares function
  4. Eratosthenes eleği ve SPF tablosu: Wikipedia — Sieve of Eratosthenes

Resumen del Problema

Se define

$$S(r)=\sum_{x^2+y^2+z^2=r^2}\left(|x|+|y|+|z|\right),$$

donde la suma recorre todos los puntos enteros de la esfera de radio \(r\). En el problema 360 el radio buscado es \(r=10^{10}=2^{10}5^{10}\). Enumerar directamente todos los triples \((x,y,z)\) es inviable, así que la solución transforma el problema tridimensional en un barrido unidimensional apoyado en la función de representaciones como suma de dos cuadrados.

Enfoque Matemático

Las soluciones del repositorio usan la misma cadena de ideas:

$$\text{(i) la simetría reduce todo a }r_2(r^2-x^2), \qquad \text{(ii) }10^{10}=2^{10}5^{10}\text{ permite separar las potencias de }2\text{ y }5.$$

Paso 1: Reducir la suma en tres variables a una sola coordenada

Las coordenadas \(x\), \(y\) y \(z\) aparecen de manera totalmente simétrica. Por tanto, la contribución total de \(|x|\), \(|y|\) y \(|z|\) es la misma; basta calcular una de ellas y multiplicar por \(3\).

Si fijamos un \(x>0\), cada par \((y,z)\in\mathbb Z^2\) que satisface

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

produce dos puntos de la esfera, uno con \(x\) y otro con \(-x\). Así, la contribución debida a la coordenada \(x\) es

$$2x\,r_2(r^2-x^2),$$

donde \(r_2(n)\) es el número de pares enteros ordenados \((u,v)\) tales que \(u^2+v^2=n\). Multiplicando por \(3\), obtenemos

$$\boxed{S(r)=6\sum_{x=1}^{r} x\,r_2(r^2-x^2).}$$

El caso \(x=0\) no aporta nada porque está multiplicado por \(x\). El extremo \(x=r\) sí es especial: entonces \(r^2-x^2=0\) y \(r_2(0)=1\) por el único par \((0,0)\).

Paso 2: Aplicar la función de suma de dos cuadrados

Para cada \(x\), escribimos

$$r^2-x^2=(r-x)(r+x)=ab,$$

con \(a=r-x\) y \(b=r+x\). El problema queda reducido a evaluar \(r_2(ab)\).

Si

$$n=2^\alpha \prod_{p\equiv 1 \pmod{4}} p^{e_p}\prod_{q\equiv 3 \pmod{4}} q^{f_q},$$

la fórmula clásica para el número de representaciones ordenadas con signos es

$$r_2(n)= \begin{cases} 0, & \text{si algún } f_q \text{ es impar},\\ 4\displaystyle\prod_{p\equiv 1 \pmod{4}}(e_p+1), & \text{en caso contrario}. \end{cases}$$

En consecuencia, los primos \(3 \pmod{4}\) solo actúan como prueba de paridad: un exponente impar anula el término. Los primos \(1 \pmod{4}\) aportan factores \(e_p+1\). La potencia de \(2\) no cambia el valor de \(r_2(n)\), por eso el código no necesita almacenar un exponente especial para el \(2\).

Paso 3: Por qué bastan esas tablas cuando \(r=5^b\)

El código resuelve primero el caso \(r=5^b\). Para \(x\in\{1,\dots,r-1\}\), define

$$a=r-x,\qquad b=r+x.$$

Cualquier divisor común de \(a\) y \(b\) divide a la vez \(a+b=2r\) y \(b-a=2x\), así que

$$\gcd(a,b)\mid 2r=2\cdot 5^b.$$

Ese hecho concentra toda la simplificación. Significa que cualquier primo impar \(p\neq 5\) puede dividir como mucho a uno de los factores \(a\) o \(b\). De ahí salen tres consecuencias:

\(1.\) Un primo \(q\equiv 3 \pmod{4}\) tiene exponente impar en \(ab\) exactamente cuando lo tiene en \(a\) o en \(b\). Basta conocer si cada factor es “bad”.

\(2.\) Un primo \(p\equiv 1 \pmod{4}\), \(p\neq 5\), contribuye de manera independiente desde \(a\) y \(b\), porque no puede aparecer en ambos a la vez.

\(3.\) El primo \(5\) debe tratarse aparte, ya que es el único primo congruente con \(1 \pmod{4}\) que sí puede estar en los dos factores.

Por eso las implementaciones precalculan, para todo \(m\le 2r\):

$$\texttt{bad}(m)= \begin{cases} 1,& \text{si algún } q\equiv 3 \pmod{4} \text{ aparece en }m\text{ con exponente impar},\\ 0,& \text{en caso contrario}, \end{cases}$$

$$\nu_5(m)=v_5(m),$$

$$g(m)=\prod_{\substack{p\equiv 1 \pmod{4}\\ p\neq 5}}(e_p(m)+1).$$

Si \(a\) o \(b\) es bad, entonces \(r_2(ab)=0\). En caso contrario, el exponente de \(5\) en \(ab\) es \(\nu_5(a)+\nu_5(b)\), y los demás primos \(1 \pmod{4}\) se separan limpiamente, de modo que

$$\boxed{r_2(ab)=4\,g(a)g(b)\bigl(\nu_5(a)+\nu_5(b)+1\bigr).}$$

Esa es exactamente la expresión utilizada por los programas en C++, Python y Java.

Paso 4: Recuperar el radio \(2^a5^b\)

Ahora tomemos \(R=2^a5^b\). Como \(R^2\) es múltiplo de \(4\), toda solución de

$$x^2+y^2+z^2=R^2$$

debe tener las tres coordenadas pares. En efecto, módulo \(4\) cada cuadrado vale \(0\) o \(1\), y la única forma de sumar \(0\) con tres de esos residuos es \(0+0+0\). Repitiendo el argumento \(a\) veces, toda solución se escribe como

$$x=2^a x',\qquad y=2^a y',\qquad z=2^a z',$$

con

$$x'^2+y'^2+z'^2=5^{2b}.$$

Esto da una biyección entre los puntos enteros de la esfera de radio \(2^a5^b\) y los de radio \(5^b\). Además, la suma \(|x|+|y|+|z|\) se escala por el mismo factor \(2^a\), así que

$$\boxed{S(2^a5^b)=2^aS(5^b).}$$

Por eso el código calcula primero el caso \(5^b\) y después multiplica por \(2^a\) mediante un desplazamiento binario.

Paso 5: Ejemplo trabajado con \(r=5\)

El pequeño caso \(r=5\) permite comprobar toda la derivación a mano:

$$S(5)=6\sum_{x=1}^{5}x\,r_2(25-x^2).$$

Veamos los términos:

\(x=1\): \(25-1=24=2^3\cdot 3\), así que un primo \(3 \pmod{4}\) aparece con exponente impar y \(r_2(24)=0\).

\(x=2\): \(25-4=21=3\cdot 7\), de nuevo imposible; \(r_2(21)=0\).

\(x=3\): \(25-9=16\), luego \(r_2(16)=4\) por \((\pm4,0)\) y \((0,\pm4)\).

\(x=4\): \(25-16=9\), luego \(r_2(9)=4\) por \((\pm3,0)\) y \((0,\pm3)\).

\(x=5\): \(25-25=0\), así que \(r_2(0)=1\).

En consecuencia,

$$S(5)=6\cdot 3\cdot 4 + 6\cdot 4\cdot 4 + 6\cdot 5 = 72+96+30=198.$$

La ley de escala da inmediatamente

$$S(10)=2S(5)=396,$$

valor que coincide con el checkpoint del código.

Cómo funciona el código

La solución en C++ construye primero una tabla SPF de factores primos mínimos hasta \(2r\), con \(r=5^{10}\), mediante build_spf. Luego build_factor_tables factoriza cada entero \(1\le m\le 2r\) y guarda las tres piezas necesarias: la marca bad, el exponente de \(5\) y el producto multiplicativo de los primos \(1 \pmod{4}\) distintos de \(5\).

El sumatorio principal recorre \(x=1,\dots,r-1\), calcula \(a=r-x\) y \(b=r+x\), descarta de inmediato los casos bad y, en los demás, evalúa

$$\texttt{ways}=4\,g(a)g(b)\bigl(\nu_5(a)+\nu_5(b)+1\bigr).$$

La cantidad añadida al total es

$$6x\cdot \texttt{ways}.$$

Al terminar, añade el término de borde \(6r\) correspondiente a \(x=r\), y después aplica el factor \(2^{10}\) mediante un desplazamiento a la izquierda.

Los archivos de Python y Java implementan la misma fórmula. La versión Java usa BigInteger para la suma final. La versión C++ usa unsigned __int128 y paraleliza el rango de \(x\) entre varios hilos. Además incluye comprobaciones explícitas: \(S(5)=198\), \(S(10)=396\), \(S(25)=5766\) y una verificación por fuerza bruta \(S(45)=34518\).

Complejidad

Sea \(r=5^{10}\) y \(\texttt{limit}=2r\). Construir la criba SPF cuesta \(O(\texttt{limit}\log\log \texttt{limit})\) tiempo y \(O(\texttt{limit})\) memoria. Llenar las tablas de factores eliminando repetidamente factores SPF es casi lineal en promedio, del mismo orden \(O(\texttt{limit}\log\log \texttt{limit})\). El barrido final sobre \(x=1,\dots,r-1\) es \(O(r)\) con solo accesos O(1) a tablas en cada iteración. La memoria total es \(O(r)\), dominada por los arreglos precalculados. La paralelización en C++ mejora el tiempo real de ejecución, pero no cambia la complejidad asintótica.

Notas y Referencias

  1. Página del problema: https://projecteuler.net/problem=360
  2. Teorema de la suma de dos cuadrados: Wikipedia — Teorema de Fermat sobre la suma de dos cuadrados
  3. Función \(r_2(n)\): Wikipedia — Sum of two squares function
  4. Criba de Eratóstenes y SPF: Wikipedia — Criba de Eratóstenes

问题概述

定义

$$S(r)=\sum_{x^2+y^2+z^2=r^2}\left(|x|+|y|+|z|\right),$$

求和对象是半径为 \(r\) 的球面上所有整数点。Project Euler 360 要求的半径是 \(r=10^{10}=2^{10}5^{10}\)。如果直接枚举所有三元组 \((x,y,z)\),规模完全不可接受,所以实现把三维问题改写成一维扫描,并把核心工作交给“两平方和表示数”函数。

数学方法

仓库里的几份解法都围绕同一条思路展开:

$$\text{(i) 用对称性把问题压到 }r_2(r^2-x^2), \qquad \text{(ii) 用 }10^{10}=2^{10}5^{10}\text{ 把 }2\text{ 的部分和 }5\text{ 的部分分开。}$$

步骤 1:把三变量求和化成单坐标求和

\(|x|\)、\(|y|\)、\(|z|\) 三项在球面对称下地位完全相同,所以三部分总贡献相等。我们只需先计算 \(|x|\) 的贡献,再乘以 \(3\)。

固定一个正的 \(x\) 后,所有满足

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

的整数对 \((y,z)\in\mathbb Z^2\),在球面上都会对应两个点,一个取坐标 \(x\),另一个取坐标 \(-x\)。因此 \(x\)-坐标对总和的贡献是

$$2x\,r_2(r^2-x^2),$$

其中 \(r_2(n)\) 表示有多少个有序整数对 \((u,v)\) 满足 \(u^2+v^2=n\)。再乘上三个坐标的对称因子,就得到

$$\boxed{S(r)=6\sum_{x=1}^{r} x\,r_2(r^2-x^2).}$$

\(x=0\) 的项本来就乘着 \(x\),所以不贡献任何值。端点 \(x=r\) 需要单独注意,因为这时 \(r^2-x^2=0\),而 \(r_2(0)=1\),对应唯一的 \((0,0)\)。

步骤 2:使用两平方和函数

对每个 \(x\),写成

$$r^2-x^2=(r-x)(r+x)=ab,$$

这里 \(a=r-x\),\(b=r+x\)。这样问题就变成了如何快速求 \(r_2(ab)\)。

$$n=2^\alpha \prod_{p\equiv 1 \pmod{4}} p^{e_p}\prod_{q\equiv 3 \pmod{4}} q^{f_q},$$

则经典结论告诉我们,带符号、带顺序的表示个数满足

$$r_2(n)= \begin{cases} 0, & \text{如果某个 } f_q \text{ 为奇数},\\ 4\displaystyle\prod_{p\equiv 1 \pmod{4}}(e_p+1), & \text{否则}. \end{cases}$$

这说明 \(3 \pmod{4}\) 的素数只起“筛掉不可能情况”的作用,只要有奇次幂就没有表示;而 \(1 \pmod{4}\) 的素数则提供乘法因子 \(e_p+1\)。至于 \(2\) 的幂次,并不会改变 \(r_2(n)\) 的值,所以实现里不需要专门为 \(2\) 再建一个表。

步骤 3:为什么在 \(r=5^b\) 时只需要这些预处理表

代码先解决 \(r=5^b\) 的情形。对 \(x\in\{1,\dots,r-1\}\),设

$$a=r-x,\qquad b=r+x.$$

\(a\) 和 \(b\) 的任意公因子同时整除 \(a+b=2r\) 与 \(b-a=2x\),因此

$$\gcd(a,b)\mid 2r=2\cdot 5^b.$$

这一步非常关键。它说明任何奇素数 \(p\neq 5\) 最多只能出现在 \(a\) 或 \(b\) 其中一个因子里,于是:

\(1.\) 对于 \(q\equiv 3 \pmod{4}\) 的素数,它在 \(ab\) 中的指数为奇数,当且仅当它在 \(a\) 或 \(b\) 中的指数为奇数。因此只要知道每个因子是否“bad”就够了。

\(2.\) 对于 \(p\equiv 1 \pmod{4}\)、且 \(p\neq 5\) 的素数,它们对 \(a\) 和 \(b\) 的贡献可以完全分开相乘,因为同一个这样的素数不可能同时整除两边。

\(3.\) 素数 \(5\) 需要单独处理,因为它是唯一一个既满足 \(1 \pmod{4}\) 又可能同时出现在 \(a\) 和 \(b\) 中的素数。

因此,程序会对所有 \(m\le 2r\) 预处理三张表:

$$\texttt{bad}(m)= \begin{cases} 1,& \text{若 }m\text{ 中存在 } q\equiv 3 \pmod{4}\text{ 的奇次幂},\\ 0,& \text{否则}, \end{cases}$$

$$\nu_5(m)=v_5(m),$$

$$g(m)=\prod_{\substack{p\equiv 1 \pmod{4}\\ p\neq 5}}(e_p(m)+1).$$

如果 \(a\) 或 \(b\) 是 bad,那么 \(r_2(ab)=0\)。否则 \(5\) 在 \(ab\) 中的指数就是 \(\nu_5(a)+\nu_5(b)\),其余 \(1 \pmod{4}\) 的素数贡献可以直接拆开,于是

$$\boxed{r_2(ab)=4\,g(a)g(b)\bigl(\nu_5(a)+\nu_5(b)+1\bigr).}$$

这正是 C++、Python 和 Java 解法里使用的公式。

步骤 4:恢复到半径 \(2^a5^b\)

现在考虑一般半径 \(R=2^a5^b\)。因为 \(R^2\) 被 \(4\) 整除,任何满足

$$x^2+y^2+z^2=R^2$$

的解都必须让 \(x,y,z\) 全部是偶数。原因是模 \(4\) 下每个平方只能是 \(0\) 或 \(1\),三个这样的数之和要得到 \(0\),唯一可能就是 \(0+0+0\)。把这个论证重复 \(a\) 次,就得到所有解都可以唯一写成

$$x=2^a x',\qquad y=2^a y',\qquad z=2^a z',$$

并且满足

$$x'^2+y'^2+z'^2=5^{2b}.$$

这给出了半径 \(2^a5^b\) 与半径 \(5^b\) 两个球面整数点集合之间的双射。与此同时,\(|x|+|y|+|z|\) 也整体放大了 \(2^a\) 倍,因此

$$\boxed{S(2^a5^b)=2^aS(5^b).}$$

代码先算 \(5^b\) 的情形、最后再左移 \(a\) 位,正是这个结论的直接实现。

步骤 5:以 \(r=5\) 为例做一次完整验算

小半径 \(r=5\) 可以把整套推导手算一遍:

$$S(5)=6\sum_{x=1}^{5}x\,r_2(25-x^2).$$

逐项来看:

\(x=1\):\(25-1=24=2^3\cdot 3\),存在 \(3 \pmod{4}\) 的素数奇次幂,所以 \(r_2(24)=0\)。

\(x=2\):\(25-4=21=3\cdot 7\),同样不可能表示,\(r_2(21)=0\)。

\(x=3\):\(25-9=16\),故 \(r_2(16)=4\),对应 \((\pm4,0)\)、\((0,\pm4)\)。

\(x=4\):\(25-16=9\),故 \(r_2(9)=4\),对应 \((\pm3,0)\)、\((0,\pm3)\)。

\(x=5\):\(25-25=0\),故 \(r_2(0)=1\)。

于是

$$S(5)=6\cdot 3\cdot 4 + 6\cdot 4\cdot 4 + 6\cdot 5 = 72+96+30=198.$$

再利用缩放关系,立刻得到

$$S(10)=2S(5)=396,$$

这与实现中的 checkpoint 完全一致。

代码如何对应这个公式

C++ 程序先用 build_spf 构建到 \(2r\) 为止的最小素因子表,其中 \(r=5^{10}\)。随后 build_factor_tables 把每个 \(1\le m\le 2r\) 分解质因数,并保存上面三种信息:bad 标记、\(5\) 的指数、以及除了 \(5\) 之外所有 \(1 \pmod{4}\) 素数贡献出的乘法因子。

主循环遍历 \(x=1,\dots,r-1\),算出 \(a=r-x\)、\(b=r+x\),如果有 bad 就直接跳过;否则计算

$$\texttt{ways}=4\,g(a)g(b)\bigl(\nu_5(a)+\nu_5(b)+1\bigr).$$

然后把

$$6x\cdot \texttt{ways}$$

加入总和。循环结束后,再补上端点 \(x=r\) 的贡献 \(6r\),最后通过左移把结果乘上 \(2^{10}\)。

Python 和 Java 文件实现的是同一套算术。Java 用 BigInteger 保存最终大整数。C++ 用 unsigned __int128 保存结果,并把主循环的 \(x\) 区间切成多个线程并行处理。C++ 里还写了显式校验:\(S(5)=198\)、\(S(10)=396\)、\(S(25)=5766\),以及通过暴力枚举验证的 \(S(45)=34518\)。

复杂度分析

设 \(r=5^{10}\),\(\texttt{limit}=2r\)。构建 SPF 筛需要 \(O(\texttt{limit}\log\log \texttt{limit})\) 时间和 \(O(\texttt{limit})\) 空间。利用 SPF 反复剥离质因子来填充三张表,平均下来仍然接近线性,约为 \(O(\texttt{limit}\log\log \texttt{limit})\)。最后对 \(x=1,\dots,r-1\) 的扫描是 \(O(r)\),每一步只做常数次表访问。总空间复杂度是 \(O(r)\),主要由预处理数组占据。C++ 的多线程只改善实际运行时间,不改变渐近复杂度。

参考资料

  1. 题目页面: https://projecteuler.net/problem=360
  2. 两平方和定理: Wikipedia — 费马两平方和定理
  3. \(r_2(n)\) 函数: Wikipedia — Sum of two squares function
  4. 埃拉托斯特尼筛法与 SPF: Wikipedia — 埃拉托斯特尼筛法

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

Рассматривается функция

$$S(r)=\sum_{x^2+y^2+z^2=r^2}\left(|x|+|y|+|z|\right),$$

где суммирование ведется по всем целочисленным точкам на сфере радиуса \(r\). В задаче 360 нужно вычислить это значение для \(r=10^{10}=2^{10}5^{10}\). Полный перебор всех троек \((x,y,z)\) слишком дорог, поэтому решение сводит трехмерную сумму к одномерному проходу и использует классическую функцию числа представлений в виде суммы двух квадратов.

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

Во всех файлах решения используется одна и та же схема:

$$\text{(i) симметрия сводит задачу к }r_2(r^2-x^2), \qquad \text{(ii) }10^{10}=2^{10}5^{10}\text{ позволяет отделить степени }2\text{ от степеней }5.$$

Шаг 1: Свести сумму по трем координатам к одной координате

Координаты \(x\), \(y\) и \(z\) входят симметрично, значит суммарные вклады \(|x|\), \(|y|\) и \(|z|\) одинаковы. Поэтому достаточно посчитать вклад \(|x|\) и затем умножить его на \(3\).

Если фиксировать положительное \(x\), то каждый целочисленный пар \((y,z)\in\mathbb Z^2\), удовлетворяющий

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

дает две точки на сфере: с координатой \(x\) и с координатой \(-x\). Следовательно, вклад \(x\)-координаты равен

$$2x\,r_2(r^2-x^2),$$

где \(r_2(n)\) обозначает число упорядоченных целочисленных пар \((u,v)\), для которых \(u^2+v^2=n\). Умножая на \(3\), получаем

$$\boxed{S(r)=6\sum_{x=1}^{r} x\,r_2(r^2-x^2).}$$

Случай \(x=0\) не влияет на результат, потому что стоит множитель \(x\). А вот точка \(x=r\) особая: тогда \(r^2-x^2=0\), и \(r_2(0)=1\) благодаря единственной паре \((0,0)\).

Шаг 2: Использовать функцию суммы двух квадратов

Для каждого \(x\) запишем

$$r^2-x^2=(r-x)(r+x)=ab,$$

где \(a=r-x\), \(b=r+x\). Теперь нужно быстро находить \(r_2(ab)\).

Если

$$n=2^\alpha \prod_{p\equiv 1 \pmod{4}} p^{e_p}\prod_{q\equiv 3 \pmod{4}} q^{f_q},$$

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

$$r_2(n)= \begin{cases} 0, & \text{если некоторый } f_q \text{ нечетен},\\ 4\displaystyle\prod_{p\equiv 1 \pmod{4}}(e_p+1), & \text{иначе}. \end{cases}$$

Итак, простые \(3 \pmod{4}\) играют роль фильтра по четности: нечетная степень сразу делает представление невозможным. Простые \(1 \pmod{4}\) дают множители \(e_p+1\). Степень двойки на значение \(r_2(n)\) не влияет, поэтому отдельная таблица для простого \(2\) не нужна.

Шаг 3: Почему для случая \(r=5^b\) достаточно именно таких таблиц

Сначала код решает задачу для радиусов вида \(r=5^b\). Для \(x\in\{1,\dots,r-1\}\) введем

$$a=r-x,\qquad b=r+x.$$

Любой общий делитель чисел \(a\) и \(b\) делит одновременно \(a+b=2r\) и \(b-a=2x\), а значит

$$\gcd(a,b)\mid 2r=2\cdot 5^b.$$

Именно здесь появляется главное упрощение. Из этого следует, что любой нечетный простой \(p\neq 5\) может делить не более одного из чисел \(a\) и \(b\). Поэтому:

\(1.\) Для простого \(q\equiv 3 \pmod{4}\) степень в \(ab\) нечетна тогда и только тогда, когда она нечетна в \(a\) или в \(b\). Достаточно хранить для каждого множителя только флаг bad.

\(2.\) Для простого \(p\equiv 1 \pmod{4}\), \(p\neq 5\), вклад распадается независимо по \(a\) и \(b\), потому что такой простой не может входить в оба множителя одновременно.

\(3.\) Простое число \(5\) нужно обрабатывать отдельно, так как это единственный простой \(1 \pmod{4}\), который может появляться и в \(a\), и в \(b\).

Поэтому решения заранее вычисляют для каждого \(m\le 2r\):

$$\texttt{bad}(m)= \begin{cases} 1,& \text{если в }m\text{ есть простой } q\equiv 3 \pmod{4}\text{ в нечетной степени},\\ 0,& \text{иначе}, \end{cases}$$

$$\nu_5(m)=v_5(m),$$

$$g(m)=\prod_{\substack{p\equiv 1 \pmod{4}\\ p\neq 5}}(e_p(m)+1).$$

Если хотя бы один из множителей bad, то \(r_2(ab)=0\). Иначе степень пятерки в произведении \(ab\) равна \(\nu_5(a)+\nu_5(b)\), а остальные простые \(1 \pmod{4}\) раскладываются независимо, поэтому

$$\boxed{r_2(ab)=4\,g(a)g(b)\bigl(\nu_5(a)+\nu_5(b)+1\bigr).}$$

Именно это выражение используется в реализациях на C++, Python и Java.

Шаг 4: Вернуть общий радиус \(2^a5^b\)

Пусть теперь \(R=2^a5^b\). Так как \(R^2\) делится на \(4\), любая точка, удовлетворяющая

$$x^2+y^2+z^2=R^2,$$

обязана иметь все три координаты четными. Причина проста: по модулю \(4\) каждый квадрат равен либо \(0\), либо \(1\), и единственный способ получить сумму \(0\) из трех таких остатков это \(0+0+0\). Повторяя это рассуждение \(a\) раз, получаем, что всякое решение имеет вид

$$x=2^a x',\qquad y=2^a y',\qquad z=2^a z',$$

где

$$x'^2+y'^2+z'^2=5^{2b}.$$

Это дает биекцию между целочисленными точками на сферах радиусов \(2^a5^b\) и \(5^b\). При этом сумма модулей координат умножается на тот же коэффициент \(2^a\), следовательно

$$\boxed{S(2^a5^b)=2^aS(5^b).}$$

Этим и объясняется завершающий сдвиг влево в коде после вычисления случая \(5^b\).

Шаг 5: Разобранный пример для \(r=5\)

Малый контрольный случай \(r=5\) можно полностью проверить вручную:

$$S(5)=6\sum_{x=1}^{5}x\,r_2(25-x^2).$$

Разберем слагаемые:

\(x=1\): \(25-1=24=2^3\cdot 3\), есть простой \(3 \pmod{4}\) в нечетной степени, значит \(r_2(24)=0\).

\(x=2\): \(25-4=21=3\cdot 7\), снова представления нет, поэтому \(r_2(21)=0\).

\(x=3\): \(25-9=16\), значит \(r_2(16)=4\), соответствуют пары \((\pm4,0)\) и \((0,\pm4)\).

\(x=4\): \(25-16=9\), значит \(r_2(9)=4\), соответствуют пары \((\pm3,0)\) и \((0,\pm3)\).

\(x=5\): \(25-25=0\), значит \(r_2(0)=1\).

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

$$S(5)=6\cdot 3\cdot 4 + 6\cdot 4\cdot 4 + 6\cdot 5 = 72+96+30=198.$$

Далее по закону масштабирования сразу получаем

$$S(10)=2S(5)=396,$$

и это совпадает с контрольным значением в реализации.

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

Сначала C++-версия строит таблицу SPF до \(2r\), где \(r=5^{10}\), с помощью build_spf. Затем функция build_factor_tables факторизует все числа \(1\le m\le 2r\) и сохраняет три величины: флаг bad, показатель степени пятерки и мультипликативный множитель от всех простых \(1 \pmod{4}\), кроме \(5\).

Основной цикл проходит по \(x=1,\dots,r-1\), вычисляет \(a=r-x\) и \(b=r+x\), сразу пропускает bad-случаи и иначе считает

$$\texttt{ways}=4\,g(a)g(b)\bigl(\nu_5(a)+\nu_5(b)+1\bigr).$$

После этого в сумму добавляется

$$6x\cdot \texttt{ways}.$$

По окончании цикла добавляется крайний вклад \(6r\) для \(x=r\), а затем результат умножается на \(2^{10}\) сдвигом влево.

Файлы Python и Java реализуют ту же арифметику. В Java итог хранится в BigInteger. В C++ используется unsigned __int128, а диапазон по \(x\) распараллеливается между потоками. Кроме того, в C++ есть явные проверки: \(S(5)=198\), \(S(10)=396\), \(S(25)=5766\), а также brute-force проверка \(S(45)=34518\).

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

Пусть \(r=5^{10}\), а \(\texttt{limit}=2r\). Построение SPF-решета требует \(O(\texttt{limit}\log\log \texttt{limit})\) времени и \(O(\texttt{limit})\) памяти. Заполнение факторных таблиц с помощью последовательного удаления SPF-факторов в среднем также близко к линейному и имеет порядок \(O(\texttt{limit}\log\log \texttt{limit})\). Финальный проход по \(x=1,\dots,r-1\) занимает \(O(r)\) времени при константном числе обращений к таблицам на шаг. Общий расход памяти равен \(O(r)\) и определяется в основном предварительно вычисленными массивами. Параллелизация в C++ уменьшает только реальное время выполнения, но не меняет асимптотику.

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

  1. Страница задачи: https://projecteuler.net/problem=360
  2. Теорема о сумме двух квадратов: Wikipedia — Теорема Ферма о сумме двух квадратов
  3. Функция \(r_2(n)\): Wikipedia — Sum of two squares function
  4. Решето Эратосфена и SPF: Wikipedia — Решето Эратосфена

ملخص المسألة

نعرف

$$S(r)=\sum_{x^2+y^2+z^2=r^2}\left(|x|+|y|+|z|\right),$$

حيث يؤخذ المجموع على جميع النقاط الصحيحة الواقعة على سطح الكرة ذات نصف القطر \(r\). في المسألة 360 نريد القيمة عند \(r=10^{10}=2^{10}5^{10}\). التعداد المباشر لكل الثلاثيات \((x,y,z)\) غير عملي تمامًا، ولذلك يحوّل الحل المسألة من فضاء ثلاثي الأبعاد إلى مسح أحادي البعد يعتمد على دالة عدد التمثيلات كمجموع مربعين.

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

تستعمل ملفات الحل في المستودع الفكرة نفسها:

$$\text{(i) التماثل يختزل المسألة إلى }r_2(r^2-x^2), \qquad \text{(ii) التحليل }10^{10}=2^{10}5^{10}\text{ يسمح بفصل قوى }2\text{ عن قوى }5.$$

الخطوة 1: اختزال المجموع ثلاثي المتغيرات إلى إحداثي واحد

المقادير \(|x|\) و\(|y|\) و\(|z|\) متناظرة تمامًا، ولذلك فإن مساهمة كل منها في المجموع الكلي متساوية. لذا يكفي أن نحسب مساهمة \(|x|\) فقط ثم نضرب الناتج في \(3\).

إذا ثبتنا \(x>0\)، فإن كل زوج صحيح \((y,z)\in\mathbb Z^2\) يحقق

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

يعطي نقطتين على الكرة: واحدة عند الإحداثي \(x\) وأخرى عند \(-x\). لذلك تكون مساهمة إحداثي \(x\) هي

$$2x\,r_2(r^2-x^2),$$

حيث \(r_2(n)\) هو عدد الأزواج الصحيحة المرتبة \((u,v)\) التي تحقق \(u^2+v^2=n\). وبعد الضرب في \(3\) نحصل على

$$\boxed{S(r)=6\sum_{x=1}^{r} x\,r_2(r^2-x^2).}$$

الحد الموافق لـ \(x=0\) لا يضيف شيئًا لأنه مضروب في \(x\). أما الحالة \(x=r\) فهي خاصة، لأن \(r^2-x^2=0\)، ومن ثم \(r_2(0)=1\) بسبب الزوج الوحيد \((0,0)\).

الخطوة 2: استخدام دالة مجموع مربعين

لكل قيمة \(x\) نكتب

$$r^2-x^2=(r-x)(r+x)=ab,$$

حيث \(a=r-x\) و\(b=r+x\). وهكذا تصبح المهمة هي حساب \(r_2(ab)\).

إذا كان

$$n=2^\alpha \prod_{p\equiv 1 \pmod{4}} p^{e_p}\prod_{q\equiv 3 \pmod{4}} q^{f_q},$$

فإن الصيغة الكلاسيكية لعدد التمثيلات المرتبة مع الإشارات هي

$$r_2(n)= \begin{cases} 0, & \text{إذا وُجد } f_q \text{ فردي},\\ 4\displaystyle\prod_{p\equiv 1 \pmod{4}}(e_p+1), & \text{في غير ذلك}. \end{cases}$$

إذن فالأوليات الموافقة لـ \(3 \pmod{4}\) تعمل فقط كاختبار زوجية: وجود أس فردي يعني عدم وجود تمثيل. أما الأوليات الموافقة لـ \(1 \pmod{4}\) فتعطي عوامل من الشكل \(e_p+1\). وأما أس \(2\) فلا يغيّر قيمة \(r_2(n)\)، ولهذا لا يحتاج الكود إلى جدول مستقل لقوى \(2\).

الخطوة 3: لماذا تكفي هذه الجداول عندما يكون \(r=5^b\)

يعالج الكود أولًا أنصاف الأقطار من الشكل \(r=5^b\). ولـ \(x\in\{1,\dots,r-1\}\) نعرّف

$$a=r-x,\qquad b=r+x.$$

أي قاسم مشترك بين \(a\) و\(b\) يقسم في الوقت نفسه \(a+b=2r\) و\(b-a=2x\)، ولذلك

$$\gcd(a,b)\mid 2r=2\cdot 5^b.$$

هذه هي نقطة التبسيط الأساسية. فهي تعني أن أي أولي فردي \(p\neq 5\) لا يمكن أن يقسم العاملين \(a\) و\(b\) معًا. ومن هنا نحصل على النتائج التالية:

\(1.\) إذا كان \(q\equiv 3 \pmod{4}\)، فإن أسه في \(ab\) يكون فرديًا إذا وفقط إذا كان أسه فرديًا في \(a\) أو في \(b\). لذا يكفي أن نعرف هل كل عامل bad أم لا.

\(2.\) إذا كان \(p\equiv 1 \pmod{4}\) و\(p\neq 5\)، فإن مساهمته تتفكك بشكل مستقل بين \(a\) و\(b\)، لأنه لا يستطيع أن يظهر في العاملين معًا.

\(3.\) الأولي \(5\) يحتاج معالجة خاصة، لأنه الأولي الوحيد الموافق لـ \(1 \pmod{4}\) الذي يمكن أن يظهر في \(a\) و\(b\) معًا.

لهذا السبب تسبق الحلول حساب ثلاث جداول لكل \(m\le 2r\):

$$\texttt{bad}(m)= \begin{cases} 1,& \text{إذا احتوى }m\text{ على أولي } q\equiv 3 \pmod{4}\text{ بأس فردي},\\ 0,& \text{خلاف ذلك}, \end{cases}$$

$$\nu_5(m)=v_5(m),$$

$$g(m)=\prod_{\substack{p\equiv 1 \pmod{4}\\ p\neq 5}}(e_p(m)+1).$$

إذا كان أحد العاملين \(a\) أو \(b\) من النوع bad فحينها \(r_2(ab)=0\). أما إذا كانا صالحين، فإن أس \(5\) في \(ab\) يساوي \(\nu_5(a)+\nu_5(b)\)، وباقي أوليات \(1 \pmod{4}\) تتفكك نظيفًا، ومن ثم

$$\boxed{r_2(ab)=4\,g(a)g(b)\bigl(\nu_5(a)+\nu_5(b)+1\bigr).}$$

وهذه هي الصيغة نفسها المستخدمة في حلول C++ وPython وJava.

الخطوة 4: استرجاع نصف القطر \(2^a5^b\)

لنأخذ الآن نصف القطر العام \(R=2^a5^b\). بما أن \(R^2\) يقبل القسمة على \(4\)، فإن أي حل للمعادلة

$$x^2+y^2+z^2=R^2$$

لا بد أن تكون فيه الإحداثيات الثلاثة جميعًا زوجية. والسبب أن كل مربع يساوي \(0\) أو \(1\) بترديد \(4\)، والطريقة الوحيدة للحصول على مجموع يساوي \(0\) من ثلاثة مثل هذه البواقي هي \(0+0+0\). وإذا كررنا هذا الاستدلال \(a\) مرة نحصل على أن كل حل يكتب على الصورة

$$x=2^a x',\qquad y=2^a y',\qquad z=2^a z',$$

بحيث

$$x'^2+y'^2+z'^2=5^{2b}.$$

وهذا يعطي تقابلًا واحدًا لواحد بين نقاط الشبكة على الكرة ذات نصف القطر \(2^a5^b\) ونقاط الشبكة على الكرة ذات نصف القطر \(5^b\). كما أن المقدار \(|x|+|y|+|z|\) يتمدد بالعامل نفسه \(2^a\)، وبالتالي

$$\boxed{S(2^a5^b)=2^aS(5^b).}$$

ولهذا السبب يحسب التنفيذ أولًا حالة \(5^b\)، ثم يضرب في \(2^a\) بإزاحة بتية في النهاية.

الخطوة 5: مثال محلول عند \(r=5\)

الحالة الصغيرة \(r=5\) تسمح بالتحقق اليدوي الكامل من الاشتقاق:

$$S(5)=6\sum_{x=1}^{5}x\,r_2(25-x^2).$$

نحسب الحدود واحدًا واحدًا:

\(x=1\): لدينا \(25-1=24=2^3\cdot 3\)، وفيه أولي من نوع \(3 \pmod{4}\) بأس فردي، ولذلك \(r_2(24)=0\).

\(x=2\): لدينا \(25-4=21=3\cdot 7\)، ومرة أخرى لا يوجد تمثيل، أي \(r_2(21)=0\).

\(x=3\): لدينا \(25-9=16\)، ومن ثم \(r_2(16)=4\) بسبب الأزواج \((\pm4,0)\) و\((0,\pm4)\).

\(x=4\): لدينا \(25-16=9\)، ومن ثم \(r_2(9)=4\) بسبب الأزواج \((\pm3,0)\) و\((0,\pm3)\).

\(x=5\): لدينا \(25-25=0\)، وبالتالي \(r_2(0)=1\).

إذًا

$$S(5)=6\cdot 3\cdot 4 + 6\cdot 4\cdot 4 + 6\cdot 5 = 72+96+30=198.$$

ثم يعطي قانون التحجيم مباشرة

$$S(10)=2S(5)=396,$$

وهو نفس مقدار التحقق الموجود في الكود.

كيف يظهر ذلك في الكود

يبني حل C++ أولًا جدول أصغر عامل أولي حتى \(2r\)، حيث \(r=5^{10}\)، باستعمال الدالة build_spf. بعد ذلك تقوم build_factor_tables بتحليل كل عدد \(1\le m\le 2r\) إلى عوامله الأولية وتخزين المعلومات الثلاث المطلوبة: علم bad، وأس \(5\)، والعامل الضربي الناتج عن جميع الأوليات الموافقة لـ \(1 \pmod{4}\) ما عدا \(5\).

الحلقة الرئيسية تمر على \(x=1,\dots,r-1\)، وتحسب \(a=r-x\) و\(b=r+x\)، ثم تتجاوز الحالة فورًا إذا كان أحدهما bad، وإلا تحسب

$$\texttt{ways}=4\,g(a)g(b)\bigl(\nu_5(a)+\nu_5(b)+1\bigr).$$

بعد ذلك تضاف الكمية

$$6x\cdot \texttt{ways}$$

إلى المجموع. وبعد انتهاء الحلقة تضاف مساهمة الطرف \(6r\) الموافقة لـ \(x=r\)، ثم يطبق عامل \(2^{10}\) بإزاحة إلى اليسار.

ملفا Python وJava يطبقان الصيغة الحسابية نفسها. في Java تُخزن النتيجة النهائية في BigInteger. أما C++ فيستخدم unsigned __int128، ويقسم مجال \(x\) بين عدة خيوط لتنفيذ الحلقة الرئيسية بالتوازي. كما يحتوي كود C++ على نقاط تحقق صريحة: \(S(5)=198\)، و\(S(10)=396\)، و\(S(25)=5766\)، إضافة إلى تحقق مباشر بالقوة الغاشمة أن \(S(45)=34518\).

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

لنضع \(r=5^{10}\) و\(\texttt{limit}=2r\). بناء غربال SPF يحتاج إلى زمن \(O(\texttt{limit}\log\log \texttt{limit})\) وذاكرة \(O(\texttt{limit})\). أما ملء جداول العوامل عبر نزع عوامل SPF بالتتابع فهو شبه خطي في المتوسط، وبالمرتبة نفسها تقريبًا \(O(\texttt{limit}\log\log \texttt{limit})\). وبعد ذلك تأتي الحلقة النهائية على \(x=1,\dots,r-1\)، وهي \(O(r)\) مع عدد ثابت من عمليات الوصول إلى الجداول في كل تكرار. الاستهلاك الكلي للذاكرة هو \(O(r)\)، وتسيطر عليه المصفوفات المسبقة الحساب. التوازي في نسخة C++ يحسن زمن التنفيذ العملي فقط، ولا يغير الرتبة التقاربية.

مراجع

  1. صفحة المسألة: https://projecteuler.net/problem=360
  2. نظرية مجموع مربعين: Wikipedia — نظرية فيرما حول مجموع مربعين
  3. دالة \(r_2(n)\): Wikipedia — Sum of two squares function
  4. غربال إراتوستينس وجداول SPF: Wikipedia — غربال إراتوستينس