Problem Summary

Define

$$G(N)=\sum_{j=1}^{N}\sum_{i=1}^{j}\gcd(i,j).$$

The problem asks for \(G(10^{11}) \bmod 998244353\). A direct double loop is impossible at that scale, so the solution rewrites the gcd contribution as a divisor sum, introduces the summatory totient function, and then evaluates the remaining expressions by floor-division blocks.

Mathematical Approach

It is convenient to isolate the inner sum

$$T(j)=\sum_{i=1}^{j}\gcd(i,j),$$

and also to define the totient prefix

$$\Phi(x)=\sum_{k=1}^{x}\varphi(k).$$

The entire method is built around turning \(T(j)\) into something that depends on divisors of \(j\), and then exploiting the fact that floor quotients stay constant on long intervals.

Step 1: Rewrite the inner gcd sum as a divisor sum

Fix \(j\). Group the indices \(i\) by the value \(d=\gcd(i,j)\). If \(d\) is fixed, then we can write

$$i=d\,a,\qquad j=d\,b,\qquad \gcd(a,b)=1.$$

Here \(b=j/d\), and the admissible values of \(a\) are exactly the integers \(1\le a\le b\) that are coprime to \(b\). Their number is \(\varphi(b)\).

So each divisor \(d\mid j\) contributes the value \(d\) exactly \(\varphi(j/d)\) times, which gives

$$T(j)=\sum_{d\mid j} d\,\varphi\!\left(\frac{j}{d}\right).$$

This identity removes the explicit gcd from the inner loop and replaces it with a clean divisor formula.

Step 2: Reverse the order of summation

Now sum \(T(j)\) over all \(j\le N\):

$$G(N)=\sum_{j=1}^{N}\sum_{d\mid j} d\,\varphi\!\left(\frac{j}{d}\right).$$

Write \(j=d\,k\). Every pair \((d,k)\) with \(d\,k\le N\) appears exactly once, so

$$G(N)=\sum_{d=1}^{N} d\sum_{k\le N/d}\varphi(k)=\sum_{d=1}^{N} d\,\Phi\!\left(\left\lfloor\frac{N}{d}\right\rfloor\right).$$

At this point the whole problem is reduced to answering many queries for the prefix sum \(\Phi(x)\).

Step 3: Derive a recurrence for the totient prefix

The classical identity

$$\sum_{d\mid n}\varphi(d)=n$$

holds for every positive integer \(n\). Summing it over \(1\le n\le x\) gives

$$\sum_{n=1}^{x} n=\sum_{n=1}^{x}\sum_{d\mid n}\varphi(d).$$

Swap the order of summation. A fixed \(d\) divides exactly \(\left\lfloor x/d\right\rfloor\) integers up to \(x\), so

$$\frac{x(x+1)}{2}=\sum_{d=1}^{x}\varphi(d)\left\lfloor\frac{x}{d}\right\rfloor.$$

Rewrite the right-hand side by grouping equal quotients. This yields

$$\frac{x(x+1)}{2}=\sum_{m=1}^{x}\Phi\!\left(\left\lfloor\frac{x}{m}\right\rfloor\right).$$

Isolating the \(m=1\) term gives the recurrence

$$\Phi(x)=\frac{x(x+1)}{2}-\sum_{m=2}^{x}\Phi\!\left(\left\lfloor\frac{x}{m}\right\rfloor\right).$$

This is exactly the large-argument formula used by the implementation.

Step 4: Compress the recurrence with floor-division blocks

The quantity \(\left\lfloor x/m\right\rfloor\) does not change at every index. If

$$q=\left\lfloor\frac{x}{\ell}\right\rfloor,\qquad r=\left\lfloor\frac{x}{q}\right\rfloor,$$

then every \(m\in[\ell,r]\) has the same quotient \(q\). Therefore the recurrence becomes

$$\Phi(x)=\frac{x(x+1)}{2}-\sum_{\text{blocks }[\ell,r]} (r-\ell+1)\,\Phi(q).$$

Instead of iterating over all \(m\), we only iterate over the distinct quotient blocks. For large \(x\), there are only \(O(\sqrt{x})\) such blocks.

The C++, Python, and Java implementations precompute \(\varphi(n)\) and its prefix values for all \(n\le 5\times 10^6\), and only use this recursive block formula beyond that cutoff.

Step 5: Apply the same block idea to the outer sum

The transformed target expression

$$G(N)=\sum_{d=1}^{N} d\,\Phi\!\left(\left\lfloor\frac{N}{d}\right\rfloor\right)$$

has the same floor structure. If \(\left\lfloor N/d\right\rfloor=q\) on a block \(d\in[\ell,r]\), then the entire block contributes

$$\left(\sum_{d=\ell}^{r} d\right)\Phi(q).$$

The arithmetic progression sum is

$$\sum_{d=\ell}^{r} d=\frac{(\ell+r)(r-\ell+1)}{2}.$$

So the final evaluation is again a loop over quotient blocks rather than a loop over all \(d\le N\).

Worked Example: \(N=10\)

The implementations verify the small checkpoint \(G(10)=122\). The block formula reproduces it directly.

First compute the needed totient prefixes:

$$\Phi(1)=1,\qquad \Phi(2)=2,\qquad \Phi(3)=4,\qquad \Phi(5)=10,\qquad \Phi(10)=32.$$

Now write

$$G(10)=\sum_{d=1}^{10} d\,\Phi\!\left(\left\lfloor\frac{10}{d}\right\rfloor\right).$$

The quotients \(\left\lfloor 10/d\right\rfloor\) are constant on the blocks

$$[1,1],\qquad [2,2],\qquad [3,3],\qquad [4,5],\qquad [6,10],$$

with quotient values \(10,5,3,2,1\), respectively. Hence

$$\begin{aligned} G(10)&=1\cdot \Phi(10)+2\cdot \Phi(5)+3\cdot \Phi(3)+(4+5)\Phi(2)+(6+7+8+9+10)\Phi(1)\\ &=1\cdot 32+2\cdot 10+3\cdot 4+9\cdot 2+40\cdot 1\\ &=32+20+12+18+40=122. \end{aligned}$$

This small case shows exactly why grouping by equal floor quotients removes the need for a linear scan.

How the Code Works

The C++, Python, and Java implementations all follow the same pipeline. First they build \(\varphi(n)\) up to \(5\times 10^6\) with a linear sieve and turn it into a prefix table, so every small \(\Phi(x)\) query becomes a constant-time lookup.

For larger arguments, the implementation evaluates \(\Phi(x)\) from the recurrence above, grouping equal values of \(\left\lfloor x/m\right\rfloor\) into blocks and caching each large result the first time it is computed. The same large quotient can appear repeatedly, so memoization is essential.

After that, the implementation computes \(G(N)\) with another quotient-block loop over \(d\). On each block it evaluates the arithmetic progression sum, multiplies it by the already known value of \(\Phi(q)\), and accumulates everything modulo \(998244353\). Because the modulus is odd, every division by \(2\) is handled as multiplication by the modular inverse of \(2\).

Complexity Analysis

Let \(L=5\times 10^6\). The totient sieve and prefix construction cost \(O(L)\) time and \(O(L)\) memory. Each large \(\Phi(x)\) query is memoized once and processed by quotient blocks, so its work is proportional to the number of distinct values of \(\left\lfloor x/k\right\rfloor\), namely \(O(\sqrt{x})\). The outer sum is also traversed by quotient blocks, so the full computation is far below \(O(N)\) and is practical for \(N=10^{11}\).

Footnotes and References

  1. Problem page: https://projecteuler.net/problem=625
  2. Greatest common divisor: Wikipedia - Greatest common divisor
  3. Euler's totient function: Wikipedia - Euler's totient function
  4. Dirichlet convolution: Wikipedia - Dirichlet convolution
  5. Linear sieve: cp-algorithms - Linear Sieve

Problemzusammenfassung

Definiert ist

$$G(N)=\sum_{j=1}^{N}\sum_{i=1}^{j}\gcd(i,j).$$

