Problem Summary

For each positive integer \(n\), define

$$A(n)=\#\left\{(x_1,\dots,x_8)\in(\mathbb{Z}/n\mathbb{Z})^8:\ x_1^2+x_2^2+x_3^2+x_4^2 \equiv x_5^2+x_6^2+x_7^2+x_8^2 \pmod n\right\}.$$

The task is to evaluate

$$S(N)=\sum_{n=1}^{N} A(n),\qquad S(12\,345\,678)\bmod 1\,001\,961\,001.$$

A direct search would require checking huge numbers of residue tuples for every modulus. The successful approach is to turn \(A(n)\) into a multiplicative arithmetic function and then compute its prime-power factors explicitly.

Mathematical Approach

The key observation is that the eight-variable congruence can be rewritten as a second moment of a four-variable residue count. After that, the Chinese remainder theorem and quadratic Gauss sums do the heavy lifting.

Step 1: Rewrite the congruence as a residue-distribution square

For a fixed modulus \(n\), let

$$r_n(t)=\#\left\{(a,b,c,d)\in(\mathbb{Z}/n\mathbb{Z})^4:\ a^2+b^2+c^2+d^2\equiv t\pmod n\right\}.$$

If two quadruples produce the same residue \(t\), then together they form one admissible octuple. Therefore

$$A(n)=\sum_{t\bmod n} r_n(t)^2.$$

So the original count is the second moment of the distribution of four-square residues modulo \(n\).

Step 2: Use the Chinese remainder theorem to prove multiplicativity

Suppose \(\gcd(m,n)=1\). A residue class modulo \(mn\) corresponds uniquely to a pair of residue classes modulo \(m\) and modulo \(n\). Under this identification, a quadruple modulo \(mn\) is just an independent pair of quadruples modulo \(m\) and modulo \(n\), so

$$r_{mn}(t_m,t_n)=r_m(t_m)\,r_n(t_n).$$

Squaring and summing over all residue pairs gives

$$A(mn)=A(m)\,A(n).$$

Hence it is enough to determine \(A(p^e)\) for prime powers, and then multiply the local factors together.

Step 3: Convert the residue count to quadratic Gauss sums

Introduce the additive characters

$$e_n(x)=\exp\!\left(\frac{2\pi i x}{n}\right),\qquad G_n(u)=\sum_{x\bmod n} e_n(ux^2).$$

By orthogonality of additive characters,

$$r_n(t)=\frac{1}{n}\sum_{u\bmod n} G_n(u)^4\,e_n(-ut).$$

Now take the second moment in \(t\). Parseval's identity for the finite Fourier transform yields

$$A(n)=\sum_{t\bmod n} r_n(t)^2=\frac{1}{n}\sum_{u\bmod n}\lvert G_n(u)\rvert^8.$$

This formula is the decisive reduction: instead of counting octuples directly, we only need the sizes of local quadratic Gauss sums.

Step 4: Evaluate the odd-prime-power factor

Let \(n=p^e\) with \(p\) odd. Write \(u=p^jv\) where \(0\le j<e\) and \(p\nmid v\). Then

$$G_{p^e}(u)=p^j G_{p^{e-j}}(v),\qquad \lvert G_{p^{e-j}}(v)\rvert=p^{(e-j)/2}.$$

So every frequency with valuation \(j\) contributes

$$\lvert G_{p^e}(u)\rvert^8=p^{4e+4j}.$$

There are \((p-1)p^{e-j-1}\) such frequencies, while \(u=0\) contributes \(p^{8e}\). Substituting into the character formula gives

$$A(p^e)=p^{7e}+(p-1)\sum_{j=0}^{e-1} p^{4e+3j-1}.$$

The same formula can be written recursively as

$$A(1)=1,\qquad A(p^e)=p^7A(p^{e-1})+(p-1)p^{4e-1}.$$

This is exactly the update rule used for odd prime powers in the implementations.

Step 5: Handle powers of \(2\) separately

The modulus \(2^e\) needs its own branch because quadratic Gauss sums over powers of \(2\) behave differently. For odd frequency modulo \(2^m\) with \(m\ge 2\), one has

$$\lvert G_{2^m}(v)\rvert=2^{(m+1)/2},$$

whereas the remaining modulus \(2\) case vanishes. Therefore only valuations \(j=0,1,\dots,e-2\) contribute when \(u=2^jv\) with \(v\) odd. The resulting local factor is

$$A(2^e)=2^{7e}+\sum_{j=0}^{e-2} 2^{4e+3+3j}.$$

In the iterative form used by the code, this becomes

$$A(2)=2^7,\qquad A(2^e)=2^7A(2^{e-1})+2^{4e+3}\quad (e\ge 2).$$

Step 6: Reconstruct the value for a general modulus

If

$$n=\prod_{i=1}^{k} p_i^{e_i},$$

then multiplicativity gives

$$A(n)=\prod_{i=1}^{k} A(p_i^{e_i}).$$

So the entire problem reduces to

$$\boxed{S(12\,345\,678)=\sum_{n=1}^{12\,345\,678}\prod_{p^e\parallel n}A(p^e)\pmod{1\,001\,961\,001}.}$$

Worked Example: \(n=10\) and the checkpoint up to \(10\)

Because \(10=2\cdot 5\), multiplicativity immediately gives

$$A(10)=A(2)\,A(5).$$

From the local formulas,

$$A(2)=2^7=128,$$

$$A(5)=5^7+4\cdot 5^3=78\,625.$$

Hence

$$A(10)=128\cdot 78\,625=10\,064\,000.$$

Applying the same formulas to all \(n\le 10\) gives

$$1,\ 128,\ 2\,241,\ 18\,432,\ 78\,625,\ 286\,848,\ 825\,601,\ 2\,392\,064,\ 4\,905\,441,\ 10\,064\,000,$$

and therefore

$$S(10)=18\,573\,381.$$

This is the small checkpoint verified by the implementation.

How the Code Works

The C++, Python, and Java implementations all follow the same structure. They first build a smallest-prime-factor table up to \(12\,345\,678\). That table makes it possible to factor every \(n\) in a fast incremental scan.

For each \(n\), the implementation strips equal prime factors to recover the prime-power decomposition. Every local factor \(A(p^e)\) is then evaluated modulo \(1\,001\,961\,001\): odd primes use the recurrence from Step 4, and the factor \(2^e\) uses the special recurrence from Step 5.

The local factors are multiplied to obtain \(A(n)\bmod 1\,001\,961\,001\), and a running sum accumulates

$$\sum_{n=1}^{12\,345\,678} A(n)\pmod{1\,001\,961\,001}.$$

The C++ implementation also performs tiny exact sanity checks on small moduli before the full computation, including \(A(4)=18\,432\) and \(S(10)=18\,573\,381\).

Complexity Analysis

Let \(N=12\,345\,678\). Building the smallest-prime-factor table with a linear sieve costs \(O(N)\) time and \(O(N)\) memory. Factoring all numbers up to \(N\) from that table requires total work proportional to the total number of prime factors encountered, which is \(O(N\log\log N)\) on average. The whole method is therefore near-linear in practice and uses \(O(N)\) memory.

Footnotes and References

  1. Problem page: Project Euler 875
  2. Chinese remainder theorem: Wikipedia — Chinese remainder theorem
  3. Gauss sum: Wikipedia — Gauss sum
  4. Parseval's theorem: Wikipedia — Parseval's theorem
  5. Four-square theorem: Wikipedia — Lagrange's four-square theorem

