Problem Summary

For each positive integer \(n\), consider the graph whose vertices are the positive divisors of \(n\). There is an upward edge from \(b\) to \(a\) exactly when \(a/b\) is prime, so each edge multiplies a divisor by one prime factor and nothing else. The quantity attached to \(n\) is the total edge weight

$$T(n)=\sum_{\substack{b\mid n,\ a\mid n\\ a/b\text{ prime}}}\bigl(\varphi(a)-\varphi(b)\bigr).$$

The problem asks for the prefix sum

$$S(N)=\sum_{n=1}^{N} T(n)\pmod{715827883},\qquad N=10^{12}.$$

A direct computation would have to build a divisor graph for every \(n\le N\), which is far too expensive. The implementations succeed because the graph sum collapses to a clean prime-power formula, and the global prefix can then be reorganized into prime-power loops plus a prime-prefix query structure.

Mathematical Approach

The key is to understand one fixed \(n\) exactly, and only then sum over all \(n\le N\).

The divisor graph splits by prime direction

Write

$$n=\prod_i p_i^{e_i}.$$

Every allowed edge in the graph is obtained by choosing one prime \(p\mid n\) and moving from a divisor \(d\) to \(pd\), provided \(pd\mid n\). So the whole graph decomposes into independent families of edges, one family for each prime dividing \(n\).

Fix one prime \(p\) with exact exponent \(p^e\parallel n\), and write

$$n=p^e m,\qquad p\nmid m.$$

Then every edge whose ratio is \(p\) has the form

$$p^r d \longrightarrow p^{r+1} d,$$

with \(0\le r<e\) and \(d\mid m\). So to find the total contribution of the \(p\)-direction, we only need to sum the totient differences on these edges.

Summing one prime direction exactly

Because \(p\nmid d\), Euler's totient function satisfies

$$\varphi(p^r d)=\varphi(p^r)\varphi(d).$$

There are two cases.

If \(r=0\), then

$$\varphi(pd)-\varphi(d)=(p-1)\varphi(d)-\varphi(d)=(p-2)\varphi(d).$$

If \(r\ge 1\), then

$$\varphi(p^{r+1}d)-\varphi(p^r d)=p^r(p-1)\varphi(d)-p^{r-1}(p-1)\varphi(d)=p^{r-1}(p-1)^2\varphi(d).$$

Therefore the full contribution of all \(p\)-edges is

$$\sum_{d\mid m}\varphi(d)\left[(p-2)+\sum_{r=1}^{e-1}p^{r-1}(p-1)^2\right].$$

The bracket simplifies to

$$\begin{aligned} (p-2)+(p-1)^2\sum_{r=1}^{e-1}p^{r-1} &=(p-2)+(p-1)^2\frac{p^{e-1}-1}{p-1} \\ &=(p-2)+(p-1)(p^{e-1}-1) \\ &=p^e-p^{e-1}-1. \end{aligned}$$

Now use the classical identity

$$\sum_{d\mid m}\varphi(d)=m.$$

So the entire \(p\)-direction contributes

$$m\bigl(p^e-p^{e-1}-1\bigr)=n-\frac{n}{p}-\frac{n}{p^e}.$$

Summing over the distinct prime powers of \(n\) gives the closed form used by the code:

$$\boxed{T(n)=\sum_{p^e\parallel n}\left(n-\frac{n}{p}-\frac{n}{p^e}\right).}$$

Worked example: \(n=12\)

The divisors of \(12\) are \(1,2,3,4,6,12\). The prime-ratio edges are

$$1\to2,\quad 1\to3,\quad 2\to4,\quad 2\to6,\quad 3\to6,\quad 4\to12,\quad 6\to12.$$

Their weights are

$$0,\ 1,\ 1,\ 1,\ 0,\ 2,\ 2,$$

because

$$\varphi(1)=1,\ \varphi(2)=1,\ \varphi(3)=2,\ \varphi(4)=2,\ \varphi(6)=2,\ \varphi(12)=4.$$

So \(T(12)=7\).

The closed form gives exactly the same result. Since \(12=2^2\cdot 3\),

$$T(12)=\left(12-\frac{12}{2}-\frac{12}{2^2}\right)+\left(12-\frac{12}{3}-\frac{12}{3}\right)=3+4=7.$$

This example shows what the formula is really doing: it adds the total contribution of each prime direction in the divisor graph.

Reordering the prefix sum by exact prime powers

To compute \(S(N)\), rewrite every \(n\) contributing to the term for \(p^e\) as

$$n=p^e m,\qquad p\nmid m.$$

For such an \(n\), the corresponding summand equals

$$n-\frac{n}{p}-\frac{n}{p^e}=m\bigl(p^e-p^{e-1}-1\bigr).$$

Hence

$$S(N)=\sum_{p^e\le N}\bigl(p^e-p^{e-1}-1\bigr)\sum_{\substack{m\le N/p^e\\ p\nmid m}} m.$$

The inner sum is simple once the multiples of \(p\) are removed. Define the triangular sum

$$\operatorname{Tri}(x)=\frac{x(x+1)}{2}.$$

Then

$$\sum_{\substack{m\le x\\ p\nmid m}} m=\operatorname{Tri}(x)-p\,\operatorname{Tri}\!\left(\left\lfloor\frac{x}{p}\right\rfloor\right),$$

because the excluded multiples are \(p,2p,\dots,p\lfloor x/p\rfloor\).

So the whole prefix sum becomes

$$\boxed{S(N)=\sum_{p^e\le N}\bigl(p^e-p^{e-1}-1\bigr)\left[\operatorname{Tri}\!\left(\left\lfloor\frac{N}{p^e}\right\rfloor\right)-p\,\operatorname{Tri}\!\left(\left\lfloor\frac{N}{p^{e+1}}\right\rfloor\right)\right].}$$

Why primes above \(\sqrt N\) can be batched together

The implementations split the summation at \(\sqrt N\).

If \(p\le \sqrt N\), the code can iterate over \(p,p^2,p^3,\dots\) directly.

If \(p>\sqrt N\), then \(p^2>N\), so only the exponent \(e=1\) can occur. In that case

$$x=\left\lfloor\frac{N}{p}\right\rfloor<p,$$

which means there is no positive multiple of \(p\) inside \(\{1,\dots,x\}\). Therefore

$$\sum_{\substack{m\le x\\ p\nmid m}}m=\operatorname{Tri}(x),$$

and the contribution of a large prime is simply

$$\bigl(p-2\bigr)\operatorname{Tri}\!\left(\left\lfloor\frac{N}{p}\right\rfloor\right).$$

Now group large primes by the constant quotient

$$q=\left\lfloor\frac{N}{p}\right\rfloor.$$

All primes in the interval

$$L=\left\lfloor\frac{N}{q+1}\right\rfloor+1,\qquad R=\left\lfloor\frac{N}{q}\right\rfloor$$

share this same \(q\), so their total contribution is

$$\operatorname{Tri}(q)\sum_{L\le p\le R}(p-2).$$

This reduces the large-prime side to interval queries for

$$\pi(x)=\#\{p\le x\},\qquad P(x)=\sum_{p\le x}p.$$

Once those two prime-prefix functions are available, every large-prime block can be evaluated as

$$\operatorname{Tri}(q)\Bigl(P(R)-P(L-1)-2\bigl(\pi(R)-\pi(L-1)\bigr)\Bigr).$$

How the Code Works

Validation comes first

The reference implementation does not assume the closed form blindly. It first builds the divisor graph directly for small \(n\), evaluates the edge sum from the definition, and checks that this agrees with the prime-power formula for all \(n\) in a small validation range. It also verifies several prefix values, including \(S(10)=26\) and \(S(100)=5282\), before moving on to the full \(N=10^{12}\) computation.

Compressed prime-prefix preprocessing

To answer many prime-count and prime-sum queries quickly, the implementation stores the distinct floor-quotients

$$\left\lfloor\frac{N}{1}\right\rfloor,\left\lfloor\frac{N}{2}\right\rfloor,\left\lfloor\frac{N}{3}\right\rfloor,\dots$$

instead of every integer up to \(N\). There are only \(O(\sqrt N)\) distinct values of this form. On that compressed coordinate set it initializes counts and sums of all integers at least \(2\), then removes composite contributions prime by prime. After that preprocessing, the code can query \(\pi(x)\) and \(P(x)\) for every relevant \(x\) in constant time.

Accumulating the answer

The final sum is assembled in two phases. First, every prime \(p\le\sqrt N\) is traversed through its powers \(p,p^2,p^3,\dots\), and the exact prime-power formula is added. Second, primes above \(\sqrt N\) are processed in quotient blocks with constant \(q=\lfloor N/p\rfloor\), using prime counts and prime sums on intervals.

All arithmetic is reduced modulo \(715827883\). The C++ implementation uses 128-bit intermediates before taking residues, the Python implementation relies on arbitrary-precision integers, and the Java entry point delegates the heavy computation to the compiled C++ solver and returns its result.

Complexity Analysis

The distinct values of \(\lfloor N/t\rfloor\) form an \(O(\sqrt N)\)-sized state space, and both the large-prime grouping and the query structure are built around that compression. Memory usage is therefore \(O(\sqrt N)\).

Time is dominated by the compressed prime-prefix sieve together with the sweep over primes up to \(\sqrt N\) and their powers. The important point is qualitative: the algorithm is far below linear in \(N\) and avoids any attempt to enumerate divisor graphs for all \(n\le 10^{12}\). That is what makes the computation practical.

Footnotes and References

  1. Problem page: https://projecteuler.net/problem=931
  2. Euler's totient function: Wikipedia - Euler's totient function
  3. Divisor: Wikipedia - Divisor
  4. Prime-counting function: Wikipedia - Prime-counting function
  5. Dirichlet hyperbola method: Wikipedia - Dirichlet hyperbola method

Problemzusammenfassung

Für jede positive ganze Zahl \(n\) betrachtet man einen Graphen, dessen Knoten die positiven Teiler von \(n\) sind. Es gibt genau dann eine aufsteigende Kante von \(b\) nach \(a\), wenn \(a/b\) eine Primzahl ist. Jede Kante entspricht also dem Multiplizieren eines Teilers mit genau einem Primfaktor. Die zu \(n\) gehörige Größe ist die gesamte Kantensumme