Gesucht ist \(G(10^{11}) \bmod 998244353\). Eine direkte Doppelschleife ist dafür viel zu langsam, also formt die Lösung die ggT-Summe in eine Teilersumme um, führt die summatorische Phi-Funktion ein und wertet die verbleibenden Summen mit Quotientenblöcken aus.

Mathematischer Ansatz

Wir isolieren zunächst die innere Summe

$$T(j)=\sum_{i=1}^{j}\gcd(i,j),$$

und definieren außerdem die Phi-Präfixsumme

$$\Phi(x)=\sum_{k=1}^{x}\varphi(k).$$

Die Methode besteht darin, \(T(j)\) so umzuformen, dass nur noch Teiler von \(j\) auftreten, und anschließend auszunutzen, dass ganzzahlige Quotienten auf langen Intervallen konstant bleiben.

Schritt 1: Die innere ggT-Summe als Teilersumme umschreiben

Fixiere \(j\). Wir gruppieren alle \(i\) nach \(d=\gcd(i,j)\). Für festes \(d\) gilt dann

$$i=d\,a,\qquad j=d\,b,\qquad \gcd(a,b)=1.$$

Dabei ist \(b=j/d\), und die zulässigen Werte von \(a\) sind genau die Zahlen \(1\le a\le b\), die zu \(b\) teilerfremd sind. Davon gibt es \(\varphi(b)\).

Somit trägt jeder Teiler \(d\mid j\) den Wert \(d\) genau \(\varphi(j/d)\)-mal bei, also

$$T(j)=\sum_{d\mid j} d\,\varphi\!\left(\frac{j}{d}\right).$$

Damit verschwindet der explizite ggT aus der inneren Schleife und wird durch eine saubere Formel über Teiler ersetzt.

Schritt 2: Die Summationsreihenfolge vertauschen

Nun summieren wir \(T(j)\) über alle \(j\le N\):

$$G(N)=\sum_{j=1}^{N}\sum_{d\mid j} d\,\varphi\!\left(\frac{j}{d}\right).$$

Schreibe \(j=d\,k\). Jedes Paar \((d,k)\) mit \(d\,k\le N\) tritt genau einmal auf, also

$$G(N)=\sum_{d=1}^{N} d\sum_{k\le N/d}\varphi(k)=\sum_{d=1}^{N} d\,\Phi\!\left(\left\lfloor\frac{N}{d}\right\rfloor\right).$$

Damit reduziert sich das gesamte Problem auf viele Abfragen der Präfixsumme \(\Phi(x)\).

Schritt 3: Eine Rekursion für die Phi-Präfixsumme herleiten

Die klassische Identität

$$\sum_{d\mid n}\varphi(d)=n$$

gilt für jedes positive \(n\). Summiert man sie für \(1\le n\le x\), erhält man

$$\sum_{n=1}^{x} n=\sum_{n=1}^{x}\sum_{d\mid n}\varphi(d).$$

Nach Vertauschen der Summen sieht man: Ein fixes \(d\) teilt genau \(\left\lfloor x/d\right\rfloor\) Zahlen bis \(x\). Also

$$\frac{x(x+1)}{2}=\sum_{d=1}^{x}\varphi(d)\left\lfloor\frac{x}{d}\right\rfloor.$$

Gruppiert man die rechte Seite nach gleichen Quotienten, so folgt

$$\frac{x(x+1)}{2}=\sum_{m=1}^{x}\Phi\!\left(\left\lfloor\frac{x}{m}\right\rfloor\right).$$

Durch Isolieren des Terms \(m=1\) erhält man die Rekursion

$$\Phi(x)=\frac{x(x+1)}{2}-\sum_{m=2}^{x}\Phi\!\left(\left\lfloor\frac{x}{m}\right\rfloor\right).$$

Genau diese Formel verwenden die Implementierungen für große Argumente.

Schritt 4: Die Rekursion mit Quotientenblöcken komprimieren

Der Ausdruck \(\left\lfloor x/m\right\rfloor\) ändert sich nicht bei jedem Index. Setzt man

$$q=\left\lfloor\frac{x}{\ell}\right\rfloor,\qquad r=\left\lfloor\frac{x}{q}\right\rfloor,$$

dann gilt für jedes \(m\in[\ell,r]\) derselbe Quotient \(q\). Daher wird die Rekursion zu

$$\Phi(x)=\frac{x(x+1)}{2}-\sum_{\text{Blöcke }[\ell,r]} (r-\ell+1)\,\Phi(q).$$

Statt alle \(m\) einzeln zu durchlaufen, betrachtet man nur die verschiedenen Quotientenblöcke. Für große \(x\) gibt es davon nur \(O(\sqrt{x})\).

Die C++-, Python- und Java-Implementierungen berechnen \(\varphi(n)\) und die Präfixsummen für alle \(n\le 5\times 10^6\) vor und greifen erst oberhalb dieser Grenze auf die rekursive Blockformel zurück.

Schritt 5: Dieselbe Blocktechnik auf die äußere Summe anwenden

Auch die transformierte Zielformel

$$G(N)=\sum_{d=1}^{N} d\,\Phi\!\left(\left\lfloor\frac{N}{d}\right\rfloor\right)$$

hat dieselbe Struktur. Wenn \(\left\lfloor N/d\right\rfloor=q\) auf einem Block \(d\in[\ell,r]\) konstant ist, dann ist der gesamte Blockbeitrag

$$\left(\sum_{d=\ell}^{r} d\right)\Phi(q).$$

Die arithmetische Reihe summiert sich zu

$$\sum_{d=\ell}^{r} d=\frac{(\ell+r)(r-\ell+1)}{2}.$$

Damit besteht die Endauswertung erneut nur aus einer Schleife über Quotientenblöcke statt über alle \(d\le N\).

Durchgerechnetes Beispiel: \(N=10\)

Die Implementierungen prüfen den kleinen Kontrollwert \(G(10)=122\). Mit der Blockformel lässt sich das direkt nachvollziehen.

Benötigt werden zunächst nur diese Phi-Präfixwerte:

$$\Phi(1)=1,\qquad \Phi(2)=2,\qquad \Phi(3)=4,\qquad \Phi(5)=10,\qquad \Phi(10)=32.$$

Dann gilt

$$G(10)=\sum_{d=1}^{10} d\,\Phi\!\left(\left\lfloor\frac{10}{d}\right\rfloor\right).$$

Die Quotienten \(\left\lfloor 10/d\right\rfloor\) sind auf den Blöcken

$$[1,1],\qquad [2,2],\qquad [3,3],\qquad [4,5],\qquad [6,10]$$

konstant und haben dort die Werte \(10,5,3,2,1\). Also

$$\begin{aligned} G(10)&=1\cdot \Phi(10)+2\cdot \Phi(5)+3\cdot \Phi(3)+(4+5)\Phi(2)+(6+7+8+9+10)\Phi(1)\\ &=1\cdot 32+2\cdot 10+3\cdot 4+9\cdot 2+40\cdot 1\\ &=32+20+12+18+40=122. \end{aligned}$$

Dieses kleine Beispiel zeigt genau, warum die Gruppierung nach gleichen Quotienten den linearen Durchlauf ersetzt.

Wie der Code arbeitet

Die C++-, Python- und Java-Implementierungen folgen derselben Pipeline. Zuerst bauen sie mit einem linearen Sieb die Werte \(\varphi(n)\) bis \(5\times 10^6\) auf und erzeugen daraus eine Präfixtabelle, sodass jede kleine Abfrage \(\Phi(x)\) in \(O(1)\) beantwortet wird.

Für größere Argumente berechnet die Implementierung \(\Phi(x)\) aus der obigen Rekursion, fasst gleiche Werte von \(\left\lfloor x/m\right\rfloor\) zu Blöcken zusammen und speichert jedes große Ergebnis beim ersten Auftreten. Diese Memoisierung ist entscheidend, weil dieselben großen Quotienten in der Rekursion und in der Endsumme mehrfach auftauchen.