Problemzusammenfassung

Für jede positive ganze Zahl \(n\) definieren wir

$$A(n)=\#\left\{(x_1,\dots,x_8)\in(\mathbb{Z}/n\mathbb{Z})^8:\ x_1^2+x_2^2+x_3^2+x_4^2 \equiv x_5^2+x_6^2+x_7^2+x_8^2 \pmod n\right\}.$$

Gesucht ist

$$S(N)=\sum_{n=1}^{N} A(n),\qquad S(12\,345\,678)\bmod 1\,001\,961\,001.$$

Eine direkte Aufzählung wäre hoffnungslos, weil für jedes Modul enorme Mengen von Restklassen-Tupeln geprüft werden müssten. Die Lösung macht aus \(A(n)\) stattdessen eine multiplikative arithmetische Funktion und berechnet ihre Primzahlpotenz-Faktoren explizit.

Mathematischer Ansatz

Der entscheidende Trick besteht darin, die Kongruenz mit acht Variablen als zweites Moment einer Vier-Quadrate-Verteilung zu schreiben. Danach zerlegt der chinesische Restsatz das Problem, und quadratische Gaußsche Summen liefern die lokalen Faktoren.

Schritt 1: Schreibe die Kongruenz als Quadrat einer Restverteilung

Für festes \(n\) setzen wir

$$r_n(t)=\#\left\{(a,b,c,d)\in(\mathbb{Z}/n\mathbb{Z})^4:\ a^2+b^2+c^2+d^2\equiv t\pmod n\right\}.$$

Wählt man zwei Vierertupel mit demselben Rest \(t\), so entsteht genau ein gültiges Achtertupel. Also gilt

$$A(n)=\sum_{t\bmod n} r_n(t)^2.$$

Damit ist die ursprüngliche Aufgabe das zweite Moment der Verteilung aller Vier-Quadrate-Residuen modulo \(n\).

Schritt 2: Multiplikativität mit dem chinesischen Restsatz

Sei \(\gcd(m,n)=1\). Eine Restklasse modulo \(mn\) entspricht dann eindeutig einem Paar von Restklassen modulo \(m\) und modulo \(n\). Unter dieser Identifikation zerfällt ein Vierertupel modulo \(mn\) in ein unabhängiges Paar von Vierertupeln, daher

$$r_{mn}(t_m,t_n)=r_m(t_m)\,r_n(t_n).$$

Nach dem Quadrieren und Summieren über alle Restpaare folgt

$$A(mn)=A(m)\,A(n).$$

Es genügt also, \(A(p^e)\) für Primzahlpotenzen zu bestimmen und die lokalen Beiträge anschließend zu multiplizieren.

Schritt 3: Übergang zu quadratischen Gaußsummen

Wir führen die additiven Charaktere ein:

$$e_n(x)=\exp\!\left(\frac{2\pi i x}{n}\right),\qquad G_n(u)=\sum_{x\bmod n} e_n(ux^2).$$

Aus der Orthogonalität additiver Charaktere erhält man

$$r_n(t)=\frac{1}{n}\sum_{u\bmod n} G_n(u)^4\,e_n(-ut).$$

Nun bildet man das zweite Moment in \(t\). Parsevals Identität für die endliche Fourier-Transformation liefert

$$A(n)=\sum_{t\bmod n} r_n(t)^2=\frac{1}{n}\sum_{u\bmod n}\lvert G_n(u)\rvert^8.$$

Das ist die zentrale Reduktion: Statt Achtertupel direkt zu zählen, müssen nur noch die Größen lokaler quadratischer Gaußsummen bestimmt werden.

Schritt 4: Der lokale Faktor für ungerade Primzahlen

Sei \(n=p^e\) mit ungerader Primzahl \(p\). Schreibe \(u=p^jv\) mit \(0\le j<e\) und \(p\nmid v\). Dann gilt

$$G_{p^e}(u)=p^j G_{p^{e-j}}(v),\qquad \lvert G_{p^{e-j}}(v)\rvert=p^{(e-j)/2}.$$

Damit trägt jede Frequenz mit Bewertung \(j\) den Wert

$$\lvert G_{p^e}(u)\rvert^8=p^{4e+4j}.$$

Es gibt \((p-1)p^{e-j-1}\) solcher Frequenzen, während \(u=0\) den Beitrag \(p^{8e}\) liefert. Einsetzen ergibt

$$A(p^e)=p^{7e}+(p-1)\sum_{j=0}^{e-1} p^{4e+3j-1}.$$

Äquivalent dazu ist die Rekursion

$$A(1)=1,\qquad A(p^e)=p^7A(p^{e-1})+(p-1)p^{4e-1}.$$

Genau diese Vorschrift verwenden die Implementierungen für ungerade Primzahlpotenzen.

Schritt 5: Potenzen von \(2\) separat behandeln

Für \(2^e\) braucht man einen eigenen Zweig, weil sich quadratische Gaußsummen über Zweierpotenzen anders verhalten. Für ungerade Frequenz modulo \(2^m\) mit \(m\ge 2\) gilt

$$\lvert G_{2^m}(v)\rvert=2^{(m+1)/2},$$

während der verbleibende Modul-\(2\)-Fall verschwindet. Deshalb tragen nur die Bewertungen \(j=0,1,\dots,e-2\) bei, wenn \(u=2^jv\) mit ungeradem \(v\). Der lokale Faktor wird damit

$$A(2^e)=2^{7e}+\sum_{j=0}^{e-2} 2^{4e+3+3j}.$$

In der iterativen Form des Programms lautet das

$$A(2)=2^7,\qquad A(2^e)=2^7A(2^{e-1})+2^{4e+3}\quad (e\ge 2).$$

Schritt 6: Allgemeine Moduli aus lokalen Faktoren zusammensetzen

Ist

$$n=\prod_{i=1}^{k} p_i^{e_i},$$

so liefert die Multiplikativität

$$A(n)=\prod_{i=1}^{k} A(p_i^{e_i}).$$

Damit reduziert sich die gesamte Aufgabe auf

$$\boxed{S(12\,345\,678)=\sum_{n=1}^{12\,345\,678}\prod_{p^e\parallel n}A(p^e)\pmod{1\,001\,961\,001}.}$$

Durchgerechnetes Beispiel: \(n=10\) und der Kontrollwert bis \(10\)

Wegen \(10=2\cdot 5\) gilt sofort

$$A(10)=A(2)\,A(5).$$

Aus den lokalen Formeln folgt

$$A(2)=2^7=128,$$

$$A(5)=5^7+4\cdot 5^3=78\,625.$$

Daher

$$A(10)=128\cdot 78\,625=10\,064\,000.$$

Für alle \(n\le 10\) ergeben dieselben Formeln die Werte

$$1,\ 128,\ 2\,241,\ 18\,432,\ 78\,625,\ 286\,848,\ 825\,601,\ 2\,392\,064,\ 4\,905\,441,\ 10\,064\,000,$$

und damit

$$S(10)=18\,573\,381.$$

Genau dieser kleine Kontrollwert wird von der Implementierung überprüft.

Wie der Code arbeitet

Die C++-, Python- und Java-Implementierungen folgen derselben Struktur. Zuerst bauen sie bis \(12\,345\,678\) eine Tabelle der kleinsten Primfaktoren auf. Dadurch lässt sich jedes \(n\) in einem schnellen Durchlauf zerlegen.