$$T(n)=\sum_{\substack{b\mid n,\ a\mid n\\ a/b\text{ prime}}}\bigl(\varphi(a)-\varphi(b)\bigr).$$

Gesucht ist dann die Präfixsumme

$$S(N)=\sum_{n=1}^{N} T(n)\pmod{715827883},\qquad N=10^{12}.$$

Eine direkte Berechnung würde für jedes \(n\le N\) den gesamten Teilergraphen erzeugen müssen und ist damit aussichtslos. Die Implementierungen funktionieren, weil sich die Graphensumme zunächst zu einer sauberen Formel über exakte Primpotenzen verdichtet und die globale Summe danach in Primpotenz-Schleifen und Primpräfix-Abfragen zerlegt werden kann.

Mathematischer Ansatz

Der erste Schritt ist, \(T(n)\) für ein festes \(n\) exakt zu verstehen. Danach wird erst über alle \(n\le N\) summiert.

Der Teilergraph zerfällt nach Primrichtungen

Schreibe

$$n=\prod_i p_i^{e_i}.$$

Jede erlaubte Kante entsteht, indem man eine Primzahl \(p\mid n\) auswählt und von einem Teiler \(d\) zu \(pd\) geht, solange \(pd\mid n\) bleibt. Der gesamte Graph zerfällt also in unabhängige Kantenfamilien, eine für jede Primzahl, die \(n\) teilt.

Fixiere eine Primzahl \(p\) mit exakter Potenz \(p^e\parallel n\) und schreibe

$$n=p^e m,\qquad p\nmid m.$$

Dann hat jede Kante mit Quotient \(p\) die Form

$$p^r d \longrightarrow p^{r+1} d,$$

wobei \(0\le r<e\) und \(d\mid m\) gilt. Um also den gesamten Beitrag der \(p\)-Richtung zu erhalten, genügt es, die Totientendifferenzen auf genau diesen Kanten zu summieren.

Eine Primrichtung exakt aufsummieren

Da \(p\nmid d\), gilt wegen der Multiplikativität von \(\varphi\)

$$\varphi(p^r d)=\varphi(p^r)\varphi(d).$$

Dabei entstehen zwei Fälle.

Für \(r=0\) ist

$$\varphi(pd)-\varphi(d)=(p-1)\varphi(d)-\varphi(d)=(p-2)\varphi(d).$$

Für \(r\ge 1\) erhält man

$$\varphi(p^{r+1}d)-\varphi(p^r d)=p^r(p-1)\varphi(d)-p^{r-1}(p-1)\varphi(d)=p^{r-1}(p-1)^2\varphi(d).$$

Somit ist der Gesamtbeitrag aller \(p\)-Kanten

$$\sum_{d\mid m}\varphi(d)\left[(p-2)+\sum_{r=1}^{e-1}p^{r-1}(p-1)^2\right].$$

Die Klammer vereinfacht sich zu

$$\begin{aligned} (p-2)+(p-1)^2\sum_{r=1}^{e-1}p^{r-1} &=(p-2)+(p-1)^2\frac{p^{e-1}-1}{p-1} \\ &=(p-2)+(p-1)(p^{e-1}-1) \\ &=p^e-p^{e-1}-1. \end{aligned}$$

Nun benutzt man die klassische Identität

$$\sum_{d\mid m}\varphi(d)=m.$$

Damit trägt die gesamte \(p\)-Richtung genau

$$m\bigl(p^e-p^{e-1}-1\bigr)=n-\frac{n}{p}-\frac{n}{p^e}$$

bei. Summiert man über alle exakten Primpotenzen von \(n\), erhält man die geschlossene Form

$$\boxed{T(n)=\sum_{p^e\parallel n}\left(n-\frac{n}{p}-\frac{n}{p^e}\right).}$$

Durchgerechnetes Beispiel: \(n=12\)

Die Teiler von \(12\) sind \(1,2,3,4,6,12\). Die Kanten mit Primquotient sind

$$1\to2,\quad 1\to3,\quad 2\to4,\quad 2\to6,\quad 3\to6,\quad 4\to12,\quad 6\to12.$$

Ihre Gewichte lauten

$$0,\ 1,\ 1,\ 1,\ 0,\ 2,\ 2,$$

denn

$$\varphi(1)=1,\ \varphi(2)=1,\ \varphi(3)=2,\ \varphi(4)=2,\ \varphi(6)=2,\ \varphi(12)=4.$$

Also ist \(T(12)=7\).

Die geschlossene Form liefert dasselbe Ergebnis. Wegen \(12=2^2\cdot 3\) gilt

$$T(12)=\left(12-\frac{12}{2}-\frac{12}{2^2}\right)+\left(12-\frac{12}{3}-\frac{12}{3}\right)=3+4=7.$$

Man sieht daran sehr klar, dass die Formel nichts anderes tut, als die gesamte Kantenmasse jeder Primrichtung zusammenzufassen.

Die Präfixsumme nach exakten Primpotenzen umordnen

Für die Berechnung von \(S(N)\) schreibt man jedes \(n\), das zum Term von \(p^e\) gehört, in der Form

$$n=p^e m,\qquad p\nmid m.$$

Dann ist der zugehörige Summand

$$n-\frac{n}{p}-\frac{n}{p^e}=m\bigl(p^e-p^{e-1}-1\bigr).$$

Folglich

$$S(N)=\sum_{p^e\le N}\bigl(p^e-p^{e-1}-1\bigr)\sum_{\substack{m\le N/p^e\\ p\nmid m}} m.$$

Die innere Summe wird einfach, sobald man die Vielfachen von \(p\) abzieht. Mit

$$\operatorname{Tri}(x)=\frac{x(x+1)}{2}$$

gilt

$$\sum_{\substack{m\le x\\ p\nmid m}} m=\operatorname{Tri}(x)-p\,\operatorname{Tri}\!\left(\left\lfloor\frac{x}{p}\right\rfloor\right),$$

weil die ausgeschlossenen Vielfachen genau \(p,2p,\dots,p\lfloor x/p\rfloor\) sind.

Somit wird die gesamte Präfixsumme zu

$$\boxed{S(N)=\sum_{p^e\le N}\bigl(p^e-p^{e-1}-1\bigr)\left[\operatorname{Tri}\!\left(\left\lfloor\frac{N}{p^e}\right\rfloor\right)-p\,\operatorname{Tri}\!\left(\left\lfloor\frac{N}{p^{e+1}}\right\rfloor\right)\right].}$$

Warum Primzahlen oberhalb von \(\sqrt N\) blockweise behandelt werden können

Die Implementierungen teilen die Summe bei \(\sqrt N\).

Für \(p\le\sqrt N\) kann man direkt über \(p,p^2,p^3,\dots\) iterieren.

Für \(p>\sqrt N\) gilt dagegen \(p^2>N\), also kommt nur der Exponent \(e=1\) überhaupt vor. Dann ist

$$x=\left\lfloor\frac{N}{p}\right\rfloor<p,$$

also liegt kein positives Vielfaches von \(p\) mehr in \(\{1,\dots,x\}\). Daher folgt

$$\sum_{\substack{m\le x\\ p\nmid m}}m=\operatorname{Tri}(x),$$

und der Beitrag einer großen Primzahl vereinfacht sich zu

$$\bigl(p-2\bigr)\operatorname{Tri}\!\left(\left\lfloor\frac{N}{p}\right\rfloor\right).$$

Nun gruppiert man große Primzahlen nach dem konstanten Quotienten

$$q=\left\lfloor\frac{N}{p}\right\rfloor.$$

Alle Primzahlen im Intervall

$$L=\left\lfloor\frac{N}{q+1}\right\rfloor+1,\qquad R=\left\lfloor\frac{N}{q}\right\rfloor$$

haben dasselbe \(q\). Ihr Gesamtbeitrag ist also

$$\operatorname{Tri}(q)\sum_{L\le p\le R}(p-2).$$

Damit reduziert sich die große-Primzahlen-Seite auf Intervallabfragen für

$$\pi(x)=\#\{p\le x\},\qquad P(x)=\sum_{p\le x}p.$$

Sobald diese beiden Primpräfixfunktionen verfügbar sind, kann jeder Block als

$$\operatorname{Tri}(q)\Bigl(P(R)-P(L-1)-2\bigl(\pi(R)-\pi(L-1)\bigr)\Bigr)$$

berechnet werden.

Wie der Code funktioniert

Zuerst wird die Formel validiert

Die Referenzimplementierung vertraut der geschlossenen Form nicht blind. Sie baut für kleine \(n\) den Teilergraphen direkt auf, berechnet die Kantensumme aus der Definition und prüft, dass diese mit der Primpotenzformel auf einem kleinen Testbereich übereinstimmt. Zusätzlich werden einige Präfixwerte überprüft, darunter \(S(10)=26\) und \(S(100)=5282\), bevor die eigentliche Rechnung für \(N=10^{12}\) startet.

Komprimierte Primpräfix-Vorbereitung

Um viele Primzahlanzahl- und Primzahlsummen-Abfragen schnell beantworten zu können, speichert die Implementierung nicht alle Werte bis \(N\), sondern nur die verschiedenen ganzzahligen Quotienten

$$\left\lfloor\frac{N}{1}\right\rfloor,\left\lfloor\frac{N}{2}\right\rfloor,\left\lfloor\frac{N}{3}\right\rfloor,\dots$$

Davon gibt es nur \(O(\sqrt N)\) verschiedene. Auf diesem komprimierten Koordinatensystem werden zunächst Anzahl und Summe aller ganzen Zahlen ab \(2\) initialisiert und anschließend die zusammengesetzten Beiträge Primzahl für Primzahl entfernt. Danach stehen \(\pi(x)\) und \(P(x)\) für alle relevanten \(x\) in konstanter Zeit zur Verfügung.

Die Antwort in zwei Phasen aufsummieren