Anschließend wird \(G(N)\) über eine zweite Quotientenblock-Schleife berechnet. Pro Block wird die Summe der arithmetischen Reihe ausgewertet, mit dem bereits bekannten Wert \(\Phi(q)\) multipliziert und alles modulo \(998244353\) aufsummiert. Da der Modulus ungerade ist, werden die Divisionen durch \(2\) als Multiplikation mit dem modularen Inversen von \(2\) umgesetzt.

Komplexitätsanalyse

Mit \(L=5\times 10^6\) kostet das Sieb samt Präfixtabelle \(O(L)\) Zeit und \(O(L)\) Speicher. Jede große Abfrage \(\Phi(x)\) wird genau einmal memoisiert und über Quotientenblöcke ausgewertet, also proportional zur Anzahl verschiedener Werte von \(\left\lfloor x/k\right\rfloor\), d. h. in \(O(\sqrt{x})\). Auch die äußere Summe läuft nur über solche Blöcke. Insgesamt liegt die Laufzeit damit weit unter \(O(N)\) und ist für \(N=10^{11}\) praktikabel.

Fußnoten und Quellen

  1. Problemseite: https://projecteuler.net/problem=625
  2. Größter gemeinsamer Teiler: Wikipedia - Greatest common divisor
  3. Eulersche Phi-Funktion: Wikipedia - Euler's totient function
  4. Dirichlet-Faltung: Wikipedia - Dirichlet convolution
  5. Lineares Sieb: cp-algorithms - Linear Sieve

Problem Özeti

Tanım

$$G(N)=\sum_{j=1}^{N}\sum_{i=1}^{j}\gcd(i,j).$$

İstenen değer \(G(10^{11}) \bmod 998244353\)'dir. Bu boyutta doğrudan çift döngü kurmak mümkün olmadığı için çözüm, gcd toplamını bölenler toplamına dönüştürür, Euler phi önek toplamını kullanır ve kalan ifadeleri tam bölüm bloklarıyla hesaplar.

Matematiksel Yaklaşım

Önce iç toplamı

$$T(j)=\sum_{i=1}^{j}\gcd(i,j)$$

şeklinde ayıralım. Ayrıca

$$\Phi(x)=\sum_{k=1}^{x}\varphi(k)$$

tanımıyla phi önek toplamını kullanalım. Yöntemin özü, \(T(j)\)'yi yalnızca \(j\)'nin bölenlerine bağlı bir forma çevirmek ve ardından \(\lfloor x/k\rfloor\) ifadesinin uzun aralıklarda sabit kalmasını kullanmaktır.

Adım 1: İç gcd toplamını bölenler toplamı olarak yaz

Sabit bir \(j\) için indisleri \(d=\gcd(i,j)\) değerine göre gruplayalım. \(d\) sabitken

$$i=d\,a,\qquad j=d\,b,\qquad \gcd(a,b)=1$$

yazılabilir. Burada \(b=j/d\)'dir ve geçerli \(a\) değerleri, \(1\le a\le b\) aralığında olup \(b\) ile aralarında asal olan sayılardır. Bunların sayısı \(\varphi(b)\) eder.

Dolayısıyla her \(d\mid j\) böleni, \(d\) değerini tam olarak \(\varphi(j/d)\) kez katkı yapar:

$$T(j)=\sum_{d\mid j} d\,\varphi\!\left(\frac{j}{d}\right).$$

Böylece iç döngüdeki açık gcd hesabı ortadan kalkar ve yerini temiz bir bölen formülü alır.

Adım 2: Toplam sırasını değiştir

Şimdi \(T(j)\)'yi tüm \(j\le N\) için toplayalım:

$$G(N)=\sum_{j=1}^{N}\sum_{d\mid j} d\,\varphi\!\left(\frac{j}{d}\right).$$

\(j=d\,k\) yazarsak, \(d\,k\le N\) koşulunu sağlayan her \((d,k)\) çifti tam bir kez görünür. Bu da

$$G(N)=\sum_{d=1}^{N} d\sum_{k\le N/d}\varphi(k)=\sum_{d=1}^{N} d\,\Phi\!\left(\left\lfloor\frac{N}{d}\right\rfloor\right)$$

sonucunu verir. Artık problem, çok sayıda \(\Phi(x)\) sorgusunu hızlı cevaplamaya indirgenmiştir.

Adım 3: Phi önek toplamı için bir bağıntı çıkar

Klasik özdeşlik

$$\sum_{d\mid n}\varphi(d)=n$$

her pozitif \(n\) için doğrudur. Bunu \(1\le n\le x\) aralığında toplarsak

$$\sum_{n=1}^{x} n=\sum_{n=1}^{x}\sum_{d\mid n}\varphi(d)$$

elde edilir. Toplamların yerini değiştirince, sabit bir \(d\) sayısının \(x\)'e kadar tam olarak \(\left\lfloor x/d\right\rfloor\) sayıyı böldüğünü görürüz. Böylece

$$\frac{x(x+1)}{2}=\sum_{d=1}^{x}\varphi(d)\left\lfloor\frac{x}{d}\right\rfloor$$

olur. Sağ tarafı eşit bölüm değerlerine göre gruplayınca

$$\frac{x(x+1)}{2}=\sum_{m=1}^{x}\Phi\!\left(\left\lfloor\frac{x}{m}\right\rfloor\right)$$

yazabiliriz. \(m=1\) terimini ayırınca şu bağıntı çıkar:

$$\Phi(x)=\frac{x(x+1)}{2}-\sum_{m=2}^{x}\Phi\!\left(\left\lfloor\frac{x}{m}\right\rfloor\right).$$

Uygulama büyük argümanlar için tam olarak bu formülü kullanır.

Adım 4: Bağıntıyı bölüm bloklarıyla sıkıştır

\(\left\lfloor x/m\right\rfloor\) ifadesi her adımda değişmez. Eğer

$$q=\left\lfloor\frac{x}{\ell}\right\rfloor,\qquad r=\left\lfloor\frac{x}{q}\right\rfloor$$

ise, \(m\in[\ell,r]\) aralığındaki tüm değerler aynı \(q\) bölümünü verir. Bu nedenle bağıntı

$$\Phi(x)=\frac{x(x+1)}{2}-\sum_{\text{bloklar }[\ell,r]} (r-\ell+1)\,\Phi(q)$$

şekline gelir. Yani tüm \(m\) değerlerini tek tek dolaşmak yerine yalnızca farklı bölüm bloklarını gezeriz. Büyük \(x\) için bu blokların sayısı \(O(\sqrt{x})\) mertebesindedir.

C++, Python ve Java implementasyonları \(\varphi(n)\) ile önek toplamını \(n\le 5\times 10^6\) için önceden hesaplar; bu sınırın üstünde ise rekürsif blok formülüne geçer.

Adım 5: Aynı blok fikrini dış toplama uygula

Dönüştürülmüş hedef ifade

$$G(N)=\sum_{d=1}^{N} d\,\Phi\!\left(\left\lfloor\frac{N}{d}\right\rfloor\right)$$

aynı tam bölüm yapısına sahiptir. Eğer \(\left\lfloor N/d\right\rfloor=q\) değeri \(d\in[\ell,r]\) aralığında sabitse, bu bloğun katkısı

$$\left(\sum_{d=\ell}^{r} d\right)\Phi(q)$$

olur. Aritmetik dizi toplamı da

$$\sum_{d=\ell}^{r} d=\frac{(\ell+r)(r-\ell+1)}{2}$$

şeklindedir. Böylece son hesap yine doğrusal tarama yerine bloklar üzerinde yapılır.

Çözümlü Örnek: \(N=10\)

Implementasyonlar küçük kontrol değeri olarak \(G(10)=122\) sonucunu doğrular. Bunu formülden doğrudan görelim.

Gerekli phi önekleri şunlardır:

$$\Phi(1)=1,\qquad \Phi(2)=2,\qquad \Phi(3)=4,\qquad \Phi(5)=10,\qquad \Phi(10)=32.$$

Sonra