Für jedes \(n\) werden gleiche Primfaktoren zusammengefasst, sodass die Primzahlpotenz-Zerlegung entsteht. Jeder lokale Faktor \(A(p^e)\) wird dann modulo \(1\,001\,961\,001\) berechnet: ungerade Primzahlen verwenden die Rekursion aus Schritt 4, und \(2^e\) benutzt die spezielle \(2\)-adische Rekursion aus Schritt 5.

Die lokalen Faktoren werden multipliziert, um \(A(n)\bmod 1\,001\,961\,001\) zu erhalten, und anschließend wird laufend

$$\sum_{n=1}^{12\,345\,678} A(n)\pmod{1\,001\,961\,001}$$

akkumuliert. Die C++-Implementierung enthält zusätzlich kleine exakte Selbsttests vor dem vollständigen Lauf, unter anderem \(A(4)=18\,432\) und \(S(10)=18\,573\,381\).

Komplexitätsanalyse

Sei \(N=12\,345\,678\). Die SPF-Tabelle wird mit einem linearen Sieb in \(O(N)\) Zeit und \(O(N)\) Speicher erzeugt. Das Faktorisieren aller Zahlen bis \(N\) benötigt insgesamt Arbeit proportional zur Gesamtzahl der auftretenden Primfaktoren, also im Mittel \(O(N\log\log N)\). Insgesamt ist das Verfahren praktisch nahezu linear und verwendet \(O(N)\) Speicher.

Fußnoten und Referenzen

  1. Problemseite: Project Euler 875
  2. Chinesischer Restsatz: Wikipedia — Chinese remainder theorem
  3. Gaußsche Summe: Wikipedia — Gauss sum
  4. Parsevals Identität: Wikipedia — Parseval's theorem
  5. Vier-Quadrate-Satz: Wikipedia — Lagrange's four-square theorem

Problem Özeti

Her pozitif \(n\) tamsayısı için

$$A(n)=\#\left\{(x_1,\dots,x_8)\in(\mathbb{Z}/n\mathbb{Z})^8:\ x_1^2+x_2^2+x_3^2+x_4^2 \equiv x_5^2+x_6^2+x_7^2+x_8^2 \pmod n\right\}$$

tanımını yapalım. İstenen değer

$$S(N)=\sum_{n=1}^{N} A(n),\qquad S(12\,345\,678)\bmod 1\,001\,961\,001.$$

Doğrudan tarama, her modül için çok büyük sayıda artık sınıfı sekizlisini kontrol etmeyi gerektirir. Başarılı yaklaşım, \(A(n)\)'yi çarpımsal bir aritmetik fonksiyona dönüştürmek ve asal kuvvet katkılarını kapalı biçimde hesaplamaktır.

Matematiksel Yaklaşım

Ana fikir, sekiz değişkenli kongruansı önce dört değişkenli bir artık dağılımının ikinci momenti olarak yazmaktır. Sonrasında Çin kalan teoremi problemi ayırır, kuadratik Gauss toplamları da yerel faktörleri verir.

Adım 1: Kongruansı artık dağılımının karesi olarak yaz

Sabit bir \(n\) için

$$r_n(t)=\#\left\{(a,b,c,d)\in(\mathbb{Z}/n\mathbb{Z})^4:\ a^2+b^2+c^2+d^2\equiv t\pmod n\right\}$$

olsun. Aynı artık \(t\)'yi üreten iki dörtlü, birlikte tam bir geçerli sekizli oluşturur. Bu nedenle

$$A(n)=\sum_{t\bmod n} r_n(t)^2.$$

Böylece asıl sayım, modulo \(n\) dört kare toplamlarının dağılımının ikinci momentine dönüşür.

Adım 2: Çin kalan teoremi ile çarpımsallığı kanıtla

\(\gcd(m,n)=1\) olsun. O zaman modulo \(mn\) bir artık sınıfı, modulo \(m\) ve modulo \(n\) iki bağımsız artık sınıfına bire bir karşılık gelir. Bu özdeşleştirme altında modulo \(mn\) bir dörtlü de modulo \(m\) ve modulo \(n\) iki bağımsız dörtlüye ayrılır. Dolayısıyla

$$r_{mn}(t_m,t_n)=r_m(t_m)\,r_n(t_n).$$

Karesini alıp tüm artık çiftleri üzerinde toplayınca

$$A(mn)=A(m)\,A(n)$$

elde edilir. Demek ki yalnızca asal kuvvetler için \(A(p^e)\)'yi bulmak yeterlidir.

Adım 3: Sayımı kuadratik Gauss toplamlarına indir

Şu toplamsal karakterleri tanımlayalım:

$$e_n(x)=\exp\!\left(\frac{2\pi i x}{n}\right),\qquad G_n(u)=\sum_{x\bmod n} e_n(ux^2).$$

Toplamsal karakterlerin ortogonalliğinden

$$r_n(t)=\frac{1}{n}\sum_{u\bmod n} G_n(u)^4\,e_n(-ut)$$

gelir. Şimdi \(t\) üzerinde ikinci momenti alalım. Sonlu Fourier dönüşümü için Parseval özdeşliği

$$A(n)=\sum_{t\bmod n} r_n(t)^2=\frac{1}{n}\sum_{u\bmod n}\lvert G_n(u)\rvert^8$$

sonucunu verir. Kritik sadeleşme budur: sekizlileri doğrudan saymak yerine yerel kuadratik Gauss toplamlarının büyüklüklerini hesaplamak yeterli olur.

Adım 4: Tek sayılı asallar için asal kuvvet faktörü

\(n=p^e\) ve \(p\) tek asal olsun. \(u=p^jv\) biçiminde, \(0\le j<e\) ve \(p\nmid v\) olacak şekilde yazalım. O zaman

$$G_{p^e}(u)=p^j G_{p^{e-j}}(v),\qquad \lvert G_{p^{e-j}}(v)\rvert=p^{(e-j)/2}.$$

Dolayısıyla \(j\) değerlemeli her frekansın katkısı

$$\lvert G_{p^e}(u)\rvert^8=p^{4e+4j}$$

olur. Böyle frekansların sayısı \((p-1)p^{e-j-1}\)'dir; ayrıca \(u=0\) terimi \(p^{8e}\) getirir. Bunları yerine koyunca

$$A(p^e)=p^{7e}+(p-1)\sum_{j=0}^{e-1} p^{4e+3j-1}$$

elde edilir. Aynı ifade yinelemeli biçimde

$$A(1)=1,\qquad A(p^e)=p^7A(p^{e-1})+(p-1)p^{4e-1}$$

şeklindedir. Uygulamanın tek sayılı asal kuvvetler için kullandığı güncelleme tam olarak budur.

Adım 5: \(2\) kuvvetlerini ayrı ele al

\(2^e\) için ayrı bir dal gerekir; çünkü ikinin kuvvetleri üzerindeki kuadratik Gauss toplamları farklı davranır. \(m\ge 2\) iken modulo \(2^m\) tek frekanslar için

$$\lvert G_{2^m}(v)\rvert=2^{(m+1)/2}$$

olur; ama geriye kalan modulo \(2\) durumu sıfırlanır. Bu nedenle yalnızca \(u=2^jv\) ve \(v\) tek iken \(j=0,1,\dots,e-2\) değerlemeleri katkı verir. Yerel faktör

$$A(2^e)=2^{7e}+\sum_{j=0}^{e-2} 2^{4e+3+3j}$$