Die Endsumme wird in zwei Phasen aufgebaut. Zuerst werden alle Primzahlen \(p\le\sqrt N\) zusammen mit ihren Potenzen \(p,p^2,p^3,\dots\) durchlaufen und der exakte Primpotenzterm addiert. Danach werden Primzahlen oberhalb von \(\sqrt N\) blockweise über konstante Quotienten \(q=\lfloor N/p\rfloor\) verarbeitet, wobei Primzahlanzahlen und Primzahlsummen auf Intervallen benutzt werden.

Alle Rechnungen werden modulo \(715827883\) reduziert. Die C++-Implementierung verwendet vor der Reduktion 128-Bit-Zwischenergebnisse, die Python-Implementierung verlässt sich auf beliebig große ganze Zahlen, und der Java-Einstiegspunkt delegiert die eigentliche Schwerarbeit an den kompilierten C++-Solver und gibt dessen Ergebnis zurück.

Komplexitätsanalyse

Die verschiedenen Werte von \(\lfloor N/t\rfloor\) bilden einen Zustandsraum der Größe \(O(\sqrt N)\). Sowohl die Gruppierung der großen Primzahlen als auch die Abfragestruktur bauen genau auf dieser Kompression auf. Der Speicherverbrauch ist deshalb \(O(\sqrt N)\).

Die Laufzeit wird von der komprimierten Primpräfix-Siebung sowie vom Durchlauf über die Primzahlen bis \(\sqrt N\) und deren Potenzen dominiert. Entscheidend ist die Größenordnung: Das Verfahren liegt deutlich unter linear in \(N\) und vermeidet vollständig den Versuch, für alle \(n\le 10^{12}\) einzelne Teilergraphen zu erzeugen. Genau dadurch wird die Rechnung praktikabel.

Fußnoten und Referenzen

  1. Problemseite: https://projecteuler.net/problem=931
  2. Eulersche Phi-Funktion: Wikipedia - Euler's totient function
  3. Teiler: Wikipedia - Divisor
  4. Primzahlzählfunktion: Wikipedia - Prime-counting function
  5. Dirichletsche Hyperbelmethode: Wikipedia - Dirichlet hyperbola method

Problem Özeti

Her pozitif \(n\) için, köşeleri \(n\)'in pozitif bölenleri olan bir grafik düşünülür. \(b\)'den \(a\)'ya tam olarak \(a/b\) bir asal sayı olduğunda bir yukarı kenar vardır. Yani her kenar, bir böleni yalnızca tek bir asal çarpanla büyütür. \(n\)'e ait nicelik şu kenar toplamıdır:

$$T(n)=\sum_{\substack{b\mid n,\ a\mid n\\ a/b\text{ prime}}}\bigl(\varphi(a)-\varphi(b)\bigr).$$

Sorulan şey ise şu önek toplamıdır:

$$S(N)=\sum_{n=1}^{N} T(n)\pmod{715827883},\qquad N=10^{12}.$$

Doğrudan yaklaşım, \(n\le N\) olan her sayı için bütün bölen grafiğini kurmayı gerektirir ve bu mümkün değildir. Uygulamalar başarılı olur, çünkü grafik üzerindeki toplam önce tam asal kuvvetler cinsinden kapalı bir forma indirgenir; ardından da küresel toplam, asal kuvvet döngüleri ve asal önek sorguları haline getirilir.

Matematiksel Yaklaşım

Ana fikir, önce sabit bir \(n\) için \(T(n)\)'i tam olarak çözmek, sonra bunu tüm \(n\le N\) üzerine toplamaktır.

Bölen grafiği asal yönlere ayrılır

Şöyle yazalım:

$$n=\prod_i p_i^{e_i}.$$

Grafikteki her geçerli kenar, \(p\mid n\) olan bir asal \(p\) seçip bir bölen \(d\)'den \(pd\)'ye gitmekle oluşur; tabii \(pd\mid n\) kalmalıdır. Dolayısıyla tüm grafik, \(n\)'i bölen her asal için bir tane olmak üzere bağımsız kenar ailelerine ayrılır.

Şimdi tam üs olarak \(p^e\parallel n\) olan bir asal \(p\) sabitleyelim ve

$$n=p^e m,\qquad p\nmid m$$

yazalım. O zaman oranı \(p\) olan her kenar

$$p^r d \longrightarrow p^{r+1} d$$

biçimindedir; burada \(0\le r<e\) ve \(d\mid m\) olur. Demek ki \(p\)-yönünün toplam katkısını bulmak için yalnızca bu kenarlar üzerindeki totient farklarını toplamak yeterlidir.

Tek bir asal yönünün toplamı

\(p\nmid d\) olduğundan \(\varphi\) fonksiyonu için

$$\varphi(p^r d)=\varphi(p^r)\varphi(d)$$

geçerlidir. Burada iki durum vardır.

\(r=0\) ise

$$\varphi(pd)-\varphi(d)=(p-1)\varphi(d)-\varphi(d)=(p-2)\varphi(d).$$

\(r\ge 1\) ise

$$\varphi(p^{r+1}d)-\varphi(p^r d)=p^r(p-1)\varphi(d)-p^{r-1}(p-1)\varphi(d)=p^{r-1}(p-1)^2\varphi(d).$$

Bu yüzden bütün \(p\)-kenarlarının toplam katkısı

$$\sum_{d\mid m}\varphi(d)\left[(p-2)+\sum_{r=1}^{e-1}p^{r-1}(p-1)^2\right]$$

olur. Köşeli parantez sadeleşince

$$\begin{aligned} (p-2)+(p-1)^2\sum_{r=1}^{e-1}p^{r-1} &=(p-2)+(p-1)^2\frac{p^{e-1}-1}{p-1} \\ &=(p-2)+(p-1)(p^{e-1}-1) \\ &=p^e-p^{e-1}-1 \end{aligned}$$

elde edilir. Şimdi klasik özdeşliği kullanırız:

$$\sum_{d\mid m}\varphi(d)=m.$$

Böylece tüm \(p\)-yönünün katkısı tam olarak

$$m\bigl(p^e-p^{e-1}-1\bigr)=n-\frac{n}{p}-\frac{n}{p^e}$$

olur. Tüm farklı asal kuvvetler üzerinde toplayınca kodun kullandığı kapalı form çıkar:

$$\boxed{T(n)=\sum_{p^e\parallel n}\left(n-\frac{n}{p}-\frac{n}{p^e}\right).}$$

Çalışılmış örnek: \(n=12\)

\(12\)'nin bölenleri \(1,2,3,4,6,12\)'dir. Asal oranlı kenarlar şunlardır:

$$1\to2,\quad 1\to3,\quad 2\to4,\quad 2\to6,\quad 3\to6,\quad 4\to12,\quad 6\to12.$$

Bu kenarların ağırlıkları

$$0,\ 1,\ 1,\ 1,\ 0,\ 2,\ 2$$

olur; çünkü

$$\varphi(1)=1,\ \varphi(2)=1,\ \varphi(3)=2,\ \varphi(4)=2,\ \varphi(6)=2,\ \varphi(12)=4.$$

Dolayısıyla \(T(12)=7\).

Kapalı form aynı sonucu verir. \(12=2^2\cdot 3\) olduğundan

$$T(12)=\left(12-\frac{12}{2}-\frac{12}{2^2}\right)+\left(12-\frac{12}{3}-\frac{12}{3}\right)=3+4=7.$$

Bu örnek, formülün gerçekte ne yaptığını gösterir: bölen grafiğinde her asal yönün toplam etkisini tek terimde toplar.

Önek toplamı tam asal kuvvetlere göre yeniden sıralamak

\(S(N)\)'i hesaplamak için, \(p^e\) terimine katkı veren her \(n\)'i

$$n=p^e m,\qquad p\nmid m$$

şeklinde yazarız. Bu durumda ilgili katkı

$$n-\frac{n}{p}-\frac{n}{p^e}=m\bigl(p^e-p^{e-1}-1\bigr)$$

olur. O halde

$$S(N)=\sum_{p^e\le N}\bigl(p^e-p^{e-1}-1\bigr)\sum_{\substack{m\le N/p^e\\ p\nmid m}} m.$$

İç toplam, \(p\)'nin katlarını çıkarınca basitleşir. Üçgensel toplamı

$$\operatorname{Tri}(x)=\frac{x(x+1)}{2}$$

diye tanımlarsak

$$\sum_{\substack{m\le x\\ p\nmid m}} m=\operatorname{Tri}(x)-p\,\operatorname{Tri}\!\left(\left\lfloor\frac{x}{p}\right\rfloor\right)$$

elde edilir; çünkü dışlanan sayılar tam olarak \(p,2p,\dots,p\lfloor x/p\rfloor\)'dir.

Böylece bütün önek toplamı

$$\boxed{S(N)=\sum_{p^e\le N}\bigl(p^e-p^{e-1}-1\bigr)\left[\operatorname{Tri}\!\left(\left\lfloor\frac{N}{p^e}\right\rfloor\right)-p\,\operatorname{Tri}\!\left(\left\lfloor\frac{N}{p^{e+1}}\right\rfloor\right)\right].}$$

\(\sqrt N\)'in üzerindeki asallar neden bloklanabiliyor?

Uygulamalar toplamı \(\sqrt N\) noktasında ikiye böler.

\(p\le\sqrt N\) için \(p,p^2,p^3,\dots\) doğrudan dolaşılabilir.

\(p>\sqrt N\) için ise \(p^2>N\) olur; dolayısıyla yalnızca \(e=1\) mümkündür. Bu durumda

$$x=\left\lfloor\frac{N}{p}\right\rfloor<p$$

olduğundan \(1\) ile \(x\) arasında \(p\)'nin hiçbir pozitif katı yoktur. Demek ki

$$\sum_{\substack{m\le x\\ p\nmid m}}m=\operatorname{Tri}(x),$$

ve büyük bir asalın katkısı

$$\bigl(p-2\bigr)\operatorname{Tri}\!\left(\left\lfloor\frac{N}{p}\right\rfloor\right)$$

şeklinde sadeleşir.

Şimdi büyük asalları sabit bölüm

$$q=\left\lfloor\frac{N}{p}\right\rfloor$$