$$G(10)=\sum_{d=1}^{10} d\,\Phi\!\left(\left\lfloor\frac{10}{d}\right\rfloor\right)$$

yazarız. \(\left\lfloor 10/d\right\rfloor\) değerleri şu bloklarda sabittir:

$$[1,1],\qquad [2,2],\qquad [3,3],\qquad [4,5],\qquad [6,10].$$

Bu bloklardaki bölüm değerleri sırasıyla \(10,5,3,2,1\) olduğundan

$$\begin{aligned} G(10)&=1\cdot \Phi(10)+2\cdot \Phi(5)+3\cdot \Phi(3)+(4+5)\Phi(2)+(6+7+8+9+10)\Phi(1)\\ &=1\cdot 32+2\cdot 10+3\cdot 4+9\cdot 2+40\cdot 1\\ &=32+20+12+18+40=122. \end{aligned}$$

Bu küçük örnek, eşit bölüm değerlerine göre gruplamanın neden işe yaradığını açık biçimde gösterir.

Kodun Çalışma Mantığı

C++, Python ve Java implementasyonlarının akışı aynıdır. Önce lineer elekle \(\varphi(n)\) değerleri \(5\times 10^6\)'ya kadar hesaplanır ve bunlardan bir önek tablosu oluşturulur; böylece küçük \(\Phi(x)\) sorguları sabit zamanda cevaplanır.

Daha büyük argümanlar için implementasyon, yukarıdaki bağıntıyı kullanarak \(\Phi(x)\) hesabını tam bölüm bloklarına ayırır ve bulunan her büyük sonucu ilk hesaplamadan sonra önbelleğe alır. Aynı büyük bölüm değerleri rekürsiyon içinde ve son toplamda tekrar tekrar ortaya çıktığı için bu önbellekleme kritik önemdedir.

Son aşamada \(G(N)\), \(d\) üzerinde ikinci bir blok döngüsüyle hesaplanır. Her blokta aritmetik dizi toplamı bulunur, ilgili \(\Phi(q)\) ile çarpılır ve tüm katkılar \(998244353\) modunda toplanır. Modül tek sayı olduğu için \(2\)'ye bölmeler, \(2\)'nin modüler tersiyle çarpma olarak uygulanır.

Karmaşıklık Analizi

\(L=5\times 10^6\) olsun. Elek ve önek tablo kurma maliyeti \(O(L)\) zaman ve \(O(L)\) bellektir. Her büyük \(\Phi(x)\) sorgusu yalnızca bir kez memoize edilir ve farklı \(\left\lfloor x/k\right\rfloor\) değerleri kadar blokla işlenir; bu sayı \(O(\sqrt{x})\) mertebesindedir. Dış toplam da aynı şekilde bloklar üzerinden yürütülür. Sonuç olarak toplam maliyet \(O(N)\)'den çok daha küçüktür ve \(N=10^{11}\) için pratiktir.

Dipnotlar ve Kaynakça

  1. Problem sayfası: https://projecteuler.net/problem=625
  2. En büyük ortak bölen: Wikipedia - Greatest common divisor
  3. Euler phi fonksiyonu: Wikipedia - Euler's totient function
  4. Dirichlet konvolüsyonu: Wikipedia - Dirichlet convolution
  5. Lineer elek: cp-algorithms - Linear Sieve

Resumen del Problema

Se define

$$G(N)=\sum_{j=1}^{N}\sum_{i=1}^{j}\gcd(i,j).$$

Hay que calcular \(G(10^{11}) \bmod 998244353\). Un doble bucle directo es inviable a esa escala, así que la solución convierte la suma de gcd en una suma sobre divisores, introduce el prefijo acumulado de \(\varphi\) y evalúa las expresiones restantes mediante bloques de cociente constante.

Enfoque Matemático

Primero aislamos la suma interior

$$T(j)=\sum_{i=1}^{j}\gcd(i,j),$$

y también definimos el prefijo de la función phi

$$\Phi(x)=\sum_{k=1}^{x}\varphi(k).$$

La estrategia consiste en reescribir \(T(j)\) para que solo dependa de los divisores de \(j\), y luego aprovechar que los cocientes enteros permanecen constantes en intervalos largos.

Paso 1: Reescribir la suma interior de mcd como suma sobre divisores

Fijemos \(j\). Agrupamos los índices \(i\) por el valor \(d=\gcd(i,j)\). Si \(d\) está fijado, podemos escribir

$$i=d\,a,\qquad j=d\,b,\qquad \gcd(a,b)=1.$$

Aquí \(b=j/d\), y los valores posibles de \(a\) son exactamente los enteros \(1\le a\le b\) coprimos con \(b\). Su cantidad es \(\varphi(b)\).

Por tanto, cada divisor \(d\mid j\) aporta el valor \(d\) exactamente \(\varphi(j/d)\) veces, de modo que

$$T(j)=\sum_{d\mid j} d\,\varphi\!\left(\frac{j}{d}\right).$$

Esta identidad elimina el gcd explícito del bucle interior y lo sustituye por una fórmula limpia sobre divisores.

Paso 2: Invertir el orden de sumación

Ahora sumamos \(T(j)\) sobre todos los \(j\le N\):

$$G(N)=\sum_{j=1}^{N}\sum_{d\mid j} d\,\varphi\!\left(\frac{j}{d}\right).$$

Si escribimos \(j=d\,k\), cada pareja \((d,k)\) con \(d\,k\le N\) aparece una sola vez. Por eso

$$G(N)=\sum_{d=1}^{N} d\sum_{k\le N/d}\varphi(k)=\sum_{d=1}^{N} d\,\Phi\!\left(\left\lfloor\frac{N}{d}\right\rfloor\right).$$

Todo el problema queda reducido a responder muchas consultas del prefijo \(\Phi(x)\).

Paso 3: Derivar una recurrencia para el prefijo de la función phi

La identidad clásica

$$\sum_{d\mid n}\varphi(d)=n$$

vale para todo entero positivo \(n\). Si la sumamos para \(1\le n\le x\), obtenemos

$$\sum_{n=1}^{x} n=\sum_{n=1}^{x}\sum_{d\mid n}\varphi(d).$$

Al intercambiar el orden de sumación, un divisor fijo \(d\) aparece exactamente en \(\left\lfloor x/d\right\rfloor\) múltiplos hasta \(x\). Entonces

$$\frac{x(x+1)}{2}=\sum_{d=1}^{x}\varphi(d)\left\lfloor\frac{x}{d}\right\rfloor.$$

Si agrupamos por cocientes iguales, la expresión se transforma en

$$\frac{x(x+1)}{2}=\sum_{m=1}^{x}\Phi\!\left(\left\lfloor\frac{x}{m}\right\rfloor\right).$$

Aislando el término \(m=1\), se obtiene la recurrencia

$$\Phi(x)=\frac{x(x+1)}{2}-\sum_{m=2}^{x}\Phi\!\left(\left\lfloor\frac{x}{m}\right\rfloor\right).$$

Esa es exactamente la fórmula que usan las implementaciones para argumentos grandes.

Paso 4: Comprimir la recurrencia con bloques de cociente

La cantidad \(\left\lfloor x/m\right\rfloor\) no cambia en cada índice. Si

$$q=\left\lfloor\frac{x}{\ell}\right\rfloor,\qquad r=\left\lfloor\frac{x}{q}\right\rfloor,$$

entonces todo \(m\in[\ell,r]\) comparte el mismo cociente \(q\). Por ello la recurrencia pasa a ser

$$\Phi(x)=\frac{x(x+1)}{2}-\sum_{\text{bloques }[\ell,r]} (r-\ell+1)\,\Phi(q).$$

En vez de recorrer todos los \(m\), solo recorremos los bloques de cociente distinto. Para \(x\) grande, el número de esos bloques es \(O(\sqrt{x})\).

Las implementaciones en C++, Python y Java precalculan \(\varphi(n)\) y su prefijo para todos los \(n\le 5\times 10^6\), y solo usan esta fórmula recursiva por bloques por encima de ese umbral.