olur. Kodun kullandığı yinelemeli biçim ise

$$A(2)=2^7,\qquad A(2^e)=2^7A(2^{e-1})+2^{4e+3}\quad (e\ge 2).$$

Adım 6: Genel modül için yerel faktörleri birleştir

Eğer

$$n=\prod_{i=1}^{k} p_i^{e_i}$$

ise, çarpımsallıktan

$$A(n)=\prod_{i=1}^{k} A(p_i^{e_i})$$

çıkar. Böylece tüm problem

$$\boxed{S(12\,345\,678)=\sum_{n=1}^{12\,345\,678}\prod_{p^e\parallel n}A(p^e)\pmod{1\,001\,961\,001}.}$$

hesabına indirgenmiş olur.

Çözümlü Örnek: \(n=10\) ve \(10\)'a kadar kontrol toplamı

\(10=2\cdot 5\) olduğundan

$$A(10)=A(2)\,A(5)$$

yazılır. Yerel formüllerden

$$A(2)=2^7=128,$$

$$A(5)=5^7+4\cdot 5^3=78\,625$$

bulunur. O halde

$$A(10)=128\cdot 78\,625=10\,064\,000.$$

Aynı formüllerle tüm \(n\le 10\) için şu değerler gelir:

$$1,\ 128,\ 2\,241,\ 18\,432,\ 78\,625,\ 286\,848,\ 825\,601,\ 2\,392\,064,\ 4\,905\,441,\ 10\,064\,000,$$

dolayısıyla

$$S(10)=18\,573\,381.$$

Bu, uygulamanın doğruladığı küçük kontrol değeridir.

Kod Nasıl Çalışır

C++, Python ve Java uygulamaları aynı akışı izler. Önce \(12\,345\,678\)'e kadar en küçük asal bölen tablosu kurulur. Bu tablo sayesinde her \(n\) hızlı bir tarama içinde çarpanlarına ayrılabilir.

Her \(n\) için eşit asal çarpanlar gruplanarak asal kuvvet ayrışımı çıkarılır. Sonra her yerel faktör \(A(p^e)\), \(1\,001\,961\,001\) modunda hesaplanır: tek asallar Adım 4'teki yinelemeyi, \(2^e\) ise Adım 5'teki özel \(2\)-adik yinelemeyi kullanır.

Yerel faktörler çarpılarak \(A(n)\bmod 1\,001\,961\,001\) elde edilir ve ardından

$$\sum_{n=1}^{12\,345\,678} A(n)\pmod{1\,001\,961\,001}$$

koşan toplamda biriktirilir. C++ uygulaması tam çalışmadan önce \(A(4)=18\,432\) ve \(S(10)=18\,573\,381\) gibi küçük kesin kontroller de yapar.

Karmaşıklık Analizi

\(N=12\,345\,678\) olsun. Doğrusal elek ile en küçük asal bölen tablosunu kurmak \(O(N)\) zaman ve \(O(N)\) bellek gerektirir. Bu tablodan tüm \(n\le N\) değerlerini çarpanlarına ayırmanın toplam maliyeti, karşılaşılan toplam asal çarpan sayısıyla orantılıdır ve ortalamada \(O(N\log\log N)\) düzeyindedir. Sonuç olarak yöntem pratikte neredeyse doğrusaldır ve \(O(N)\) bellek kullanır.

Dipnotlar ve Kaynaklar

  1. Problem sayfası: Project Euler 875
  2. Çin kalan teoremi: Wikipedia — Chinese remainder theorem
  3. Gauss toplamı: Wikipedia — Gauss sum
  4. Parseval özdeşliği: Wikipedia — Parseval's theorem
  5. Dört kare teoremi: Wikipedia — Lagrange's four-square theorem

Resumen del Problema

Para cada entero positivo \(n\), definimos

$$A(n)=\#\left\{(x_1,\dots,x_8)\in(\mathbb{Z}/n\mathbb{Z})^8:\ x_1^2+x_2^2+x_3^2+x_4^2 \equiv x_5^2+x_6^2+x_7^2+x_8^2 \pmod n\right\}.$$

La cantidad pedida es

$$S(N)=\sum_{n=1}^{N} A(n),\qquad S(12\,345\,678)\bmod 1\,001\,961\,001.$$

Una enumeración directa exigiría revisar cantidades enormes de octuplas de residuos para cada módulo. La solución útil consiste en convertir \(A(n)\) en una función aritmética multiplicativa y calcular explícitamente sus factores en potencias primas.

Enfoque Matemático

La observación central es que la congruencia de ocho variables puede reescribirse como el segundo momento de una distribución de residuos con cuatro variables. Después, el teorema chino del resto y las sumas de Gauss cuadráticas hacen el resto.

Paso 1: Reescribir la congruencia como un cuadrado de distribución

Para un módulo fijo \(n\), definimos

$$r_n(t)=\#\left\{(a,b,c,d)\in(\mathbb{Z}/n\mathbb{Z})^4:\ a^2+b^2+c^2+d^2\equiv t\pmod n\right\}.$$

Si dos cuádruplas producen el mismo residuo \(t\), juntas forman una octupla válida. Por tanto,

$$A(n)=\sum_{t\bmod n} r_n(t)^2.$$

Así, el problema original se convierte en el segundo momento de la distribución de las sumas de cuatro cuadrados módulo \(n\).

Paso 2: Probar la multiplicatividad con el teorema chino del resto

Supongamos \(\gcd(m,n)=1\). Un residuo módulo \(mn\) se corresponde de manera única con un par de residuos módulo \(m\) y módulo \(n\). Bajo esa identificación, una cuádrupla módulo \(mn\) es simplemente un par independiente de cuádruplas, de modo que

$$r_{mn}(t_m,t_n)=r_m(t_m)\,r_n(t_n).$$

Al elevar al cuadrado y sumar sobre todos los pares de residuos se obtiene

$$A(mn)=A(m)\,A(n).$$

Por eso basta estudiar \(A(p^e)\) para potencias primas y luego multiplicar los factores locales.

Paso 3: Pasar a sumas de Gauss cuadráticas

Introducimos los caracteres aditivos

$$e_n(x)=\exp\!\left(\frac{2\pi i x}{n}\right),\qquad G_n(u)=\sum_{x\bmod n} e_n(ux^2).$$

La ortogonalidad de los caracteres aditivos da

$$r_n(t)=\frac{1}{n}\sum_{u\bmod n} G_n(u)^4\,e_n(-ut).$$

Ahora tomamos el segundo momento en \(t\). La identidad de Parseval para la transformada discreta finita produce

$$A(n)=\sum_{t\bmod n} r_n(t)^2=\frac{1}{n}\sum_{u\bmod n}\lvert G_n(u)\rvert^8.$$

Esta es la reducción decisiva: en lugar de contar octuplas directamente, basta con conocer los tamaños de las sumas de Gauss locales.

Paso 4: Evaluar el factor de potencia prima impar

Sea \(n=p^e\) con \(p\) impar. Escribimos \(u=p^jv\), donde \(0\le j<e\) y \(p\nmid v\). Entonces

$$G_{p^e}(u)=p^j G_{p^{e-j}}(v),\qquad \lvert G_{p^{e-j}}(v)\rvert=p^{(e-j)/2}.$$

Por tanto, toda frecuencia con valoración \(j\) aporta

$$\lvert G_{p^e}(u)\rvert^8=p^{4e+4j}.$$