değerine göre gruplayalım. Şu aralıktaki tüm asallar

$$L=\left\lfloor\frac{N}{q+1}\right\rfloor+1,\qquad R=\left\lfloor\frac{N}{q}\right\rfloor$$

aynı \(q\) değerini verir. O halde bu bloğun toplamı

$$\operatorname{Tri}(q)\sum_{L\le p\le R}(p-2)$$

olur. Yani büyük-asal tarafı, şu iki önek fonksiyonunun aralık sorgularına indirgenir:

$$\pi(x)=\#\{p\le x\},\qquad P(x)=\sum_{p\le x}p.$$

Bu iki değer elde edilince her blok

$$\operatorname{Tri}(q)\Bigl(P(R)-P(L-1)-2\bigl(\pi(R)-\pi(L-1)\bigr)\Bigr)$$

olarak hesaplanır.

Kod Nasıl Çalışır

Önce doğrulama yapılır

Referans uygulama kapalı forma körü körüne güvenmez. Küçük \(n\) değerleri için bölen grafiğini doğrudan kurar, tanımdaki kenar toplamını hesaplar ve bunun asal-kuvvet formülüyle küçük bir aralıkta birebir uyuştuğunu denetler. Ayrıca tam hesaplamaya geçmeden önce \(S(10)=26\) ve \(S(100)=5282\) gibi bazı önek değerleri de doğrular.

Sıkıştırılmış asal önek ön işlemi

Çok sayıda asal sayma ve asal toplamı sorgusunu hızlı cevaplamak için uygulama \(N\)'e kadar bütün tam sayıları tutmaz; bunun yerine yalnızca şu farklı bölüm değerlerini saklar:

$$\left\lfloor\frac{N}{1}\right\rfloor,\left\lfloor\frac{N}{2}\right\rfloor,\left\lfloor\frac{N}{3}\right\rfloor,\dots$$

Bunların sayısı yalnızca \(O(\sqrt N)\)'dir. Bu sıkıştırılmış koordinat kümesinde önce \(2\) ve üzerindeki bütün tamsayıların sayıları ve toplamları başlatılır, sonra bileşik sayı katkıları asal asal ayıklanır. Bu hazırlıktan sonra ilgili tüm \(x\) değerleri için \(\pi(x)\) ve \(P(x)\) sabit zamanda sorgulanabilir.

Cevap iki fazda toplanır

Son toplam iki aşamada kurulur. İlk olarak \(p\le\sqrt N\) olan tüm asallar ve onların \(p,p^2,p^3,\dots\) kuvvetleri dolaşılıp tam asal-kuvvet terimi eklenir. İkinci olarak \(\sqrt N\)'in üzerindeki asallar, sabit \(q=\lfloor N/p\rfloor\) bölüm blokları halinde işlenir ve aralıklardaki asal sayıları ile asal toplamları kullanılır.

Tüm aritmetik \(715827883\) modunda tutulur. C++ uygulaması kalanı almadan önce 128 bit ara değerler kullanır, Python uygulaması sınırsız büyüklükte tamsayılara yaslanır, Java giriş noktası ise ağır hesabı derlenmiş C++ çözücüsüne devredip onun sonucunu döndürür.

Karmaşıklık Analizi

\(\lfloor N/t\rfloor\) değerlerinin farklı halleri \(O(\sqrt N)\) büyüklüğünde bir durum uzayı oluşturur. Hem büyük asalların bloklanması hem de sorgu yapısı bu sıkıştırmaya dayanır. Bu yüzden bellek kullanımı \(O(\sqrt N)\)'dir.

Çalışma süresinin baskın kısmı sıkıştırılmış asal-önek eleği ile \(\sqrt N\)'e kadar olan asalların ve onların kuvvetlerinin taranmasıdır. Buradaki kritik nokta nitelikseldir: algoritma \(N\)'e göre doğrusal olmaktan çok daha iyidir ve \(n\le 10^{12}\) için tek tek bölen grafikleri kurma girişiminden tamamen kaçınır. Pratik olmasını sağlayan şey budur.

Dipnotlar ve Kaynaklar

  1. Problem sayfası: https://projecteuler.net/problem=931
  2. Euler totient fonksiyonu: Wikipedia - Euler's totient function
  3. Bölen: Wikipedia - Divisor
  4. Asal sayma fonksiyonu: Wikipedia - Prime-counting function
  5. Dirichlet hiperbol yöntemi: Wikipedia - Dirichlet hyperbola method

Resumen del Problema

Para cada entero positivo \(n\), se considera el grafo cuyos vértices son los divisores positivos de \(n\). Hay una arista ascendente de \(b\) hacia \(a\) exactamente cuando \(a/b\) es primo, de modo que cada arista multiplica un divisor por un solo factor primo y nada más. La cantidad asociada a \(n\) es la suma total de pesos de aristas

$$T(n)=\sum_{\substack{b\mid n,\ a\mid n\\ a/b\text{ prime}}}\bigl(\varphi(a)-\varphi(b)\bigr).$$

El problema pide la suma prefijo

$$S(N)=\sum_{n=1}^{N} T(n)\pmod{715827883},\qquad N=10^{12}.$$

Un ataque directo tendría que construir un grafo de divisores distinto para cada \(n\le N\), lo cual es inviable. Las implementaciones funcionan porque primero reducen la suma del grafo a una fórmula cerrada sobre potencias primas exactas, y después reorganizan la suma global en bucles sobre potencias primas y consultas a prefijos de primos.

Enfoque Matemático

La idea central es entender \(T(n)\) exactamente para un \(n\) fijo y solo después sumar sobre todos los \(n\le N\).

El grafo de divisores se separa por direcciones primas

Escribamos

$$n=\prod_i p_i^{e_i}.$$

Cada arista permitida del grafo se obtiene eligiendo un primo \(p\mid n\) y pasando de un divisor \(d\) a \(pd\), siempre que \(pd\mid n\). Por tanto, el grafo completo se descompone en familias independientes de aristas, una por cada primo que divide a \(n\).

Fijemos ahora un primo \(p\) con exponente exacto \(p^e\parallel n\), y escribamos

$$n=p^e m,\qquad p\nmid m.$$

Entonces toda arista cuya razón es \(p\) tiene la forma

$$p^r d \longrightarrow p^{r+1} d,$$

con \(0\le r<e\) y \(d\mid m\). Así, para obtener la contribución total de la dirección \(p\), basta sumar las diferencias de totientes en esas aristas.

Sumar exactamente una dirección prima

Como \(p\nmid d\), la función totiente satisface

$$\varphi(p^r d)=\varphi(p^r)\varphi(d).$$

Hay dos casos.

Si \(r=0\), entonces

$$\varphi(pd)-\varphi(d)=(p-1)\varphi(d)-\varphi(d)=(p-2)\varphi(d).$$

Si \(r\ge 1\), entonces

$$\varphi(p^{r+1}d)-\varphi(p^r d)=p^r(p-1)\varphi(d)-p^{r-1}(p-1)\varphi(d)=p^{r-1}(p-1)^2\varphi(d).$$

Por consiguiente, la contribución total de todas las aristas de tipo \(p\) es

$$\sum_{d\mid m}\varphi(d)\left[(p-2)+\sum_{r=1}^{e-1}p^{r-1}(p-1)^2\right].$$

El corchete se simplifica a

$$\begin{aligned} (p-2)+(p-1)^2\sum_{r=1}^{e-1}p^{r-1} &=(p-2)+(p-1)^2\frac{p^{e-1}-1}{p-1} \\ &=(p-2)+(p-1)(p^{e-1}-1) \\ &=p^e-p^{e-1}-1. \end{aligned}$$

Ahora usamos la identidad clásica

$$\sum_{d\mid m}\varphi(d)=m.$$

Así, la dirección completa asociada a \(p\) aporta

$$m\bigl(p^e-p^{e-1}-1\bigr)=n-\frac{n}{p}-\frac{n}{p^e}.$$

Al sumar sobre todas las potencias primas exactas de \(n\), obtenemos la forma cerrada que usan las implementaciones:

$$\boxed{T(n)=\sum_{p^e\parallel n}\left(n-\frac{n}{p}-\frac{n}{p^e}\right).}$$

Ejemplo trabajado: \(n=12\)

Los divisores de \(12\) son \(1,2,3,4,6,12\). Las aristas de razón prima son

$$1\to2,\quad 1\to3,\quad 2\to4,\quad 2\to6,\quad 3\to6,\quad 4\to12,\quad 6\to12.$$

Sus pesos son

$$0,\ 1,\ 1,\ 1,\ 0,\ 2,\ 2,$$

porque

$$\varphi(1)=1,\ \varphi(2)=1,\ \varphi(3)=2,\ \varphi(4)=2,\ \varphi(6)=2,\ \varphi(12)=4.$$

Por tanto, \(T(12)=7\).

La fórmula cerrada produce exactamente lo mismo. Como \(12=2^2\cdot 3\),

$$T(12)=\left(12-\frac{12}{2}-\frac{12}{2^2}\right)+\left(12-\frac{12}{3}-\frac{12}{3}\right)=3+4=7.$$

Este ejemplo deja claro lo que está ocurriendo: la fórmula suma la contribución total de cada dirección prima del grafo.

Reordenar la suma prefijo por potencias primas exactas

Para calcular \(S(N)\), reescribimos cada \(n\) que contribuye al término de \(p^e\) como

$$n=p^e m,\qquad p\nmid m.$$

En ese caso, el sumando correspondiente vale

$$n-\frac{n}{p}-\frac{n}{p^e}=m\bigl(p^e-p^{e-1}-1\bigr).$$

Por ello,

$$S(N)=\sum_{p^e\le N}\bigl(p^e-p^{e-1}-1\bigr)\sum_{\substack{m\le N/p^e\\ p\nmid m}} m.$$

La suma interior se vuelve sencilla al eliminar los múltiplos de \(p\). Definimos

$$\operatorname{Tri}(x)=\frac{x(x+1)}{2}.$$

Entonces

$$\sum_{\substack{m\le x\\ p\nmid m}} m=\operatorname{Tri}(x)-p\,\operatorname{Tri}\!\left(\left\lfloor\frac{x}{p}\right\rfloor\right),$$