Paso 5: Aplicar la misma idea de bloques a la suma exterior

La expresión transformada

$$G(N)=\sum_{d=1}^{N} d\,\Phi\!\left(\left\lfloor\frac{N}{d}\right\rfloor\right)$$

tiene la misma estructura. Si \(\left\lfloor N/d\right\rfloor=q\) es constante en un bloque \(d\in[\ell,r]\), la contribución del bloque completo es

$$\left(\sum_{d=\ell}^{r} d\right)\Phi(q).$$

La suma de la progresión aritmética es

$$\sum_{d=\ell}^{r} d=\frac{(\ell+r)(r-\ell+1)}{2}.$$

Así, la evaluación final vuelve a ser un recorrido por bloques de cociente y no por todos los \(d\le N\).

Ejemplo trabajado: \(N=10\)

Las implementaciones verifican el punto de control pequeño \(G(10)=122\). La fórmula por bloques lo reproduce de inmediato.

Primero necesitamos estos prefijos de \(\varphi\):

$$\Phi(1)=1,\qquad \Phi(2)=2,\qquad \Phi(3)=4,\qquad \Phi(5)=10,\qquad \Phi(10)=32.$$

Luego escribimos

$$G(10)=\sum_{d=1}^{10} d\,\Phi\!\left(\left\lfloor\frac{10}{d}\right\rfloor\right).$$

Los cocientes \(\left\lfloor 10/d\right\rfloor\) son constantes en los bloques

$$[1,1],\qquad [2,2],\qquad [3,3],\qquad [4,5],\qquad [6,10],$$

con valores \(10,5,3,2,1\), respectivamente. Por lo tanto,

$$\begin{aligned} G(10)&=1\cdot \Phi(10)+2\cdot \Phi(5)+3\cdot \Phi(3)+(4+5)\Phi(2)+(6+7+8+9+10)\Phi(1)\\ &=1\cdot 32+2\cdot 10+3\cdot 4+9\cdot 2+40\cdot 1\\ &=32+20+12+18+40=122. \end{aligned}$$

Este caso pequeño muestra con claridad por qué agrupar por cocientes iguales elimina la necesidad de un barrido lineal.

Cómo funciona el código

Las implementaciones en C++, Python y Java siguen la misma secuencia. Primero construyen \(\varphi(n)\) hasta \(5\times 10^6\) con una criba lineal y lo convierten en una tabla de prefijos, de modo que cualquier consulta pequeña de \(\Phi(x)\) se resuelve en tiempo constante.

Para argumentos grandes, la implementación calcula \(\Phi(x)\) con la recurrencia anterior, agrupa los valores iguales de \(\left\lfloor x/m\right\rfloor\) en bloques y guarda cada resultado grande la primera vez que aparece. Esa memoización es esencial, porque los mismos cocientes grandes se repiten muchas veces.

Después, la implementación evalúa \(G(N)\) con otro recorrido por bloques sobre \(d\). En cada bloque calcula la suma de la progresión aritmética, la multiplica por el valor ya conocido de \(\Phi(q)\) y acumula todo módulo \(998244353\). Como el módulo es impar, cada división entre \(2\) se realiza multiplicando por el inverso modular de \(2\).

Análisis de Complejidad

Sea \(L=5\times 10^6\). La criba y la construcción del prefijo cuestan \(O(L)\) tiempo y \(O(L)\) memoria. Cada consulta grande \(\Phi(x)\) se memoiza una sola vez y se procesa por bloques de cociente, así que su coste es proporcional al número de valores distintos de \(\left\lfloor x/k\right\rfloor\), es decir, \(O(\sqrt{x})\). La suma exterior también se recorre por bloques, de modo que el coste total queda muy por debajo de \(O(N)\) y resulta práctico para \(N=10^{11}\).

Notas y Referencias

  1. Página del problema: https://projecteuler.net/problem=625
  2. Máximo común divisor: Wikipedia - Greatest common divisor
  3. Función phi de Euler: Wikipedia - Euler's totient function
  4. Convolución de Dirichlet: Wikipedia - Dirichlet convolution
  5. Criba lineal: cp-algorithms - Linear Sieve

问题概述

定义

$$G(N)=\sum_{j=1}^{N}\sum_{i=1}^{j}\gcd(i,j).$$

题目要求计算 \(G(10^{11}) \bmod 998244353\)。如果直接做双重循环,规模完全不可接受,因此解法必须先把 gcd 求和改写成约数求和,再把问题化成欧拉函数前缀和,最后利用整除分块处理大量相同的商。

数学方法

先把内层和单独记成

$$T(j)=\sum_{i=1}^{j}\gcd(i,j),$$

再定义欧拉函数的前缀和

$$\Phi(x)=\sum_{k=1}^{x}\varphi(k).$$

整个方法的核心是两次变形。第一步把 \(T(j)\) 写成只和 \(j\) 的约数有关的形式;第二步利用 \(\left\lfloor x/k\right\rfloor\) 在长区间上保持不变这一事实,把看似很长的求和压缩成少量分块。

步骤 1:把内层 gcd 求和改写成约数和

固定 \(j\)。把所有 \(i\) 按照 \(d=\gcd(i,j)\) 分类。如果某一类的 gcd 恰好等于 \(d\),那么可以写成

$$i=d\,a,\qquad j=d\,b,\qquad \gcd(a,b)=1.$$

这里 \(b=j/d\)。对固定的 \(d\) 来说,允许的 \(a\) 正好是区间 \(1\le a\le b\) 中与 \(b\) 互素的整数个数,而这个数量就是 \(\varphi(b)\)。

因此,每个约数 \(d\mid j\) 会把数值 \(d\) 贡献 \(\varphi(j/d)\) 次,从而得到

$$T(j)=\sum_{d\mid j} d\,\varphi\!\left(\frac{j}{d}\right).$$

这一步非常关键,因为它把原本需要显式计算 gcd 的内层循环,转成了标准的约数型表达式。

步骤 2:交换求和顺序

现在把 \(T(j)\) 对所有 \(j\le N\) 累加:

$$G(N)=\sum_{j=1}^{N}\sum_{d\mid j} d\,\varphi\!\left(\frac{j}{d}\right).$$

令 \(j=d\,k\)。只要满足 \(d\,k\le N\),这样的 \((d,k)\) 配对就会且只会出现一次,所以

$$G(N)=\sum_{d=1}^{N} d\sum_{k\le N/d}\varphi(k)=\sum_{d=1}^{N} d\,\Phi\!\left(\left\lfloor\frac{N}{d}\right\rfloor\right).$$

这样一来,原问题就变成了大量求 \(\Phi(x)\) 的问题。只要 \(\Phi(x)\) 能快速得到,整个题就可以在可接受的时间内完成。

步骤 3:推导欧拉函数前缀和的递推式

我们使用经典恒等式

$$\sum_{d\mid n}\varphi(d)=n,$$

它对每个正整数 \(n\) 都成立。把它从 \(1\) 累加到 \(x\),有

$$\sum_{n=1}^{x} n=\sum_{n=1}^{x}\sum_{d\mid n}\varphi(d).$$

交换求和顺序后可以看到:固定一个 \(d\),在 \(1\) 到 \(x\) 中被它整除的数恰好有 \(\left\lfloor x/d\right\rfloor\) 个,因此

$$\frac{x(x+1)}{2}=\sum_{d=1}^{x}\varphi(d)\left\lfloor\frac{x}{d}\right\rfloor.$$

再把右边按相同的商重新分组,可以写成

$$\frac{x(x+1)}{2}=\sum_{m=1}^{x}\Phi\!\left(\left\lfloor\frac{x}{m}\right\rfloor\right).$$

把 \(m=1\) 的那一项移到左边,就得到递推式

$$\Phi(x)=\frac{x(x+1)}{2}-\sum_{m=2}^{x}\Phi\!\left(\left\lfloor\frac{x}{m}\right\rfloor\right).$$

这正是实现中处理大参数时所使用的公式。