Hay \((p-1)p^{e-j-1}\) frecuencias de ese tipo, mientras que \(u=0\) aporta \(p^{8e}\). Sustituyendo en la fórmula principal se obtiene

$$A(p^e)=p^{7e}+(p-1)\sum_{j=0}^{e-1} p^{4e+3j-1}.$$

De forma recursiva, esto equivale a

$$A(1)=1,\qquad A(p^e)=p^7A(p^{e-1})+(p-1)p^{4e-1}.$$

Esa es exactamente la regla local que usan las implementaciones para potencias primas impares.

Paso 5: Tratar por separado las potencias de \(2\)

El caso \(2^e\) necesita una rama distinta porque las sumas de Gauss cuadráticas sobre potencias de \(2\) tienen un tamaño diferente. Para frecuencia impar módulo \(2^m\) con \(m\ge 2\), se tiene

$$\lvert G_{2^m}(v)\rvert=2^{(m+1)/2},$$

mientras que el caso residual con módulo \(2\) se anula. Por eso solo contribuyen las valoraciones \(j=0,1,\dots,e-2\) cuando \(u=2^jv\) con \(v\) impar. El factor local queda

$$A(2^e)=2^{7e}+\sum_{j=0}^{e-2} 2^{4e+3+3j}.$$

En la forma iterativa usada por el código, esto se escribe como

$$A(2)=2^7,\qquad A(2^e)=2^7A(2^{e-1})+2^{4e+3}\quad (e\ge 2).$$

Paso 6: Combinar los factores locales para un módulo general

Si

$$n=\prod_{i=1}^{k} p_i^{e_i},$$

entonces la multiplicatividad implica

$$A(n)=\prod_{i=1}^{k} A(p_i^{e_i}).$$

Así, todo el problema se reduce a

$$\boxed{S(12\,345\,678)=\sum_{n=1}^{12\,345\,678}\prod_{p^e\parallel n}A(p^e)\pmod{1\,001\,961\,001}.}$$

Ejemplo Resuelto: \(n=10\) y la comprobación acumulada hasta \(10\)

Como \(10=2\cdot 5\), la multiplicatividad da inmediatamente

$$A(10)=A(2)\,A(5).$$

De las fórmulas locales obtenemos

$$A(2)=2^7=128,$$

$$A(5)=5^7+4\cdot 5^3=78\,625.$$

Por lo tanto,

$$A(10)=128\cdot 78\,625=10\,064\,000.$$

Aplicando las mismas fórmulas a todos los \(n\le 10\) se obtienen

$$1,\ 128,\ 2\,241,\ 18\,432,\ 78\,625,\ 286\,848,\ 825\,601,\ 2\,392\,064,\ 4\,905\,441,\ 10\,064\,000,$$

y en consecuencia

$$S(10)=18\,573\,381.$$

Ese es el pequeño punto de control que verifica la implementación.

Cómo Funciona el Código

Las implementaciones en C++, Python y Java siguen la misma estrategia. Primero construyen una tabla del menor factor primo hasta \(12\,345\,678\). Esa tabla permite factorizar cada \(n\) mediante un recorrido rápido.

Para cada \(n\), el programa agrupa factores primos iguales para reconstruir la descomposición en potencias primas. Después evalúa cada factor local \(A(p^e)\) módulo \(1\,001\,961\,001\): las potencias de primos impares usan la recurrencia del Paso 4, y la potencia \(2^e\) usa la recurrencia especial del Paso 5.

Los factores locales se multiplican para producir \(A(n)\bmod 1\,001\,961\,001\), y luego se acumula

$$\sum_{n=1}^{12\,345\,678} A(n)\pmod{1\,001\,961\,001}.$$

La implementación en C++ también incluye pequeñas comprobaciones exactas antes del recorrido completo, entre ellas \(A(4)=18\,432\) y \(S(10)=18\,573\,381\).

Análisis de Complejidad

Sea \(N=12\,345\,678\). Construir la tabla del menor factor primo con una criba lineal cuesta \(O(N)\) tiempo y \(O(N)\) memoria. Factorizar todos los enteros hasta \(N\) a partir de esa tabla requiere trabajo total proporcional al número total de factores primos encontrados, lo que es \(O(N\log\log N)\) en promedio. En conjunto, el método es casi lineal en la práctica y utiliza \(O(N)\) memoria.

Notas y Referencias

  1. Página del problema: Project Euler 875
  2. Teorema chino del resto: Wikipedia — Chinese remainder theorem
  3. Suma de Gauss: Wikipedia — Gauss sum
  4. Teorema de Parseval: Wikipedia — Parseval's theorem
  5. Teorema de los cuatro cuadrados: Wikipedia — Lagrange's four-square theorem

问题概述

对每个正整数 \(n\),定义

$$A(n)=\#\left\{(x_1,\dots,x_8)\in(\mathbb{Z}/n\mathbb{Z})^8:\ x_1^2+x_2^2+x_3^2+x_4^2 \equiv x_5^2+x_6^2+x_7^2+x_8^2 \pmod n\right\}.$$

题目要求计算

$$S(N)=\sum_{n=1}^{N} A(n),\qquad S(12\,345\,678)\bmod 1\,001\,961\,001.$$

如果直接枚举,每个模数都要检查海量的八元组剩余类,完全不可行。真正有效的方法是把 \(A(n)\) 改写成一个乘法性算术函数,再把每个素数幂的局部贡献单独算出来。

数学方法

核心思路是先把八变量同余改写成四平方剩余分布的二次矩,然后再利用中国剩余定理与二次 Gauss 和,把全局计数拆成局部素数幂问题。

步骤 1:把原问题改写成剩余分布的平方和

固定模数 \(n\),定义

$$r_n(t)=\#\left\{(a,b,c,d)\in(\mathbb{Z}/n\mathbb{Z})^4:\ a^2+b^2+c^2+d^2\equiv t\pmod n\right\}.$$

如果两个四元组产生相同的剩余 \(t\),那么它们拼在一起正好构成一个满足条件的八元组。因此

$$A(n)=\sum_{t\bmod n} r_n(t)^2.$$

也就是说,原来的八变量计数,等价于“模 \(n\) 下四平方和的剩余分布”的二次矩。

步骤 2:利用中国剩余定理证明乘法性

设 \(\gcd(m,n)=1\)。根据中国剩余定理,模 \(mn\) 的一个剩余类与模 \(m\) 和模 \(n\) 的一对剩余类一一对应。在这种对应下,模 \(mn\) 的一个四元组恰好等价于模 \(m\) 与模 \(n\) 的两个独立四元组,所以

$$r_{mn}(t_m,t_n)=r_m(t_m)\,r_n(t_n).$$

对所有剩余对平方并求和,立刻得到

$$A(mn)=A(m)\,A(n).$$

因此只要弄清楚每个素数幂 \(p^e\) 的值 \(A(p^e)\),一般 \(n\) 的结果就能通过乘法性恢复出来。

步骤 3:把计数转换成二次 Gauss 和

引入加性特征

$$e_n(x)=\exp\!\left(\frac{2\pi i x}{n}\right),\qquad G_n(u)=\sum_{x\bmod n} e_n(ux^2).$$

由加性特征的正交性可得

$$r_n(t)=\frac{1}{n}\sum_{u\bmod n} G_n(u)^4\,e_n(-ut).$$

接着对 \(t\) 取二次矩。有限 Fourier 变换的 Parseval 恒等式给出