porque los términos excluidos son exactamente \(p,2p,\dots,p\lfloor x/p\rfloor\).

Así, la suma prefijo completa queda

$$\boxed{S(N)=\sum_{p^e\le N}\bigl(p^e-p^{e-1}-1\bigr)\left[\operatorname{Tri}\!\left(\left\lfloor\frac{N}{p^e}\right\rfloor\right)-p\,\operatorname{Tri}\!\left(\left\lfloor\frac{N}{p^{e+1}}\right\rfloor\right)\right].}$$

Por qué los primos mayores que \(\sqrt N\) se agrupan por bloques

Las implementaciones dividen la suma en \(\sqrt N\).

Si \(p\le\sqrt N\), se puede iterar directamente sobre \(p,p^2,p^3,\dots\).

Si \(p>\sqrt N\), entonces \(p^2>N\), así que solo puede aparecer el exponente \(e=1\). En ese caso

$$x=\left\lfloor\frac{N}{p}\right\rfloor<p,$$

lo que significa que no existe ningún múltiplo positivo de \(p\) dentro de \(\{1,\dots,x\}\). Por tanto

$$\sum_{\substack{m\le x\\ p\nmid m}}m=\operatorname{Tri}(x),$$

y la contribución de un primo grande se reduce a

$$\bigl(p-2\bigr)\operatorname{Tri}\!\left(\left\lfloor\frac{N}{p}\right\rfloor\right).$$

Ahora agrupamos los primos grandes por el cociente constante

$$q=\left\lfloor\frac{N}{p}\right\rfloor.$$

Todos los primos del intervalo

$$L=\left\lfloor\frac{N}{q+1}\right\rfloor+1,\qquad R=\left\lfloor\frac{N}{q}\right\rfloor$$

comparten ese mismo \(q\), así que su contribución total es

$$\operatorname{Tri}(q)\sum_{L\le p\le R}(p-2).$$

Con ello, la parte de primos grandes se reduce a consultas por intervalos de

$$\pi(x)=\#\{p\le x\},\qquad P(x)=\sum_{p\le x}p.$$

Una vez disponibles esas dos funciones prefijo, cada bloque se evalúa como

$$\operatorname{Tri}(q)\Bigl(P(R)-P(L-1)-2\bigl(\pi(R)-\pi(L-1)\bigr)\Bigr).$$

Cómo Funciona el Código

Primero se valida la identidad

La implementación de referencia no da por buena la fórmula cerrada sin comprobarla. Para valores pequeños de \(n\), construye el grafo de divisores de manera directa, evalúa la suma de aristas desde la definición y verifica que coincide con la fórmula por potencias primas en un rango de prueba. También comprueba varios valores prefijo, entre ellos \(S(10)=26\) y \(S(100)=5282\), antes de pasar al cálculo completo para \(N=10^{12}\).

Preprocesamiento comprimido de prefijos de primos

Para responder rápidamente muchas consultas de conteo de primos y suma de primos, la implementación no almacena todos los enteros hasta \(N\), sino solo los cocientes enteros distintos

$$\left\lfloor\frac{N}{1}\right\rfloor,\left\lfloor\frac{N}{2}\right\rfloor,\left\lfloor\frac{N}{3}\right\rfloor,\dots$$

Solo hay \(O(\sqrt N)\) valores distintos de esta forma. Sobre ese conjunto comprimido se inicializan los conteos y sumas de todos los enteros al menos \(2\), y luego se eliminan las contribuciones compuestas primo por primo. Después de ese preprocesamiento, \(\pi(x)\) y \(P(x)\) pueden consultarse en tiempo constante para todos los \(x\) relevantes.

Acumulación final en dos fases

La respuesta se arma en dos etapas. Primero se recorren todos los primos \(p\le\sqrt N\) junto con sus potencias \(p,p^2,p^3,\dots\), y se añade el término exacto de potencias primas. Después, los primos mayores que \(\sqrt N\) se procesan por bloques de cociente constante \(q=\lfloor N/p\rfloor\), usando conteos y sumas de primos sobre intervalos.

Toda la aritmética se reduce módulo \(715827883\). La versión en C++ usa enteros intermedios de 128 bits antes de reducir, la versión en Python se apoya en enteros de precisión arbitraria, y la entrada en Java delega el cálculo pesado al solucionador compilado en C++ y devuelve su resultado.

Análisis de Complejidad

Los valores distintos de \(\lfloor N/t\rfloor\) forman un espacio de estados de tamaño \(O(\sqrt N)\). Tanto la agrupación de primos grandes como la estructura de consultas se construyen sobre esa compresión. Por ello, el uso de memoria es \(O(\sqrt N)\).

El tiempo está dominado por el tamiz comprimido de prefijos de primos y por el recorrido de los primos hasta \(\sqrt N\) junto con sus potencias. Lo importante es la escala: el algoritmo queda muy por debajo de lineal en \(N\) y evita por completo cualquier intento de construir grafos de divisores para todos los \(n\le 10^{12}\). Eso es lo que lo hace práctico.

Notas y Referencias

  1. Página del problema: https://projecteuler.net/problem=931
  2. Función phi de Euler: Wikipedia - Euler's totient function
  3. Divisor: Wikipedia - Divisor
  4. Función contadora de primos: Wikipedia - Prime-counting function
  5. Método de la hipérbola de Dirichlet: Wikipedia - Dirichlet hyperbola method

问题概述

对每个正整数 \(n\),考虑一个以 \(n\) 的所有正因子为顶点的图。如果两个因子 \(b\) 与 \(a\) 满足 \(a/b\) 是质数,那么就从 \(b\) 向 \(a\) 连一条向上的边。也就是说,每一条边都只是在原来的因子上额外乘进一个质因子。与 \(n\) 对应的量定义为所有这些边权之和:

$$T(n)=\sum_{\substack{b\mid n,\ a\mid n\\ a/b\text{ prime}}}\bigl(\varphi(a)-\varphi(b)\bigr).$$

题目最终要求的是前缀和

$$S(N)=\sum_{n=1}^{N} T(n)\pmod{715827883},\qquad N=10^{12}.$$

如果直接做,就意味着要为每个 \(n\le N\) 单独建立它的因子图并枚举边,这显然不可行。程序之所以能处理 \(10^{12}\),是因为它先把单个 \(T(n)\) 的图论定义压缩成一个精确的素幂公式,再把整个前缀和改写成“小质数的素幂遍历 + 大质数的前缀查询”。

数学方法

真正的突破口是:先把固定 \(n\) 的结构完全看清楚,再去做对所有 \(n\le N\) 的求和。

因子图会按质数方向自然分解

把 \(n\) 写成

$$n=\prod_i p_i^{e_i}.$$

图中的每一条合法边都来自下面这个操作:选定一个质数 \(p\mid n\),然后把某个因子 \(d\) 变成 \(pd\),前提是 \(pd\) 仍然是 \(n\) 的因子。因此,整张图可以按“是哪一个质数提供了这条边”拆成若干互不干扰的方向。

固定其中一个质数 \(p\),并设它在 \(n\) 中的精确指数为 \(p^e\parallel n\)。写成

$$n=p^e m,\qquad p\nmid m.$$

那么所有比值等于 \(p\) 的边都恰好具有形式

$$p^r d \longrightarrow p^{r+1} d,$$

其中 \(0\le r<e\),并且 \(d\mid m\)。所以只要把这些边上的 \(\varphi\) 差值全都加起来,就得到了质数 \(p\) 这一整个方向的总贡献。

把单个质数方向完整求和

由于这里 \(p\nmid d\),欧拉函数满足

$$\varphi(p^r d)=\varphi(p^r)\varphi(d).$$

于是会出现两个明显不同的情况。

当 \(r=0\) 时,边是 \(d\to pd\),所以

$$\varphi(pd)-\varphi(d)=(p-1)\varphi(d)-\varphi(d)=(p-2)\varphi(d).$$

当 \(r\ge 1\) 时,边是 \(p^r d\to p^{r+1}d\),于是

$$\varphi(p^{r+1}d)-\varphi(p^r d)=p^r(p-1)\varphi(d)-p^{r-1}(p-1)\varphi(d)=p^{r-1}(p-1)^2\varphi(d).$$

因此,所有 \(p\)-方向边的总贡献等于

$$\sum_{d\mid m}\varphi(d)\left[(p-2)+\sum_{r=1}^{e-1}p^{r-1}(p-1)^2\right].$$

括号里的表达式可以化简为

$$\begin{aligned} (p-2)+(p-1)^2\sum_{r=1}^{e-1}p^{r-1} &=(p-2)+(p-1)^2\frac{p^{e-1}-1}{p-1} \\ &=(p-2)+(p-1)(p^{e-1}-1) \\ &=p^e-p^{e-1}-1. \end{aligned}$$

接着使用经典恒等式

$$\sum_{d\mid m}\varphi(d)=m.$$

于是整个 \(p\)-方向的贡献恰好是

$$m\bigl(p^e-p^{e-1}-1\bigr)=n-\frac{n}{p}-\frac{n}{p^e}.$$

对 \(n\) 的所有精确素幂 \(p^e\parallel n\) 求和,就得到程序真正依赖的闭式:

$$\boxed{T(n)=\sum_{p^e\parallel n}\left(n-\frac{n}{p}-\frac{n}{p^e}\right).}$$

完整例子:\(n=12\)

\(12\) 的正因子是 \(1,2,3,4,6,12\)。满足“商为质数”的边共有

$$1\to2,\quad 1\to3,\quad 2\to4,\quad 2\to6,\quad 3\to6,\quad 4\to12,\quad 6\to12.$$

这些边的权值分别是

$$0,\ 1,\ 1,\ 1,\ 0,\ 2,\ 2,$$

因为

$$\varphi(1)=1,\ \varphi(2)=1,\ \varphi(3)=2,\ \varphi(4)=2,\ \varphi(6)=2,\ \varphi(12)=4.$$

所以 \(T(12)=7\)。

闭式也给出同样的答案。由于 \(12=2^2\cdot 3\),所以