步骤 4:用整除分块压缩递推

虽然上面的递推仍然看起来要从 \(2\) 加到 \(x\),但 \(\left\lfloor x/m\right\rfloor\) 并不会每一步都变化。若设

$$q=\left\lfloor\frac{x}{\ell}\right\rfloor,\qquad r=\left\lfloor\frac{x}{q}\right\rfloor,$$

那么区间 \(m\in[\ell,r]\) 内的每个位置都对应同一个商 \(q\)。于是递推可以压缩成

$$\Phi(x)=\frac{x(x+1)}{2}-\sum_{\text{分块 }[\ell,r]} (r-\ell+1)\,\Phi(q).$$

换句话说,我们不必枚举每一个 \(m\),只需要枚举不同的商对应的分块即可。对大的 \(x\) 来说,这样的分块只有 \(O(\sqrt{x})\) 个。

C++、Python 和 Java 实现都先用线性筛预处理 \(n\le 5\times 10^6\) 的 \(\varphi(n)\) 及其前缀和;只有超过这个阈值的 \(\Phi(x)\) 才用上面的递推加记忆化来求。

步骤 5:把同样的分块思想用于外层求和

变形后的目标式

$$G(N)=\sum_{d=1}^{N} d\,\Phi\!\left(\left\lfloor\frac{N}{d}\right\rfloor\right)$$

也有完全相同的整除结构。如果在某个区间 \(d\in[\ell,r]\) 上,\(\left\lfloor N/d\right\rfloor=q\) 保持不变,那么整个区间的贡献就是

$$\left(\sum_{d=\ell}^{r} d\right)\Phi(q).$$

其中等差和可直接写成

$$\sum_{d=\ell}^{r} d=\frac{(\ell+r)(r-\ell+1)}{2}.$$

因此最终求 \(G(N)\) 时,也无需从 \(1\) 一直枚举到 \(N\),只需遍历这些商不变的区间。

示例:\(N=10\)

实现里验证了一个小检查点:\(G(10)=122\)。下面用上面的公式直接算出它。

先列出需要的欧拉函数前缀和值:

$$\Phi(1)=1,\qquad \Phi(2)=2,\qquad \Phi(3)=4,\qquad \Phi(5)=10,\qquad \Phi(10)=32.$$

于是

$$G(10)=\sum_{d=1}^{10} d\,\Phi\!\left(\left\lfloor\frac{10}{d}\right\rfloor\right).$$

\(\left\lfloor 10/d\right\rfloor\) 在以下区间上保持不变:

$$[1,1],\qquad [2,2],\qquad [3,3],\qquad [4,5],\qquad [6,10].$$

这些区间对应的商分别是 \(10,5,3,2,1\),所以

$$\begin{aligned} G(10)&=1\cdot \Phi(10)+2\cdot \Phi(5)+3\cdot \Phi(3)+(4+5)\Phi(2)+(6+7+8+9+10)\Phi(1)\\ &=1\cdot 32+2\cdot 10+3\cdot 4+9\cdot 2+40\cdot 1\\ &=32+20+12+18+40=122. \end{aligned}$$

这个例子说明了整除分块的真正作用:把大量重复的商合并处理,计算量立刻下降。

代码如何实现

C++、Python 和 Java 实现遵循同一条计算流程。首先用线性筛求出 \(5\times 10^6\) 以内的 \(\varphi(n)\),再构造前缀表,这样所有较小的 \(\Phi(x)\) 查询都能在常数时间内返回。

对于更大的参数,代码使用上面的递推式,通过整除分块把相同的 \(\left\lfloor x/m\right\rfloor\) 合并,并把第一次算出的较大 \(\Phi(x)\) 保存起来。由于相同的大商会在递归和最终求和中反复出现,记忆化能显著减少重复工作。

准备好 \(\Phi(x)\) 之后,代码再对外层的 \(d\) 做一次整除分块。每个分块中先算等差和,再乘以对应的 \(\Phi(q)\),最后全程按模 \(998244353\) 累加。因为模数是奇数,涉及除以 \(2\) 的地方都可以改写成乘以 \(2\) 的模逆元。

复杂度分析

记预处理上界为 \(L=5\times 10^6\)。线性筛和前缀表构造需要 \(O(L)\) 时间与 \(O(L)\) 空间。每个较大的 \(\Phi(x)\) 只会记忆化计算一次,而其耗时取决于不同商 \(\left\lfloor x/k\right\rfloor\) 的个数,也就是 \(O(\sqrt{x})\)。外层求和同样只遍历这些分块,因此总时间远小于 \(O(N)\),足以处理 \(N=10^{11}\)。

参考资料

  1. 题目页面:https://projecteuler.net/problem=625
  2. 最大公约数:Wikipedia - Greatest common divisor
  3. 欧拉函数:Wikipedia - Euler's totient function
  4. Dirichlet 卷积:Wikipedia - Dirichlet convolution
  5. 线性筛:cp-algorithms - Linear Sieve

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

Определена сумма

$$G(N)=\sum_{j=1}^{N}\sum_{i=1}^{j}\gcd(i,j).$$

Нужно вычислить \(G(10^{11}) \bmod 998244353\). Прямой двойной цикл здесь невозможен, поэтому решение переписывает сумму gcd через делители, вводит префикс функции Эйлера и затем вычисляет оставшиеся выражения с помощью блочного разбиения по целой части.

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

Сначала выделим внутреннюю сумму

$$T(j)=\sum_{i=1}^{j}\gcd(i,j),$$

а также определим префикс функции Эйлера

$$\Phi(x)=\sum_{k=1}^{x}\varphi(k).$$

Идея состоит в том, чтобы представить \(T(j)\) только через делители числа \(j\), а затем воспользоваться тем, что величина \(\left\lfloor x/k\right\rfloor\) принимает одинаковые значения на длинных интервалах.

Шаг 1: Переписать внутреннюю сумму gcd как сумму по делителям

Зафиксируем \(j\). Разобьем все \(i\) по значению \(d=\gcd(i,j)\). Если это значение равно \(d\), то можно записать

$$i=d\,a,\qquad j=d\,b,\qquad \gcd(a,b)=1.$$

Здесь \(b=j/d\), а допустимые значения \(a\) - это в точности числа \(1\le a\le b\), взаимно простые с \(b\). Их количество равно \(\varphi(b)\).

Значит, каждый делитель \(d\mid j\) вносит значение \(d\) ровно \(\varphi(j/d)\) раз, и поэтому

$$T(j)=\sum_{d\mid j} d\,\varphi\!\left(\frac{j}{d}\right).$$

Так мы убираем явный расчет gcd из внутреннего цикла и заменяем его стандартной формулой по делителям.

Шаг 2: Поменять порядок суммирования

Теперь просуммируем \(T(j)\) по всем \(j\le N\):

$$G(N)=\sum_{j=1}^{N}\sum_{d\mid j} d\,\varphi\!\left(\frac{j}{d}\right).$$

Положим \(j=d\,k\). Каждая пара \((d,k)\) с условием \(d\,k\le N\) появляется ровно один раз, поэтому

$$G(N)=\sum_{d=1}^{N} d\sum_{k\le N/d}\varphi(k)=\sum_{d=1}^{N} d\,\Phi\!\left(\left\lfloor\frac{N}{d}\right\rfloor\right).$$

Тем самым задача сводится к быстрому вычислению большого числа значений \(\Phi(x)\).

Шаг 3: Вывести рекурсию для префикса функции Эйлера

Используем классическое тождество

$$\sum_{d\mid n}\varphi(d)=n,$$

которое верно для любого положительного \(n\). Суммируя его по \(1\le n\le x\), получаем

$$\sum_{n=1}^{x} n=\sum_{n=1}^{x}\sum_{d\mid n}\varphi(d).$$

После перестановки сумм видно, что фиксированный делитель \(d\) делит ровно \(\left\lfloor x/d\right\rfloor\) чисел не больше \(x\). Следовательно,

$$\frac{x(x+1)}{2}=\sum_{d=1}^{x}\varphi(d)\left\lfloor\frac{x}{d}\right\rfloor.$$