$$A(n)=\sum_{t\bmod n} r_n(t)^2=\frac{1}{n}\sum_{u\bmod n}\lvert G_n(u)\rvert^8.$$

这是整个推导中最关键的一步:我们不再直接数八元组,而是只需要知道不同频率上的二次 Gauss 和大小。

步骤 4:求奇素数幂的局部因子

设 \(n=p^e\),其中 \(p\) 是奇素数。把频率写成 \(u=p^jv\),其中 \(0\le j<e\) 且 \(p\nmid v\)。则有

$$G_{p^e}(u)=p^j G_{p^{e-j}}(v),\qquad \lvert G_{p^{e-j}}(v)\rvert=p^{(e-j)/2}.$$

于是 valuation 为 \(j\) 的每个频率都会贡献

$$\lvert G_{p^e}(u)\rvert^8=p^{4e+4j}.$$

这类频率共有 \((p-1)p^{e-j-1}\) 个,而 \(u=0\) 单独贡献 \(p^{8e}\)。代回主公式后得到

$$A(p^e)=p^{7e}+(p-1)\sum_{j=0}^{e-1} p^{4e+3j-1}.$$

把它写成递推式会更贴近实现:

$$A(1)=1,\qquad A(p^e)=p^7A(p^{e-1})+(p-1)p^{4e-1}.$$

这正是实现中处理奇素数幂时使用的更新规律。

步骤 5:单独处理 \(2\) 的幂

\(2^e\) 需要单独分支,因为二次 Gauss 和在二的幂模数下的行为与奇素数不同。对模 \(2^m\) 的奇频率且 \(m\ge 2\),有

$$\lvert G_{2^m}(v)\rvert=2^{(m+1)/2},$$

而当剩下的模数退化成 \(2\) 时,对应项会消失。因此只有 valuation \(j=0,1,\dots,e-2\) 的频率 \(u=2^jv\)(其中 \(v\) 为奇数)会有贡献。于是局部因子变成

$$A(2^e)=2^{7e}+\sum_{j=0}^{e-2} 2^{4e+3+3j}.$$

把它写成程序最容易使用的递推形式,就是

$$A(2)=2^7,\qquad A(2^e)=2^7A(2^{e-1})+2^{4e+3}\quad (e\ge 2).$$

步骤 6:把一般模数分解为局部因子的乘积

如果

$$n=\prod_{i=1}^{k} p_i^{e_i},$$

那么乘法性立刻给出

$$A(n)=\prod_{i=1}^{k} A(p_i^{e_i}).$$

于是整个问题最终化成

$$\boxed{S(12\,345\,678)=\sum_{n=1}^{12\,345\,678}\prod_{p^e\parallel n}A(p^e)\pmod{1\,001\,961\,001}.}$$

计算示例:\(n=10\) 以及到 \(10\) 为止的校验值

由于 \(10=2\cdot 5\),由乘法性可知

$$A(10)=A(2)\,A(5).$$

根据局部公式,

$$A(2)=2^7=128,$$

$$A(5)=5^7+4\cdot 5^3=78\,625.$$

因此

$$A(10)=128\cdot 78\,625=10\,064\,000.$$

同样的方法可算出所有 \(n\le 10\) 的值:

$$1,\ 128,\ 2\,241,\ 18\,432,\ 78\,625,\ 286\,848,\ 825\,601,\ 2\,392\,064,\ 4\,905\,441,\ 10\,064\,000,$$

从而

$$S(10)=18\,573\,381.$$

这正是实现里用于确认推导无误的小型校验点。

代码如何工作

C++、Python 和 Java 实现遵循同一条主线。首先构建到 \(12\,345\,678\) 为止的最小质因子表,这样每个 \(n\) 都可以在顺序扫描时快速分解。

对于每个 \(n\),实现会把相同的质因子连续剥离,恢复出素数幂分解。然后对每个局部因子 \(A(p^e)\) 在模 \(1\,001\,961\,001\) 下计算:奇素数使用步骤 4 的递推,\(2^e\) 使用步骤 5 的特殊 \(2\)-进递推。

所有局部因子相乘得到 \(A(n)\bmod 1\,001\,961\,001\),再把它累加到

$$\sum_{n=1}^{12\,345\,678} A(n)\pmod{1\,001\,961\,001}$$

这个总和中。C++ 实现还会在完整运行前做很小的精确自检,例如验证 \(A(4)=18\,432\) 与 \(S(10)=18\,573\,381\)。

复杂度分析

记 \(N=12\,345\,678\)。用线性筛建立最小质因子表需要 \(O(N)\) 时间和 \(O(N)\) 内存。随后用这张表分解所有 \(n\le N\) 的总代价,与遇到的质因子总数成正比,平均规模为 \(O(N\log\log N)\)。因此整个方法在实践中接近线性,并且只需要 \(O(N)\) 内存。

注释与参考资料

  1. 题目页面: Project Euler 875
  2. 中国剩余定理: Wikipedia — Chinese remainder theorem
  3. Gauss 和: Wikipedia — Gauss sum
  4. Parseval 定理: Wikipedia — Parseval's theorem
  5. 四平方定理: Wikipedia — Lagrange's four-square theorem

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

Для каждого положительного целого \(n\) определим

$$A(n)=\#\left\{(x_1,\dots,x_8)\in(\mathbb{Z}/n\mathbb{Z})^8:\ x_1^2+x_2^2+x_3^2+x_4^2 \equiv x_5^2+x_6^2+x_7^2+x_8^2 \pmod n\right\}.$$

Нужно вычислить

$$S(N)=\sum_{n=1}^{N} A(n),\qquad S(12\,345\,678)\bmod 1\,001\,961\,001.$$

Прямой перебор здесь бесполезен: для каждого модуля пришлось бы рассматривать огромное число восьмёрок вычетов. Рабочий путь состоит в том, чтобы превратить \(A(n)\) в мультипликативную арифметическую функцию и явно вычислять её локальные множители на степенях простых чисел.

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

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

Шаг 1: Перепишем задачу как квадрат распределения вычетов

Для фиксированного \(n\) положим

$$r_n(t)=\#\left\{(a,b,c,d)\in(\mathbb{Z}/n\mathbb{Z})^4:\ a^2+b^2+c^2+d^2\equiv t\pmod n\right\}.$$

Если две четвёрки дают один и тот же вычет \(t\), вместе они образуют допустимую восьмёрку. Поэтому

$$A(n)=\sum_{t\bmod n} r_n(t)^2.$$

Итак, исходный счёт сводится ко второму моменту распределения вычетов суммы четырёх квадратов по модулю \(n\).

Шаг 2: Докажем мультипликативность через китайскую теорему об остатках

Пусть \(\gcd(m,n)=1\). Тогда класс вычетов по модулю \(mn\) однозначно соответствует паре классов по модулям \(m\) и \(n\). При этом четвёрка по модулю \(mn\) распадается на независимую пару четвёрок, так что

$$r_{mn}(t_m,t_n)=r_m(t_m)\,r_n(t_n).$$

После возведения в квадрат и суммирования по всем парам вычетов получаем

$$A(mn)=A(m)\,A(n).$$

Значит, достаточно найти \(A(p^e)\) для степеней простых чисел, а затем перемножить локальные множители.

Шаг 3: Перейдём к квадратичным суммам Гаусса

Введём аддитивные характеры

$$e_n(x)=\exp\!\left(\frac{2\pi i x}{n}\right),\qquad G_n(u)=\sum_{x\bmod n} e_n(ux^2).$$