$$T(12)=\left(12-\frac{12}{2}-\frac{12}{2^2}\right)+\left(12-\frac{12}{3}-\frac{12}{3}\right)=3+4=7.$$

这个例子非常有代表性:闭式并没有做什么神秘的事,它只是把因子图中每一个质数方向的总边权提前汇总好了。

把前缀和改写成对精确素幂的求和

现在开始计算

$$S(N)=\sum_{n\le N}T(n).$$

固定一个素幂 \(p^e\),所有对它有贡献的 \(n\) 都可以唯一写成

$$n=p^e m,\qquad p\nmid m.$$

对这样的 \(n\),对应项正是

$$n-\frac{n}{p}-\frac{n}{p^e}=m\bigl(p^e-p^{e-1}-1\bigr).$$

因此整个前缀和可以重排为

$$S(N)=\sum_{p^e\le N}\bigl(p^e-p^{e-1}-1\bigr)\sum_{\substack{m\le N/p^e\\ p\nmid m}} m.$$

内层和只是在求“不被 \(p\) 整除的整数之和”。定义三角和

$$\operatorname{Tri}(x)=\frac{x(x+1)}{2}.$$

那么

$$\sum_{\substack{m\le x\\ p\nmid m}} m=\operatorname{Tri}(x)-p\,\operatorname{Tri}\!\left(\left\lfloor\frac{x}{p}\right\rfloor\right),$$

因为被删掉的恰好是 \(p,2p,\dots,p\lfloor x/p\rfloor\)。于是得到最终的重排公式:

$$\boxed{S(N)=\sum_{p^e\le N}\bigl(p^e-p^{e-1}-1\bigr)\left[\operatorname{Tri}\!\left(\left\lfloor\frac{N}{p^e}\right\rfloor\right)-p\,\operatorname{Tri}\!\left(\left\lfloor\frac{N}{p^{e+1}}\right\rfloor\right)\right].}$$

为什么 \(\sqrt N\) 以上的质数可以整块处理

程序的实现把求和自然地分成“小质数部分”和“大质数部分”。

如果 \(p\le\sqrt N\),那么 \(p,p^2,p^3,\dots\) 这些幂都可能出现,所以直接逐个素幂遍历是合理的。

如果 \(p>\sqrt N\),那么 \(p^2>N\),也就是说此时只可能有指数 \(e=1\)。这时

$$x=\left\lfloor\frac{N}{p}\right\rfloor<p,$$

因此在 \(1\) 到 \(x\) 的范围里根本没有正的 \(p\) 倍数,于是

$$\sum_{\substack{m\le x\\ p\nmid m}}m=\operatorname{Tri}(x).$$

所以每个大质数的贡献会简化为

$$\bigl(p-2\bigr)\operatorname{Tri}\!\left(\left\lfloor\frac{N}{p}\right\rfloor\right).$$

接下来按固定商

$$q=\left\lfloor\frac{N}{p}\right\rfloor$$

来分组。所有满足

$$L=\left\lfloor\frac{N}{q+1}\right\rfloor+1,\qquad R=\left\lfloor\frac{N}{q}\right\rfloor$$

的区间 \([L,R]\) 中的质数,都对应同一个 \(q\)。于是这整个区间的贡献可以写成

$$\operatorname{Tri}(q)\sum_{L\le p\le R}(p-2).$$

这样一来,大质数部分就完全降成了两个前缀函数的区间查询:

$$\pi(x)=\#\{p\le x\},\qquad P(x)=\sum_{p\le x}p.$$

一旦这两个函数能快速求值,每个区间块就能写成

$$\operatorname{Tri}(q)\Bigl(P(R)-P(L-1)-2\bigl(\pi(R)-\pi(L-1)\bigr)\Bigr).$$

代码如何工作

先验证公式,再跑大规模计算

参考实现并不是直接把闭式写进去就结束了。它先对较小的 \(n\) 真的构造因子图,按定义把所有边权 \(\varphi(a)-\varphi(b)\) 加起来,然后检查这个结果是否与闭式在一段小范围内完全一致。接着它还会验证若干前缀值,包括 \(S(10)=26\) 和 \(S(100)=5282\),确认数学化简没有问题之后,才开始正式计算 \(N=10^{12}\)。

用压缩后的商值集合建立质数前缀

程序需要频繁查询“\(\le x\) 的质数个数”和“\(\le x\) 的质数和”。如果对所有整数都建表,规模会太大,所以实现只保存所有不同的

$$\left\lfloor\frac{N}{1}\right\rfloor,\left\lfloor\frac{N}{2}\right\rfloor,\left\lfloor\frac{N}{3}\right\rfloor,\dots$$

这些商值。这样的不同取值只有 \(O(\sqrt N)\) 个。然后在这个压缩坐标系上,先把所有不小于 \(2\) 的整数计数与求和初始化,再逐个质数剔除合数贡献。预处理完成后,\(\pi(x)\) 与 \(P(x)\) 对所有需要的 \(x\) 都能常数时间查询。

最终答案分两段累加

第一段处理所有 \(p\le\sqrt N\) 的质数,并沿着 \(p,p^2,p^3,\dots\) 这些素幂直接应用精确公式。第二段处理所有 \(p>\sqrt N\) 的质数,但不是逐个做,而是按商 \(q=\lfloor N/p\rfloor\) 相同的区间分块,再通过质数个数和质数前缀和来一次算完整块。

所有运算都在模 \(715827883\) 下进行。C++ 实现先用 128 位整数承接中间值再取模,Python 实现依靠任意精度整数,而 Java 入口并不重写整套算法,而是调用已编译的 C++ 求解器并返回解析后的结果。

复杂度分析

所有不同的 \(\lfloor N/t\rfloor\) 构成了一个 \(O(\sqrt N)\) 大小的状态空间,质数前缀结构和大质数分块都建立在这层压缩之上,所以空间复杂度是 \(O(\sqrt N)\)。

时间主要消耗在压缩质数前缀筛和对 \(\sqrt N\) 以内质数及其素幂的遍历上。这里最重要的不是某个过于细的渐近常数,而是整体规模已经远低于线性于 \(N\) 的做法,更不会去为所有 \(n\le 10^{12}\) 构造因子图。正因为如此,这个问题才变得可算。

注释与参考资料

  1. 题目页面:https://projecteuler.net/problem=931
  2. 欧拉函数:Wikipedia - Euler's totient function
  3. 因子:Wikipedia - Divisor
  4. 质数计数函数:Wikipedia - Prime-counting function
  5. Dirichlet 双曲线方法:Wikipedia - Dirichlet hyperbola method

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

Для каждого положительного целого \(n\) рассматривается граф, вершинами которого являются все положительные делители числа \(n\). Из \(b\) в \(a\) проводится ребро тогда и только тогда, когда \(a/b\) является простым числом. Иначе говоря, каждое ребро соответствует умножению делителя ровно на один простой множитель. Величина, связанная с \(n\), задается суммой весов всех таких ребер:

$$T(n)=\sum_{\substack{b\mid n,\ a\mid n\\ a/b\text{ prime}}}\bigl(\varphi(a)-\varphi(b)\bigr).$$

Требуется вычислить префиксную сумму

$$S(N)=\sum_{n=1}^{N} T(n)\pmod{715827883},\qquad N=10^{12}.$$

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

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

Основной ход состоит в том, чтобы сначала полностью понять \(T(n)\) для фиксированного \(n\), а уже потом суммировать по всем \(n\le N\).

Граф делителей раскладывается по простым направлениям

Запишем

$$n=\prod_i p_i^{e_i}.$$

Каждое допустимое ребро получается так: выбирается простой делитель \(p\mid n\), после чего некоторый делитель \(d\) заменяется на \(pd\), если \(pd\) по-прежнему делит \(n\). Значит, весь граф естественно распадается на независимые семейства ребер, по одному семейству на каждый простой делитель числа \(n\).

Зафиксируем простой \(p\) с точной степенью \(p^e\parallel n\) и представим \(n\) в виде

$$n=p^e m,\qquad p\nmid m.$$

Тогда любое ребро с отношением \(p\) имеет вид

$$p^r d \longrightarrow p^{r+1} d,$$

где \(0\le r<e\) и \(d\mid m\). Следовательно, чтобы найти полный вклад направления \(p\), нужно просто просуммировать разности значений функции Эйлера на этих ребрах.

Точный суммарный вклад одного простого направления

Поскольку \(p\nmid d\), выполняется мультипликативная формула

$$\varphi(p^r d)=\varphi(p^r)\varphi(d).$$

Далее возникают два случая.

Если \(r=0\), то

$$\varphi(pd)-\varphi(d)=(p-1)\varphi(d)-\varphi(d)=(p-2)\varphi(d).$$

Если \(r\ge 1\), то

$$\varphi(p^{r+1}d)-\varphi(p^r d)=p^r(p-1)\varphi(d)-p^{r-1}(p-1)\varphi(d)=p^{r-1}(p-1)^2\varphi(d).$$

Поэтому суммарный вклад всех ребер типа \(p\) равен

$$\sum_{d\mid m}\varphi(d)\left[(p-2)+\sum_{r=1}^{e-1}p^{r-1}(p-1)^2\right].$$

Скобка упрощается так:

$$\begin{aligned} (p-2)+(p-1)^2\sum_{r=1}^{e-1}p^{r-1} &=(p-2)+(p-1)^2\frac{p^{e-1}-1}{p-1} \\ &=(p-2)+(p-1)(p^{e-1}-1) \\ &=p^e-p^{e-1}-1. \end{aligned}$$

Теперь используется классическое тождество

$$\sum_{d\mid m}\varphi(d)=m.$$

Значит, все направление \(p\) дает вклад

$$m\bigl(p^e-p^{e-1}-1\bigr)=n-\frac{n}{p}-\frac{n}{p^e}.$$

После суммирования по всем точным простым степеням числа \(n\) получаем формулу, на которой и держится решение:

$$\boxed{T(n)=\sum_{p^e\parallel n}\left(n-\frac{n}{p}-\frac{n}{p^e}\right).}$$

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

Положительные делители числа \(12\) таковы: \(1,2,3,4,6,12\). Ребра, у которых отношение концов является простым числом, имеют вид