Если теперь сгруппировать слагаемые по одинаковому значению частного, получится

$$\frac{x(x+1)}{2}=\sum_{m=1}^{x}\Phi\!\left(\left\lfloor\frac{x}{m}\right\rfloor\right).$$

Выделяя член при \(m=1\), получаем рекурсию

$$\Phi(x)=\frac{x(x+1)}{2}-\sum_{m=2}^{x}\Phi\!\left(\left\lfloor\frac{x}{m}\right\rfloor\right).$$

Именно эта формула используется в реализации для больших аргументов.

Шаг 4: Сжать рекурсию блоками одинакового частного

Выражение \(\left\lfloor x/m\right\rfloor\) меняется не на каждом шаге. Если задать

$$q=\left\lfloor\frac{x}{\ell}\right\rfloor,\qquad r=\left\lfloor\frac{x}{q}\right\rfloor,$$

то для всех \(m\in[\ell,r]\) частное одинаково и равно \(q\). Поэтому рекурсия переписывается так:

$$\Phi(x)=\frac{x(x+1)}{2}-\sum_{\text{блоки }[\ell,r]} (r-\ell+1)\,\Phi(q).$$

Вместо прохода по каждому \(m\) отдельно мы перебираем только различные блоки частного. Для больших \(x\) их число равно \(O(\sqrt{x})\).

Реализации на C++, Python и Java заранее вычисляют \(\varphi(n)\) и префиксные суммы для всех \(n\le 5\times 10^6\), а рекурсивную блочную формулу применяют только выше этого порога.

Шаг 5: Применить ту же блочную идею к внешней сумме

Преобразованная целевая формула

$$G(N)=\sum_{d=1}^{N} d\,\Phi\!\left(\left\lfloor\frac{N}{d}\right\rfloor\right)$$

имеет ту же структуру. Если на блоке \(d\in[\ell,r]\) значение \(\left\lfloor N/d\right\rfloor=q\) постоянно, то вклад всего блока равен

$$\left(\sum_{d=\ell}^{r} d\right)\Phi(q).$$

Сумма арифметической прогрессии здесь равна

$$\sum_{d=\ell}^{r} d=\frac{(\ell+r)(r-\ell+1)}{2}.$$

Значит, и финальное вычисление \(G(N)\) сводится к перебору блоков одинакового частного, а не всех \(d\le N\).

Разобранный пример: \(N=10\)

В реализациях используется маленькая проверка \(G(10)=122\). Покажем, как она получается из формулы.

Нужны следующие значения префикса функции Эйлера:

$$\Phi(1)=1,\qquad \Phi(2)=2,\qquad \Phi(3)=4,\qquad \Phi(5)=10,\qquad \Phi(10)=32.$$

Тогда

$$G(10)=\sum_{d=1}^{10} d\,\Phi\!\left(\left\lfloor\frac{10}{d}\right\rfloor\right).$$

Значения \(\left\lfloor 10/d\right\rfloor\) постоянны на блоках

$$[1,1],\qquad [2,2],\qquad [3,3],\qquad [4,5],\qquad [6,10],$$

где соответствующие частные равны \(10,5,3,2,1\). Следовательно,

$$\begin{aligned} G(10)&=1\cdot \Phi(10)+2\cdot \Phi(5)+3\cdot \Phi(3)+(4+5)\Phi(2)+(6+7+8+9+10)\Phi(1)\\ &=1\cdot 32+2\cdot 10+3\cdot 4+9\cdot 2+40\cdot 1\\ &=32+20+12+18+40=122. \end{aligned}$$

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

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

Реализации на C++, Python и Java работают по одной и той же схеме. Сначала они строят значения \(\varphi(n)\) до \(5\times 10^6\) линейным решетом и превращают их в префиксную таблицу, так что все малые запросы \(\Phi(x)\) становятся обращениями к массиву за \(O(1)\).

Для больших аргументов реализация вычисляет \(\Phi(x)\) по приведенной рекурсии, объединяя одинаковые значения \(\left\lfloor x/m\right\rfloor\) в блоки и сохраняя каждый крупный результат при первом вычислении. Это важно, потому что одни и те же большие частные возникают снова и снова.

После этого \(G(N)\) вычисляется еще одним проходом по блокам по переменной \(d\). На каждом блоке находится сумма арифметической прогрессии, умножается на уже известное значение \(\Phi(q)\), и все аккумулируется по модулю \(998244353\). Так как модуль нечетный, деление на \(2\) реализуется умножением на обратный элемент к \(2\).

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

Пусть \(L=5\times 10^6\). Построение решета и префиксной таблицы требует \(O(L)\) времени и \(O(L)\) памяти. Каждый большой запрос \(\Phi(x)\) мемоизируется только один раз и обрабатывается по блокам одинакового частного, так что его стоимость определяется числом различных значений \(\left\lfloor x/k\right\rfloor\), то есть равна \(O(\sqrt{x})\). Внешняя сумма тоже обходится блоками, поэтому общая сложность значительно меньше \(O(N)\) и практична для \(N=10^{11}\).

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

  1. Страница задачи: https://projecteuler.net/problem=625
  2. Наибольший общий делитель: Wikipedia - Greatest common divisor
  3. Функция Эйлера: Wikipedia - Euler's totient function
  4. Свёртка Дирихле: Wikipedia - Dirichlet convolution
  5. Линейное решето: cp-algorithms - Linear Sieve

ملخص المسألة

نعرّف

$$G(N)=\sum_{j=1}^{N}\sum_{i=1}^{j}\gcd(i,j).$$

المطلوب هو حساب \(G(10^{11}) \bmod 998244353\). تنفيذ حل مباشر بحلقتين متداخلتين غير ممكن عمليًا عند هذا الحجم، لذلك يعيد الحل كتابة مجموع gcd على صورة مجموع على القواسم، ثم يستخدم المجموع التراكمي لدالة أويلر، ثم يضغط الحساب بتقسيمات تعتمد على ثبات القسمة الصحيحة.

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

لنكتب أولًا المجموع الداخلي على الصورة

$$T(j)=\sum_{i=1}^{j}\gcd(i,j),$$

ونعرّف أيضًا المجموع التراكمي لدالة أويلر:

$$\Phi(x)=\sum_{k=1}^{x}\varphi(k).$$

جوهر الفكرة هو تحويل \(T(j)\) إلى تعبير يعتمد فقط على قواسم \(j\)، ثم الاستفادة من أن القيمة \(\left\lfloor x/k\right\rfloor\) تبقى ثابتة على فواصل طويلة، مما يسمح بتجميع كثير من الحدود دفعة واحدة.

الخطوة 1: إعادة كتابة مجموع gcd الداخلي كمجموع على القواسم

ثبّت \(j\). نجمع الحدود بحسب القيمة \(d=\gcd(i,j)\). إذا كانت هذه القيمة ثابتة وتساوي \(d\)، فيمكن كتابة

$$i=d\,a,\qquad j=d\,b,\qquad \gcd(a,b)=1.$$

هنا \(b=j/d\)، أما القيم المسموح بها لـ \(a\) فهي الأعداد من \(1\) إلى \(b\) التي تكون أولية نسبيًا مع \(b\). وعدد هذه القيم يساوي \(\varphi(b)\).

إذن كل قاسم \(d\mid j\) يضيف القيمة \(d\) بعدد مرات يساوي \(\varphi(j/d)\)، ومن ثم نحصل على

$$T(j)=\sum_{d\mid j} d\,\varphi\!\left(\frac{j}{d}\right).$$

وهذه هي الهوية الأساسية التي تزيل حساب gcd الصريح من الحلقة الداخلية.

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

الآن نجمع \(T(j)\) على جميع \(j\le N\):

$$G(N)=\sum_{j=1}^{N}\sum_{d\mid j} d\,\varphi\!\left(\frac{j}{d}\right).$$

إذا كتبنا \(j=d\,k\)، فإن كل زوج \((d,k)\) يحقق \(d\,k\le N\) يظهر مرة واحدة فقط، ولذلك