Из ортогональности аддитивных характеров следует

$$r_n(t)=\frac{1}{n}\sum_{u\bmod n} G_n(u)^4\,e_n(-ut).$$

Теперь берём второй момент по \(t\). Формула Парсеваля для конечного преобразования Фурье даёт

$$A(n)=\sum_{t\bmod n} r_n(t)^2=\frac{1}{n}\sum_{u\bmod n}\lvert G_n(u)\rvert^8.$$

Это ключевое упрощение: вместо прямого подсчёта восьмёрок нужно знать только величины локальных квадратичных сумм Гаусса.

Шаг 4: Вычислим множитель для нечётной степени простого

Пусть \(n=p^e\), где \(p\) нечётное простое. Запишем \(u=p^jv\), где \(0\le j<e\) и \(p\nmid v\). Тогда

$$G_{p^e}(u)=p^j G_{p^{e-j}}(v),\qquad \lvert G_{p^{e-j}}(v)\rvert=p^{(e-j)/2}.$$

Следовательно, каждая частота с valuation \(j\) даёт вклад

$$\lvert G_{p^e}(u)\rvert^8=p^{4e+4j}.$$

Таких частот \((p-1)p^{e-j-1}\), а член \(u=0\) даёт \(p^{8e}\). Подставляя это в основную формулу, получаем

$$A(p^e)=p^{7e}+(p-1)\sum_{j=0}^{e-1} p^{4e+3j-1}.$$

В рекурсивной форме это выглядит так:

$$A(1)=1,\qquad A(p^e)=p^7A(p^{e-1})+(p-1)p^{4e-1}.$$

Именно это локальное обновление используется реализациями для нечётных простых степеней.

Шаг 5: Отдельно обработаем степени двойки

Для \(2^e\) нужна отдельная ветвь, потому что квадратичные суммы Гаусса по модулям \(2^m\) имеют другой размер. Для нечётной частоты по модулю \(2^m\) при \(m\ge 2\) выполняется

$$\lvert G_{2^m}(v)\rvert=2^{(m+1)/2},$$

а при оставшемся модуле \(2\) соответствующий случай зануляется. Поэтому вклад дают только valuation \(j=0,1,\dots,e-2\) у частот вида \(u=2^jv\), где \(v\) нечётно. Получаем

$$A(2^e)=2^{7e}+\sum_{j=0}^{e-2} 2^{4e+3+3j}.$$

В виде рекурсии, удобном для программы, это

$$A(2)=2^7,\qquad A(2^e)=2^7A(2^{e-1})+2^{4e+3}\quad (e\ge 2).$$

Шаг 6: Восстановим значение для общего модуля

Если

$$n=\prod_{i=1}^{k} p_i^{e_i},$$

то по мультипликативности

$$A(n)=\prod_{i=1}^{k} A(p_i^{e_i}).$$

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

$$\boxed{S(12\,345\,678)=\sum_{n=1}^{12\,345\,678}\prod_{p^e\parallel n}A(p^e)\pmod{1\,001\,961\,001}.}$$

Разобранный пример: \(n=10\) и контрольная сумма до \(10\)

Так как \(10=2\cdot 5\), по мультипликативности сразу имеем

$$A(10)=A(2)\,A(5).$$

Из локальных формул следует

$$A(2)=2^7=128,$$

$$A(5)=5^7+4\cdot 5^3=78\,625.$$

Значит,

$$A(10)=128\cdot 78\,625=10\,064\,000.$$

Те же формулы для всех \(n\le 10\) дают значения

$$1,\ 128,\ 2\,241,\ 18\,432,\ 78\,625,\ 286\,848,\ 825\,601,\ 2\,392\,064,\ 4\,905\,441,\ 10\,064\,000,$$

и потому

$$S(10)=18\,573\,381.$$

Именно этот небольшой контрольный результат подтверждается реализацией.

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

Реализации на C++, Python и Java используют одну и ту же схему. Сначала строится таблица наименьших простых делителей до \(12\,345\,678\). Благодаря ей каждое \(n\) можно быстро разложить на множители в одном последовательном проходе.

Для каждого \(n\) программа снимает одинаковые простые множители и получает разложение на простые степени. Затем каждый локальный множитель \(A(p^e)\) вычисляется по модулю \(1\,001\,961\,001\): нечётные простые используют рекурсию из шага 4, а \(2^e\) обрабатывается специальной \(2\)-адической рекурсией из шага 5.

Локальные множители перемножаются, чтобы получить \(A(n)\bmod 1\,001\,961\,001\), после чего накапливается сумма

$$\sum_{n=1}^{12\,345\,678} A(n)\pmod{1\,001\,961\,001}.$$

Реализация на C++ также содержит маленькие точные самопроверки перед полным запуском, включая \(A(4)=18\,432\) и \(S(10)=18\,573\,381\).

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

Пусть \(N=12\,345\,678\). Построение таблицы наименьших простых делителей линейным решетом требует \(O(N)\) времени и \(O(N)\) памяти. Факторизация всех чисел до \(N\) с использованием этой таблицы требует суммарной работы, пропорциональной общему числу встреченных простых множителей, то есть в среднем \(O(N\log\log N)\). В итоге метод практически квазилинеен и использует \(O(N)\) памяти.

Примечания и ссылки

  1. Страница задачи: Project Euler 875
  2. Китайская теорема об остатках: Wikipedia — Chinese remainder theorem
  3. Сумма Гаусса: Wikipedia — Gauss sum
  4. Теорема Парсеваля: Wikipedia — Parseval's theorem
  5. Теорема о четырёх квадратах: Wikipedia — Lagrange's four-square theorem

ملخص المسألة

لكل عدد صحيح موجب \(n\)، نعرّف

$$A(n)=\#\left\{(x_1,\dots,x_8)\in(\mathbb{Z}/n\mathbb{Z})^8:\ x_1^2+x_2^2+x_3^2+x_4^2 \equiv x_5^2+x_6^2+x_7^2+x_8^2 \pmod n\right\}.$$

والمطلوب هو حساب

$$S(N)=\sum_{n=1}^{N} A(n),\qquad S(12\,345\,678)\bmod 1\,001\,961\,001.$$

العد المباشر غير عملي تمامًا، لأن كل قيمة لـ \(n\) ستتطلب فحص عدد هائل من الثمانيات من الأصناف التوافقية. الفكرة الصحيحة هي تحويل \(A(n)\) إلى دالة حسابية ضربية، ثم حساب مساهمة كل قوة أولية على حدة.

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

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

الخطوة 1: إعادة كتابة التوافق على شكل مربع لتوزيع البواقي

للمودولو الثابت \(n\)، نعرّف

$$r_n(t)=\#\left\{(a,b,c,d)\in(\mathbb{Z}/n\mathbb{Z})^4:\ a^2+b^2+c^2+d^2\equiv t\pmod n\right\}.$$

إذا أعطت رباعيتان الباقي نفسه \(t\)، فإنهما معًا تشكلان ثمانية صحيحة تحقق الشرط. لذلك

$$A(n)=\sum_{t\bmod n} r_n(t)^2.$$

وهكذا يتحول العد الأصلي إلى العزم الثاني لتوزيع بواقي مجموع أربعة مربعات modulo \(n\).

الخطوة 2: إثبات الضربية باستخدام مبرهنة الباقي الصيني