$$1\to2,\quad 1\to3,\quad 2\to4,\quad 2\to6,\quad 3\to6,\quad 4\to12,\quad 6\to12.$$

Их веса равны

$$0,\ 1,\ 1,\ 1,\ 0,\ 2,\ 2,$$

потому что

$$\varphi(1)=1,\ \varphi(2)=1,\ \varphi(3)=2,\ \varphi(4)=2,\ \varphi(6)=2,\ \varphi(12)=4.$$

Следовательно, \(T(12)=7\).

Закрытая формула дает тот же ответ. Поскольку \(12=2^2\cdot 3\), имеем

$$T(12)=\left(12-\frac{12}{2}-\frac{12}{2^2}\right)+\left(12-\frac{12}{3}-\frac{12}{3}\right)=3+4=7.$$

Этот пример хорошо показывает смысл формулы: она заранее складывает суммарный вклад каждого простого направления в графе делителей.

Перестановка префиксной суммы по точным простым степеням

Теперь нужно вычислить

$$S(N)=\sum_{n\le N}T(n).$$

Для фиксированной степени \(p^e\) каждое подходящее \(n\) единственным образом записывается как

$$n=p^e m,\qquad p\nmid m.$$

Соответствующий вклад равен

$$n-\frac{n}{p}-\frac{n}{p^e}=m\bigl(p^e-p^{e-1}-1\bigr).$$

Поэтому

$$S(N)=\sum_{p^e\le N}\bigl(p^e-p^{e-1}-1\bigr)\sum_{\substack{m\le N/p^e\\ p\nmid m}} m.$$

Внутренняя сумма легко выражается через треугольные числа. Обозначим

$$\operatorname{Tri}(x)=\frac{x(x+1)}{2}.$$

Тогда

$$\sum_{\substack{m\le x\\ p\nmid m}} m=\operatorname{Tri}(x)-p\,\operatorname{Tri}\!\left(\left\lfloor\frac{x}{p}\right\rfloor\right),$$

поскольку исключаются именно кратные \(p,2p,\dots,p\lfloor x/p\rfloor\).

В итоге вся префиксная сумма принимает вид

$$\boxed{S(N)=\sum_{p^e\le N}\bigl(p^e-p^{e-1}-1\bigr)\left[\operatorname{Tri}\!\left(\left\lfloor\frac{N}{p^e}\right\rfloor\right)-p\,\operatorname{Tri}\!\left(\left\lfloor\frac{N}{p^{e+1}}\right\rfloor\right)\right].}$$

Почему простые числа больше \(\sqrt N\) удобно обрабатывать блоками

В реализации сумма естественно разбивается по границе \(\sqrt N\).

Если \(p\le\sqrt N\), то можно напрямую перебирать \(p,p^2,p^3,\dots\).

Если \(p>\sqrt N\), то уже \(p^2>N\), а значит возможен только показатель \(e=1\). Тогда

$$x=\left\lfloor\frac{N}{p}\right\rfloor<p,$$

то есть среди чисел от \(1\) до \(x\) нет ни одного положительного кратного \(p\). Следовательно,

$$\sum_{\substack{m\le x\\ p\nmid m}}m=\operatorname{Tri}(x),$$

и вклад большого простого числа упрощается до

$$\bigl(p-2\bigr)\operatorname{Tri}\!\left(\left\lfloor\frac{N}{p}\right\rfloor\right).$$

Дальше большие простые числа группируются по постоянному частному

$$q=\left\lfloor\frac{N}{p}\right\rfloor.$$

Все простые \(p\) из интервала

$$L=\left\lfloor\frac{N}{q+1}\right\rfloor+1,\qquad R=\left\lfloor\frac{N}{q}\right\rfloor$$

дают одно и то же значение \(q\). Поэтому их общий вклад равен

$$\operatorname{Tri}(q)\sum_{L\le p\le R}(p-2).$$

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

$$\pi(x)=\#\{p\le x\},\qquad P(x)=\sum_{p\le x}p.$$

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

$$\operatorname{Tri}(q)\Bigl(P(R)-P(L-1)-2\bigl(\pi(R)-\pi(L-1)\bigr)\Bigr).$$

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

Сначала проверяется сама формула

Эталонная реализация не просто подставляет закрытую формулу. Сначала она для малых \(n\) действительно строит граф делителей, считает сумму по ребрам из исходного определения и проверяет, что полученный результат совпадает с формулой по простым степеням на небольшом диапазоне. Затем подтверждаются несколько префиксных значений, в том числе \(S(10)=26\) и \(S(100)=5282\), и только после этого запускается вычисление для \(N=10^{12}\).

Сжатая предварительная обработка префиксов по простым

Алгоритму нужно много раз узнавать число простых \(\le x\) и сумму всех простых \(\le x\). Вместо хранения таблицы до самого \(N\) он использует только различные значения

$$\left\lfloor\frac{N}{1}\right\rfloor,\left\lfloor\frac{N}{2}\right\rfloor,\left\lfloor\frac{N}{3}\right\rfloor,\dots$$

Таких значений всего \(O(\sqrt N)\). На этом сжатом наборе сначала инициализируются количества и суммы всех целых чисел не меньше \(2\), а затем простое за простым вычитаются составные вклады. После подготовки значения \(\pi(x)\) и \(P(x)\) для всех нужных \(x\) извлекаются за константное время.

Итоговое накопление состоит из двух частей

В первой части перебираются все простые \(p\le\sqrt N\) вместе с их степенями \(p,p^2,p^3,\dots\), и добавляется точный вклад соответствующих простых степеней. Во второй части обрабатываются простые \(p>\sqrt N\), но уже не по одному, а блоками с одинаковым частным \(q=\lfloor N/p\rfloor\), используя подсчитанные префиксы количества простых и суммы простых на интервалах.

Вся арифметика ведется по модулю \(715827883\). В реализации на C++ промежуточные выражения удерживаются в 128-битном типе до взятия остатка, реализация на Python использует целые произвольной длины, а Java-вход не переписывает алгоритм отдельно, а вызывает скомпилированный C++-решатель и возвращает его ответ.

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

Множество различных значений \(\lfloor N/t\rfloor\) имеет размер \(O(\sqrt N)\). Именно на этой компрессии построены и структура запросов, и группировка больших простых. Поэтому потребление памяти равно \(O(\sqrt N)\).

Время в основном уходит на сжатое сито для префиксов по простым и на проход по простым числам до \(\sqrt N\) вместе с их степенями. Главный вывод здесь не в тонкой формуле, а в масштабе: алгоритм существенно лучше линейного по \(N\) и полностью избегает попытки строить граф делителей для всех \(n\le 10^{12}\). Именно это делает вычисление выполнимым.

Сноски и ссылки

  1. Страница задачи: https://projecteuler.net/problem=931
  2. Функция Эйлера: Wikipedia - Euler's totient function
  3. Делитель: Wikipedia - Divisor
  4. Функция подсчета простых: Wikipedia - Prime-counting function
  5. Метод гиперболы Дирихле: Wikipedia - Dirichlet hyperbola method

ملخص المسألة

لكل عدد صحيح موجب \(n\)، ننظر إلى رسم بياني رؤوسه هي جميع القواسم الموجبة لـ \(n\). توجد حافة صاعدة من \(b\) إلى \(a\) إذا وفقط إذا كان \(a/b\) عددًا أوليًا. أي إن كل حافة تمثل ضرب قاسم ما في عامل أولي واحد فقط. الكمية المرتبطة بـ \(n\) هي مجموع أوزان هذه الحواف:

$$T(n)=\sum_{\substack{b\mid n,\ a\mid n\\ a/b\text{ prime}}}\bigl(\varphi(a)-\varphi(b)\bigr).$$

والمطلوب في النهاية هو حساب المجموع التراكمي

$$S(N)=\sum_{n=1}^{N} T(n)\pmod{715827883},\qquad N=10^{12}.$$

الطريقة المباشرة كانت ستتطلب بناء رسم القواسم لكل \(n\le N\) على حدة، وهذا غير عملي تمامًا. الحل ينجح لأن مجموع الحواف ينهار أولًا إلى صيغة مغلقة تعتمد على القوى الأولية الدقيقة، ثم يعاد ترتيب الجمع الكلي بحيث يتحول إلى حلقات على القوى الأولية مع بنية سريعة لاستعلامات البادئات الخاصة بالأوليات.

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

الفكرة الحاسمة هي أن نفهم \(T(n)\) تمامًا لعدد ثابت \(n\)، ثم بعد ذلك فقط نجمع على جميع \(n\le N\).

رسم القواسم ينفصل بحسب اتجاه كل عدد أولي

لنكتب

$$n=\prod_i p_i^{e_i}.$$

كل حافة مسموحة في الرسم تنشأ باختيار عدد أولي \(p\mid n\) والانتقال من قاسم \(d\) إلى \(pd\)، بشرط أن يبقى \(pd\) قاسمًا لـ \(n\). لذلك يمكن تقسيم الرسم كله إلى عائلات مستقلة من الحواف، عائلة لكل أولي يقسم \(n\).

ثبّت الآن عددًا أوليًا \(p\) بحيث \(p^e\parallel n\)، واكتب

$$n=p^e m,\qquad p\nmid m.$$

عندها تكون كل حافة نسبتها \(p\) من الصورة

$$p^r d \longrightarrow p^{r+1} d,$$

حيث \(0\le r<e\) و\(d\mid m\). إذن الحصول على المساهمة الكاملة لاتجاه \(p\) لا يحتاج إلا إلى جمع فروق دالة أويلر على هذه الحواف.

حساب مساهمة اتجاه أولي واحد بدقة

بما أن \(p\nmid d\)، فإن دالة أويلر تحقق

$$\varphi(p^r d)=\varphi(p^r)\varphi(d).$$

وهنا تظهر حالتان مختلفتان.

إذا كان \(r=0\)، فإن

$$\varphi(pd)-\varphi(d)=(p-1)\varphi(d)-\varphi(d)=(p-2)\varphi(d).$$

أما إذا كان \(r\ge 1\)، فإن