$$G(N)=\sum_{d=1}^{N} d\sum_{k\le N/d}\varphi(k)=\sum_{d=1}^{N} d\,\Phi\!\left(\left\lfloor\frac{N}{d}\right\rfloor\right).$$

وهكذا تصبح المسألة كلها عبارة عن عدد كبير من الاستعلامات عن \(\Phi(x)\).

الخطوة 3: اشتقاق علاقة عودية لمجموع دالة أويلر التراكمي

نستخدم الهوية المعروفة

$$\sum_{d\mid n}\varphi(d)=n,$$

وهي صحيحة لكل \(n\) موجب. إذا جمعناها من \(1\) حتى \(x\) نحصل على

$$\sum_{n=1}^{x} n=\sum_{n=1}^{x}\sum_{d\mid n}\varphi(d).$$

بعد تبديل ترتيب الجمع نلاحظ أن القاسم الثابت \(d\) يقسم بالضبط \(\left\lfloor x/d\right\rfloor\) عددًا لا يتجاوز \(x\)، لذلك

$$\frac{x(x+1)}{2}=\sum_{d=1}^{x}\varphi(d)\left\lfloor\frac{x}{d}\right\rfloor.$$

ثم نعيد تجميع الطرف الأيمن بحسب القيم المتساوية للقسمة الصحيحة فنحصل على

$$\frac{x(x+1)}{2}=\sum_{m=1}^{x}\Phi\!\left(\left\lfloor\frac{x}{m}\right\rfloor\right).$$

وبعزل الحد الموافق لـ \(m=1\) نحصل على العلاقة العودية

$$\Phi(x)=\frac{x(x+1)}{2}-\sum_{m=2}^{x}\Phi\!\left(\left\lfloor\frac{x}{m}\right\rfloor\right).$$

هذه هي الصيغة التي تستخدمها الشيفرة عندما يكون \(x\) أكبر من حد المعالجة المسبقة.

الخطوة 4: ضغط العلاقة باستخدام كتل القسمة الصحيحة

القيمة \(\left\lfloor x/m\right\rfloor\) لا تتغير عند كل \(m\). إذا وضعنا

$$q=\left\lfloor\frac{x}{\ell}\right\rfloor,\qquad r=\left\lfloor\frac{x}{q}\right\rfloor,$$

فإن جميع القيم \(m\in[\ell,r]\) تعطي نفس الخارج \(q\). لذلك يمكن ضغط العلاقة إلى

$$\Phi(x)=\frac{x(x+1)}{2}-\sum_{\text{blocks }[\ell,r]} (r-\ell+1)\,\Phi(q).$$

بدل المرور على كل \(m\) على حدة، نمر فقط على الكتل التي تغيّر فيها القسمة الصحيحة قيمتها. وعدد هذه الكتل يساوي \(O(\sqrt{x})\) عندما يكون \(x\) كبيرًا.

التنفيذات في C++ وPython وJava تحسب \(\varphi(n)\) ومجاميعها التراكمية مسبقًا لكل \(n\le 5\times 10^6\)، ثم تستخدم هذه العلاقة العودية المجمّعة مع التخزين التذكري فوق هذا الحد.

الخطوة 5: تطبيق الفكرة نفسها على المجموع الخارجي

التعبير النهائي المحوَّل

$$G(N)=\sum_{d=1}^{N} d\,\Phi\!\left(\left\lfloor\frac{N}{d}\right\rfloor\right)$$

له البنية نفسها. فإذا كانت \(\left\lfloor N/d\right\rfloor=q\) ثابتة على كتلة \(d\in[\ell,r]\)، فإن مساهمة الكتلة كلها هي

$$\left(\sum_{d=\ell}^{r} d\right)\Phi(q).$$

ومجموع المتتالية الحسابية هنا يساوي

$$\sum_{d=\ell}^{r} d=\frac{(\ell+r)(r-\ell+1)}{2}.$$

لذلك فإن الحساب النهائي لـ \(G(N)\) يتم هو الآخر بكتل قليلة نسبيًا بدل المرور على كل \(d\) حتى \(N\).

مثال محلول: \(N=10\)

تتحقق التنفيذات من القيمة الصغيرة \(G(10)=122\). ويمكن استرجاعها مباشرة من الصيغ السابقة.

نحتاج أولًا إلى القيم التالية للمجموع التراكمي:

$$\Phi(1)=1,\qquad \Phi(2)=2,\qquad \Phi(3)=4,\qquad \Phi(5)=10,\qquad \Phi(10)=32.$$

ثم نكتب

$$G(10)=\sum_{d=1}^{10} d\,\Phi\!\left(\left\lfloor\frac{10}{d}\right\rfloor\right).$$

القيم \(\left\lfloor 10/d\right\rfloor\) تبقى ثابتة على الكتل

$$[1,1],\qquad [2,2],\qquad [3,3],\qquad [4,5],\qquad [6,10],$$

وقيم الخارج المقابلة هي \(10,5,3,2,1\). لذلك

$$\begin{aligned} G(10)&=1\cdot \Phi(10)+2\cdot \Phi(5)+3\cdot \Phi(3)+(4+5)\Phi(2)+(6+7+8+9+10)\Phi(1)\\ &=1\cdot 32+2\cdot 10+3\cdot 4+9\cdot 2+40\cdot 1\\ &=32+20+12+18+40=122. \end{aligned}$$

هذا المثال الصغير يوضح بدقة لماذا يسمح تجميع القيم ذات الخارج نفسه بالتخلص من المسح الخطي الكامل.

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

تتبع تنفيذات C++ وPython وJava المسار نفسه. في البداية تُبنى قيم \(\varphi(n)\) حتى \(5\times 10^6\) باستخدام منخل خطي، ثم تُحوَّل إلى جدول تراكمي، وبذلك تصبح الاستعلامات الصغيرة عن \(\Phi(x)\) مجرد وصول مباشر إلى الجدول.

أما القيم الأكبر فتُحسب من العلاقة العودية السابقة بعد تجميع الحدود في كتل ذات قسمة صحيحة متساوية، مع تخزين كل نتيجة كبيرة عند أول حساب لها. هذا التخزين التذكري مهم لأن القيم الكبيرة نفسها تتكرر كثيرًا داخل الاستدعاءات العودية وفي الجمع النهائي.

بعد تجهيز آلية \(\Phi(x)\)، تحسب الشيفرة \(G(N)\) من خلال مرور آخر على كتل المتغيّر \(d\). في كل كتلة يُحسب مجموع المتتالية الحسابية، ثم يُضرب في القيمة المناسبة لـ \(\Phi(q)\)، ثم تُجمع النتيجة كلها بترديد \(998244353\). ولأن الترديد فردي، فإن القسمة على \(2\) تُنفَّذ على هيئة ضرب في المعكوس الضربي لـ \(2\) modulo هذا الترديد.

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

إذا كان حد المعالجة المسبقة هو \(L=5\times 10^6\)، فإن بناء المنخل والجدول التراكمي يكلف \(O(L)\) زمنًا و\(O(L)\) ذاكرة. وكل استعلام كبير عن \(\Phi(x)\) يُحسب مرة واحدة فقط مع التخزين التذكري، وتكلفته تتناسب مع عدد القيم المختلفة لـ \(\left\lfloor x/k\right\rfloor\)، أي \(O(\sqrt{x})\). كما أن المجموع الخارجي يُعالج هو الآخر بالكتل نفسها تقريبًا، لذلك يكون الزمن الكلي أقل بكثير من \(O(N)\) ومناسبًا لحجم \(N=10^{11}\).

المراجع

  1. صفحة المسألة: https://projecteuler.net/problem=625
  2. القاسم المشترك الأكبر: Wikipedia - Greatest common divisor
  3. دالة أويلر: Wikipedia - Euler's totient function
  4. التفاف ديريشليه: Wikipedia - Dirichlet convolution
  5. المنخل الخطي: cp-algorithms - Linear Sieve