إذا كان \(\gcd(m,n)=1\)، فإن كل فئة توافقية modulo \(mn\) تقابل زوجًا وحيدًا من الفئات modulo \(m\) وmodulo \(n\). وتحت هذا التطابق، تصبح كل رباعية modulo \(mn\) مجرد زوج مستقل من رباعيتين، ولهذا

$$r_{mn}(t_m,t_n)=r_m(t_m)\,r_n(t_n).$$

وبعد التربيع والجمع على جميع أزواج البواقي نحصل على

$$A(mn)=A(m)\,A(n).$$

إذن يكفي حساب \(A(p^e)\) لقوى الأعداد الأولية، ثم ضرب العوامل المحلية معًا.

الخطوة 3: تحويل العد إلى مجاميع غاوس التربيعية

ندخل المحارف الجمعية

$$e_n(x)=\exp\!\left(\frac{2\pi i x}{n}\right),\qquad G_n(u)=\sum_{x\bmod n} e_n(ux^2).$$

وبفضل تعامد المحارف الجمعية نحصل على

$$r_n(t)=\frac{1}{n}\sum_{u\bmod n} G_n(u)^4\,e_n(-ut).$$

ثم نأخذ العزم الثاني بالنسبة إلى \(t\). تعطي هوية بارسيفال لتحويل فورييه المنتهي

$$A(n)=\sum_{t\bmod n} r_n(t)^2=\frac{1}{n}\sum_{u\bmod n}\lvert G_n(u)\rvert^8.$$

وهذه هي الخطوة الحاسمة: بدلًا من عد الثمانيات مباشرة، يكفي معرفة أحجام مجاميع غاوس التربيعية المحلية.

الخطوة 4: حساب العامل المحلي لقوى الأوليات الفردية

لنأخذ \(n=p^e\) حيث \(p\) أولي فردي. نكتب \(u=p^jv\) بحيث \(0\le j<e\) و\(p\nmid v\). عندئذٍ

$$G_{p^e}(u)=p^j G_{p^{e-j}}(v),\qquad \lvert G_{p^{e-j}}(v)\rvert=p^{(e-j)/2}.$$

وبالتالي فإن كل تردد ذي تقييم \(j\) يساهم بالمقدار

$$\lvert G_{p^e}(u)\rvert^8=p^{4e+4j}.$$

وعدد هذه الترددات هو \((p-1)p^{e-j-1}\)، بينما يعطي الحد \(u=0\) مساهمة \(p^{8e}\). بالتعويض في الصيغة الرئيسية نحصل على

$$A(p^e)=p^{7e}+(p-1)\sum_{j=0}^{e-1} p^{4e+3j-1}.$$

وبالصيغة التراجعية المكافئة،

$$A(1)=1,\qquad A(p^e)=p^7A(p^{e-1})+(p-1)p^{4e-1}.$$

وهذه بالضبط هي قاعدة التحديث التي تستخدمها التطبيقات لقوى الأوليات الفردية.

الخطوة 5: معالجة قوى العدد \(2\) على نحو منفصل

تحتاج \(2^e\) إلى فرع مستقل، لأن مجاميع غاوس التربيعية على قوى \(2\) لا تتصرف مثل الحالة الفردية. فبالنسبة إلى تردد فردي modulo \(2^m\) مع \(m\ge 2\)، لدينا

$$\lvert G_{2^m}(v)\rvert=2^{(m+1)/2},$$

أما الحالة المتبقية عندما يصبح المودولو \(2\) فتختفي. لذلك تسهم فقط التقييمات \(j=0,1,\dots,e-2\) حين يكون \(u=2^jv\) و\(v\) فرديًا. ومن ثم يصبح العامل المحلي

$$A(2^e)=2^{7e}+\sum_{j=0}^{e-2} 2^{4e+3+3j}.$$

وفي الصيغة التراجعية التي تستخدمها الشفرة،

$$A(2)=2^7,\qquad A(2^e)=2^7A(2^{e-1})+2^{4e+3}\quad (e\ge 2).$$

الخطوة 6: تركيب القيمة العامة من العوامل المحلية

إذا كان

$$n=\prod_{i=1}^{k} p_i^{e_i},$$

فإن الضربية تعطي مباشرة

$$A(n)=\prod_{i=1}^{k} A(p_i^{e_i}).$$

وبذلك تختزل المسألة كلها إلى

$$\boxed{S(12\,345\,678)=\sum_{n=1}^{12\,345\,678}\prod_{p^e\parallel n}A(p^e)\pmod{1\,001\,961\,001}.}$$

مثال محلول: \(n=10\) ونقطة التحقق حتى \(10\)

بما أن \(10=2\cdot 5\)، فإن الضربية تعطي فورًا

$$A(10)=A(2)\,A(5).$$

ومن الصيغ المحلية نجد

$$A(2)=2^7=128,$$

$$A(5)=5^7+4\cdot 5^3=78\,625.$$

إذن

$$A(10)=128\cdot 78\,625=10\,064\,000.$$

وبتطبيق الصيغ نفسها على جميع \(n\le 10\) نحصل على

$$1,\ 128,\ 2\,241,\ 18\,432,\ 78\,625,\ 286\,848,\ 825\,601,\ 2\,392\,064,\ 4\,905\,441,\ 10\,064\,000,$$

ومن ثم

$$S(10)=18\,573\,381.$$

وهذه هي قيمة التحقق الصغيرة التي تؤكدها الشفرة قبل الحساب الكامل.

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

تتبع تطبيقات C++ وPython وJava البنية نفسها. أولًا تُبنى جدولـة أصغر عامل أولي حتى \(12\,345\,678\)، مما يسمح بتفكيك كل \(n\) بسرعة أثناء المسح التسلسلي.

لكل \(n\)، تُجمع العوامل الأولية المتساوية لاستخراج التحليل إلى قوى أولية. ثم يُحسب كل عامل محلي \(A(p^e)\) modulo \(1\,001\,961\,001\): تستعمل الأوليات الفردية العلاقة التراجعية في الخطوة 4، بينما تستخدم \(2^e\) العلاقة الخاصة في الخطوة 5.

بعد ذلك تُضرب العوامل المحلية للحصول على \(A(n)\bmod 1\,001\,961\,001\)، ثم تُراكم الكمية

$$\sum_{n=1}^{12\,345\,678} A(n)\pmod{1\,001\,961\,001}.$$

كما تتضمن نسخة C++ اختبارات تحقق دقيقة صغيرة قبل التشغيل الكامل، ومنها \(A(4)=18\,432\) و\(S(10)=18\,573\,381\).

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

لنكتب \(N=12\,345\,678\). بناء جدول أصغر عامل أولي باستخدام غربال خطي يكلف \(O(N)\) زمنًا و\(O(N)\) ذاكرة. أما تحليل جميع الأعداد حتى \(N\) اعتمادًا على هذا الجدول فيتطلب عملًا إجماليًا يتناسب مع العدد الكلي للعوامل الأولية الظاهرة، أي \(O(N\log\log N)\) في المتوسط. لذلك فالطريقة شبه خطية عمليًا وتستخدم \(O(N)\) ذاكرة.

حواشٍ ومراجع

  1. صفحة المسألة: Project Euler 875
  2. مبرهنة الباقي الصيني: Wikipedia — Chinese remainder theorem
  3. مجموع غاوس: Wikipedia — Gauss sum
  4. هوية بارسيفال: Wikipedia — Parseval's theorem
  5. مبرهنة المربعات الأربعة: Wikipedia — Lagrange's four-square theorem