$$\varphi(p^{r+1}d)-\varphi(p^r d)=p^r(p-1)\varphi(d)-p^{r-1}(p-1)\varphi(d)=p^{r-1}(p-1)^2\varphi(d).$$

إذًا تصبح المساهمة الكلية لجميع حواف الاتجاه \(p\)

$$\sum_{d\mid m}\varphi(d)\left[(p-2)+\sum_{r=1}^{e-1}p^{r-1}(p-1)^2\right].$$

والقوس يبسط إلى

$$\begin{aligned} (p-2)+(p-1)^2\sum_{r=1}^{e-1}p^{r-1} &=(p-2)+(p-1)^2\frac{p^{e-1}-1}{p-1} \\ &=(p-2)+(p-1)(p^{e-1}-1) \\ &=p^e-p^{e-1}-1. \end{aligned}$$

والآن نستخدم الهوية الكلاسيكية

$$\sum_{d\mid m}\varphi(d)=m.$$

ومن ثم فإن الاتجاه الكامل المرتبط بـ \(p\) يساهم بمقدار

$$m\bigl(p^e-p^{e-1}-1\bigr)=n-\frac{n}{p}-\frac{n}{p^e}.$$

وبالجمع على جميع القوى الأولية الدقيقة في \(n\) نحصل على الصيغة المغلقة التي تعتمد عليها التطبيقات:

$$\boxed{T(n)=\sum_{p^e\parallel n}\left(n-\frac{n}{p}-\frac{n}{p^e}\right).}$$

مثال محلول: \(n=12\)

قواسم \(12\) هي \(1,2,3,4,6,12\). والحواف التي تكون نسبتها عددًا أوليًا هي

$$1\to2,\quad 1\to3,\quad 2\to4,\quad 2\to6,\quad 3\to6,\quad 4\to12,\quad 6\to12.$$

وأوزان هذه الحواف هي

$$0,\ 1,\ 1,\ 1,\ 0,\ 2,\ 2,$$

لأن

$$\varphi(1)=1,\ \varphi(2)=1,\ \varphi(3)=2,\ \varphi(4)=2,\ \varphi(6)=2,\ \varphi(12)=4.$$

ومن ثم \(T(12)=7\).

وتعطي الصيغة المغلقة النتيجة نفسها تمامًا. بما أن \(12=2^2\cdot 3\)، فإن

$$T(12)=\left(12-\frac{12}{2}-\frac{12}{2^2}\right)+\left(12-\frac{12}{3}-\frac{12}{3}\right)=3+4=7.$$

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

إعادة ترتيب المجموع التراكمي بحسب القوى الأولية الدقيقة

لحساب

$$S(N)=\sum_{n\le N}T(n),$$

نكتب كل عدد \(n\) يساهم في الحد المرتبط بـ \(p^e\) على الصورة الوحيدة

$$n=p^e m,\qquad p\nmid m.$$

وعندئذ يكون الحد المقابل مساويًا لـ

$$n-\frac{n}{p}-\frac{n}{p^e}=m\bigl(p^e-p^{e-1}-1\bigr).$$

ومن ثم نحصل على

$$S(N)=\sum_{p^e\le N}\bigl(p^e-p^{e-1}-1\bigr)\sum_{\substack{m\le N/p^e\\ p\nmid m}} m.$$

المجموع الداخلي هو ببساطة مجموع الأعداد غير القابلة للقسمة على \(p\). إذا عرفنا

$$\operatorname{Tri}(x)=\frac{x(x+1)}{2},$$

فإن

$$\sum_{\substack{m\le x\\ p\nmid m}} m=\operatorname{Tri}(x)-p\,\operatorname{Tri}\!\left(\left\lfloor\frac{x}{p}\right\rfloor\right),$$

لأن الحدود المستبعدة هي تحديدًا \(p,2p,\dots,p\lfloor x/p\rfloor\). وهكذا يصبح المجموع الكامل

$$\boxed{S(N)=\sum_{p^e\le N}\bigl(p^e-p^{e-1}-1\bigr)\left[\operatorname{Tri}\!\left(\left\lfloor\frac{N}{p^e}\right\rfloor\right)-p\,\operatorname{Tri}\!\left(\left\lfloor\frac{N}{p^{e+1}}\right\rfloor\right)\right].}$$

لماذا يمكن تجميع الأوليات الأكبر من \(\sqrt N\) في كتل؟

التنفيذ يقسم الجمع طبيعيًا عند \(\sqrt N\).

إذا كان \(p\le\sqrt N\)، فمن المناسب المرور مباشرة على \(p,p^2,p^3,\dots\).

أما إذا كان \(p>\sqrt N\)، فإن \(p^2>N\)، أي إن الأس الممكن الوحيد هو \(e=1\). وفي هذه الحالة

$$x=\left\lfloor\frac{N}{p}\right\rfloor<p,$$

وبالتالي لا يوجد أي مضاعف موجب لـ \(p\) داخل المجال \(\{1,\dots,x\}\). لذلك

$$\sum_{\substack{m\le x\\ p\nmid m}}m=\operatorname{Tri}(x),$$

فتتبسط مساهمة كل أولي كبير إلى

$$\bigl(p-2\bigr)\operatorname{Tri}\!\left(\left\lfloor\frac{N}{p}\right\rfloor\right).$$

بعد ذلك تجمع الأوليات الكبيرة بحسب خارج القسمة الثابت

$$q=\left\lfloor\frac{N}{p}\right\rfloor.$$

فكل الأوليات داخل المجال

$$L=\left\lfloor\frac{N}{q+1}\right\rfloor+1,\qquad R=\left\lfloor\frac{N}{q}\right\rfloor$$

تعطي القيمة نفسها لـ \(q\). ولذلك يكون مجموع مساهمتها

$$\operatorname{Tri}(q)\sum_{L\le p\le R}(p-2).$$

وهكذا تختزل جهة الأوليات الكبيرة كلها إلى استعلامات على دالتين بادئيتين:

$$\pi(x)=\#\{p\le x\},\qquad P(x)=\sum_{p\le x}p.$$

ومتى أصبحت هاتان الدالتان متاحتين بسرعة، يمكن حساب كل كتلة على الصورة

$$\operatorname{Tri}(q)\Bigl(P(R)-P(L-1)-2\bigl(\pi(R)-\pi(L-1)\bigr)\Bigr).$$

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

التحقق من الهوية قبل الحساب الكبير

التنفيذ المرجعي لا يفترض صحة الصيغة من دون اختبار. فهو يبني رسم القواسم فعلًا للأعداد الصغيرة، ويحسب مجموع الحواف مباشرة من التعريف، ثم يتأكد من تطابق ذلك مع صيغة القوى الأولية على مجال تحقق صغير. وبعد ذلك يفحص عدة قيم تراكمية، منها \(S(10)=26\) و\(S(100)=5282\)، قبل الانتقال إلى حساب \(N=10^{12}\).

تمهيد مضغوط لبادئات الأوليات

الخوارزمية تحتاج إلى عدد كبير من الاستعلامات عن عدد الأوليات حتى \(x\) وعن مجموعها حتى \(x\). بدل الاحتفاظ بجدول كامل حتى \(N\)، فإنها تخزن فقط القيم المختلفة لـ

$$\left\lfloor\frac{N}{1}\right\rfloor,\left\lfloor\frac{N}{2}\right\rfloor,\left\lfloor\frac{N}{3}\right\rfloor,\dots$$

وهذه القيم المختلفة لا يتجاوز عددها \(O(\sqrt N)\). وعلى هذا الجهاز المضغوط تُهيَّأ أولًا أعداد ومجاميع جميع الأعداد الصحيحة التي لا تقل عن \(2\)، ثم تُزال الإسهامات المركبة أوليًا بعد أولي. وبعد انتهاء هذه المرحلة تصبح قيم \(\pi(x)\) و\(P(x)\) متاحة بزمن ثابت لكل \(x\) مهم في الحساب.

تجميع الجواب في مرحلتين

في المرحلة الأولى تمر الخوارزمية على جميع الأوليات \(p\le\sqrt N\) ومعها القوى \(p,p^2,p^3,\dots\)، وتضيف حد القوة الأولية الدقيق. وفي المرحلة الثانية تعالج الأوليات \(p>\sqrt N\) لا واحدًا واحدًا، بل على شكل كتل ذات خارج قسمة ثابت \(q=\lfloor N/p\rfloor\)، مستعملةً عدد الأوليات ومجموعها على الفترات.

كل الحسابات تؤخذ بترديد \(715827883\). تنفيذ C++ يستخدم أعدادًا وسيطة من 128 بت قبل أخذ الباقي، وتنفيذ Python يعتمد على الأعداد الصحيحة ذات الدقة غير المحدودة، أما مدخل Java فلا يعيد كتابة الخوارزمية كاملة، بل يستدعي المحلل المترجم بلغة C++ ثم يعيد النتيجة التي ينتجها.

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

مجموعة القيم المختلفة لـ \(\lfloor N/t\rfloor\) تملك حجمًا \(O(\sqrt N)\). وكل من بنية الاستعلامات وتجميع الأوليات الكبيرة مبني على هذا الضغط. لذلك يكون استهلاك الذاكرة \(O(\sqrt N)\).

أما الزمن فتسيطر عليه عملية الغربلة المضغوطة لبادئات الأوليات، مع المرور على الأوليات حتى \(\sqrt N\) وقواها. المهم هنا ليس ثابتًا دقيقًا بقدر ما هو الترتيب العام: الخوارزمية أفضل بكثير من الحلول الخطية في \(N\)، ولا تحاول أصلًا إنشاء رسم قواسم لكل \(n\le 10^{12}\). وهذا بالتحديد ما يجعل المسألة قابلة للحساب.

الحواشي والمراجع

  1. صفحة المسألة: https://projecteuler.net/problem=931
  2. دالة أويلر: Wikipedia - Euler's totient function
  3. القاسم أو القواسم: Wikipedia - Divisor
  4. دالة عدّ الأوليات: Wikipedia - Prime-counting function
  5. طريقة القطع الزائد لديريشليه: Wikipedia - Dirichlet hyperbola method