For each integer \(M \ge 2\), define
$$R(M)=\sum_{\substack{1\le p \lt q \le M\\ p+q \ge M\\ \gcd(p,q)=1}}\frac{1}{pq}, \qquad S(N)=\sum_{M=2}^{N}R(M).$$
The checkpoints are \(S(2)=\tfrac12\), \(S(10)=6.9146825397\ldots\), and \(S(100)=58.2962380621\ldots\). The task is to evaluate \(S(10^7)\) to four decimal places. A direct evaluation over all triples \((M,p,q)\) is far too slow, so the implementation reorganizes the sum and replaces the coprimality test by Möbius inclusion-exclusion.
Fix a coprime pair \(1 \le p \lt q \le N\). Such a pair contributes to \(R(M)\) exactly when \(q \le M\) and \(M \le p+q\). Since \(M\) is also bounded by \(N\), the admissible values are
$$q \le M \le \min(N,p+q).$$
Therefore this pair appears with multiplicity
$$w_N(p,q)=\min(N,p+q)-q+1.$$
After switching the order of summation,
$$S(N)=\sum_{\substack{1\le p \lt q \le N\\ \gcd(p,q)=1}}\frac{w_N(p,q)}{pq} =\sum_{q=2}^{N}\frac{1}{q}\sum_{\substack{1\le p \lt q\\ \gcd(p,q)=1}}\frac{w_N(p,q)}{p}.$$
This is the key simplification: the problem becomes a sum of independent contributions indexed by the larger denominator \(q\).
Define the coprime reciprocal sum and the coprime counting function
$$A_q(m)=\sum_{\substack{1\le p \le m\\ \gcd(p,q)=1}}\frac{1}{p}, \qquad C_q(m)=\sum_{\substack{1\le p \le m\\ \gcd(p,q)=1}}1.$$
Using the identity
$$\mathbf{1}_{\gcd(p,q)=1}=\sum_{d \mid \gcd(p,q)}\mu(d),$$
we can rewrite both quantities in terms of divisors of \(q\). If \(H_t=\sum_{j=1}^{t}\frac{1}{j}\) denotes the \(t\)-th harmonic number, then
$$A_q(m)=\sum_{d \mid q}\frac{\mu(d)}{d}H_{\lfloor m/d \rfloor}, \qquad C_q(m)=\sum_{d \mid q}\mu(d)\left\lfloor \frac{m}{d}\right\rfloor.$$
Only squarefree divisors matter, because \(\mu(d)=0\) whenever \(d\) contains a repeated prime factor. That is why the implementation factors \(q\) into distinct primes and enumerates all squarefree divisor subsets with their Möbius signs.
If \(q \le \lfloor N/2 \rfloor\), then for every \(p \lt q\) we have
$$p+q \le (q-1)+q = 2q-1 \le N-1,$$
so the upper bound \(\min(N,p+q)\) is just \(p+q\). Hence
$$w_N(p,q)=p+1.$$
The contribution of this \(q\) becomes
$$\frac{1}{q}\sum_{\substack{1\le p \lt q\\ \gcd(p,q)=1}}\frac{p+1}{p} =\frac{1}{q}\left(\sum_{\substack{1\le p \lt q\\ \gcd(p,q)=1}}1+\sum_{\substack{1\le p \lt q\\ \gcd(p,q)=1}}\frac{1}{p}\right).$$
The first sum is \(\varphi(q)\), because among \(1,2,\dots,q-1\) there are exactly \(\varphi(q)\) integers coprime to \(q\). Therefore
$$T_N(q)=\frac{\varphi(q)+A_q(q-1)}{q}\qquad \text{for } q \le \left\lfloor \frac{N}{2}\right\rfloor.$$
Now let \(q > \lfloor N/2 \rfloor\) and set \(a=N-q\). Then \(0 \le a \lt q\), and the weight changes at \(p=a\):
$$w_N(p,q)= \begin{cases} p+1, & 1 \le p \le a,\\ a+1, & a \lt p \lt q. \end{cases}$$
So the contribution splits into two pieces:
$$T_N(q)=\frac{1}{q}\left(\sum_{\substack{1\le p \le a\\ \gcd(p,q)=1}}\left(1+\frac{1}{p}\right) +(a+1)\sum_{\substack{a \lt p \lt q\\ \gcd(p,q)=1}}\frac{1}{p}\right).$$
Writing the second reciprocal sum as \(A_q(q-1)-A_q(a)\), we obtain
$$T_N(q)=\frac{C_q(a)+A_q(a)+(a+1)\bigl(A_q(q-1)-A_q(a)\bigr)}{q}.$$
This is the exact branch used by the implementation once \(q\) passes the midpoint.
Combining the two ranges gives the complete formula
$$\boxed{ S(N)= \sum_{q=2}^{\lfloor N/2 \rfloor}\frac{\varphi(q)+A_q(q-1)}{q} + \sum_{q=\lfloor N/2 \rfloor+1}^{N} \frac{C_q(N-q)+A_q(N-q)+(N-q+1)\bigl(A_q(q-1)-A_q(N-q)\bigr)}{q} .}$$
All expensive work has now been reduced to three precomputable ingredients: Euler's totient values, harmonic prefixes, and factorizations of the integers \(2,\dots,N\).
For \(q=5\), we are still in the first branch. Since \(\varphi(5)=4\) and
$$A_5(4)=1+\frac12+\frac13+\frac14=\frac{25}{12},$$
we get
$$T_{10}(5)=\frac{4+\frac{25}{12}}{5}=\frac{73}{60}=1.216666\ldots$$
For \(q=8\), we are in the second branch with \(a=10-8=2\). Here
$$C_8(2)=1,\qquad A_8(2)=1,\qquad A_8(7)-A_8(2)=\frac13+\frac15+\frac17=\frac{71}{105},$$
so
$$T_{10}(8)=\frac{1+1+3\cdot \frac{71}{105}}{8}=\frac{141}{280}=0.503571\ldots$$
Summing all contributions from \(q=2\) to \(10\) yields
$$S(10)=\frac{3485}{504}=6.914682539682\ldots,$$
which matches the published checkpoint.
The C++, Python, and Java implementations all follow the same pipeline. First they build a smallest-prime-factor table up to \(N\), which makes factorization of each \(q\) fast and also supports a linear-time computation of \(\varphi(q)\). Next they precompute the harmonic prefix array \(H_0,H_1,\dots,H_N\). Then each value of \(q\) is processed independently: factor \(q\), enumerate the squarefree divisors induced by its distinct prime factors, evaluate the Möbius formulas for \(A_q\) and \(C_q\), apply the appropriate branch, and add the result to the global sum.
The C++ and Java implementations additionally use compensated summation to reduce floating-point drift, and the C++ version can split the outer loop into parallel chunks. None of those engineering choices changes the mathematics; they only improve numerical stability and runtime for the target size.
Building the smallest-prime-factor table, the totient table, and the harmonic prefix array costs \(O(N)\) time and \(O(N)\) memory. For each \(q\), the factorization is recovered from the precomputed table, and the Möbius expansion enumerates exactly \(2^{\omega(q)}\) squarefree divisors, where \(\omega(q)\) is the number of distinct prime factors of \(q\). Therefore the total running time is
$$O\!\left(N+\sum_{q=2}^{N}2^{\omega(q)}\right),$$
with \(O(N)\) memory usage. For the specific target \(N=10^7\), every integer has at most eight distinct prime factors, so the subset enumeration remains very small in practice. Parallel execution reduces wall-clock time but does not change the asymptotic work.
Für jedes \(M \ge 2\) sei
$$R(M)=\sum_{\substack{1\le p \lt q \le M\\ p+q \ge M\\ \gcd(p,q)=1}}\frac{1}{pq}, \qquad S(N)=\sum_{M=2}^{N}R(M).$$
Die Kontrollwerte lauten \(S(2)=\tfrac12\), \(S(10)=6.9146825397\ldots\) und \(S(100)=58.2962380621\ldots\). Gesucht ist \(S(10^7)\) auf vier Nachkommastellen. Eine direkte Auswertung aller Tripel \((M,p,q)\) ist viel zu teuer; deshalb ordnet die Lösung die Summe um und ersetzt die Teilerfremdheitsbedingung durch Möbius-Inklusion-Exklusion.
Fixiere ein teilerfremdes Paar \(1 \le p \lt q \le N\). Dieses Paar trägt genau dann zu \(R(M)\) bei, wenn \(q \le M\) und zugleich \(M \le p+q\) gilt. Wegen der zusätzlichen Schranke \(M \le N\) sind die zulässigen Werte
$$q \le M \le \min(N,p+q).$$
Die Multiplizität dieses Paares ist also
$$w_N(p,q)=\min(N,p+q)-q+1.$$
Nach dem Vertauschen der Summen erhält man
$$S(N)=\sum_{\substack{1\le p \lt q \le N\\ \gcd(p,q)=1}}\frac{w_N(p,q)}{pq} =\sum_{q=2}^{N}\frac{1}{q}\sum_{\substack{1\le p \lt q\\ \gcd(p,q)=1}}\frac{w_N(p,q)}{p}.$$
Damit zerfällt das gesamte Problem in unabhängige Beiträge, die nur noch vom größeren Nenner \(q\) abhängen.
Definiere die teilerfremde Kehrwertsumme und die zugehörige Zählfunktion durch
$$A_q(m)=\sum_{\substack{1\le p \le m\\ \gcd(p,q)=1}}\frac{1}{p}, \qquad C_q(m)=\sum_{\substack{1\le p \le m\\ \gcd(p,q)=1}}1.$$
Mit der Identität
$$\mathbf{1}_{\gcd(p,q)=1}=\sum_{d \mid \gcd(p,q)}\mu(d)$$
lassen sich beide Größen in Divisorsummen über \(q\) umschreiben. Bezeichnet \(H_t=\sum_{j=1}^{t}\frac{1}{j}\) die \(t\)-te harmonische Zahl, dann gilt
$$A_q(m)=\sum_{d \mid q}\frac{\mu(d)}{d}H_{\lfloor m/d \rfloor}, \qquad C_q(m)=\sum_{d \mid q}\mu(d)\left\lfloor \frac{m}{d}\right\rfloor.$$
Nur quadratfreie Teiler tragen bei, weil \(\mu(d)=0\) für nicht quadratfreie \(d\) ist. Genau deshalb zerlegt die Implementierung jedes \(q\) in seine verschiedenen Primteiler und erzeugt daraus alle quadratfreien Teilermengen mit dem passenden Möbius-Vorzeichen.
Ist \(q \le \lfloor N/2 \rfloor\), so gilt für jedes \(p \lt q\)
$$p+q \le (q-1)+q = 2q-1 \le N-1,$$
also ist \(\min(N,p+q)=p+q\). Damit folgt
$$w_N(p,q)=p+1.$$
Der Beitrag dieses \(q\) ist daher
$$\frac{1}{q}\sum_{\substack{1\le p \lt q\\ \gcd(p,q)=1}}\frac{p+1}{p} =\frac{1}{q}\left(\sum_{\substack{1\le p \lt q\\ \gcd(p,q)=1}}1+\sum_{\substack{1\le p \lt q\\ \gcd(p,q)=1}}\frac{1}{p}\right).$$
Die erste Summe ist \(\varphi(q)\), denn unter \(1,2,\dots,q-1\) gibt es genau \(\varphi(q)\) Zahlen, die zu \(q\) teilerfremd sind. Somit
$$T_N(q)=\frac{\varphi(q)+A_q(q-1)}{q}\qquad \text{für } q \le \left\lfloor \frac{N}{2}\right\rfloor.$$
Nun sei \(q > \lfloor N/2 \rfloor\) und \(a=N-q\). Dann gilt \(0 \le a \lt q\), und das Gewicht ändert sich bei \(p=a\):
$$w_N(p,q)= \begin{cases} p+1, & 1 \le p \le a,\\ a+1, & a \lt p \lt q. \end{cases}$$
Also zerfällt der Beitrag in zwei Teile:
$$T_N(q)=\frac{1}{q}\left(\sum_{\substack{1\le p \le a\\ \gcd(p,q)=1}}\left(1+\frac{1}{p}\right) +(a+1)\sum_{\substack{a \lt p \lt q\\ \gcd(p,q)=1}}\frac{1}{p}\right).$$
Schreibt man die zweite Kehrwertsumme als \(A_q(q-1)-A_q(a)\), dann erhält man
$$T_N(q)=\frac{C_q(a)+A_q(a)+(a+1)\bigl(A_q(q-1)-A_q(a)\bigr)}{q}.$$
Genau diese Verzweigung verwendet die Implementierung für den Bereich oberhalb der Mitte.
Die beiden Bereiche zusammen liefern
$$\boxed{ S(N)= \sum_{q=2}^{\lfloor N/2 \rfloor}\frac{\varphi(q)+A_q(q-1)}{q} + \sum_{q=\lfloor N/2 \rfloor+1}^{N} \frac{C_q(N-q)+A_q(N-q)+(N-q+1)\bigl(A_q(q-1)-A_q(N-q)\bigr)}{q} .}$$
Damit ist die ursprüngliche Dreifachsummation auf drei vorab berechenbare Bausteine reduziert: Totienten, harmonische Präfixsummen und Faktorisierungen der Zahlen \(2,\dots,N\).
Für \(q=5\) liegt man noch im ersten Fall. Es ist \(\varphi(5)=4\) und
$$A_5(4)=1+\frac12+\frac13+\frac14=\frac{25}{12},$$
daher
$$T_{10}(5)=\frac{4+\frac{25}{12}}{5}=\frac{73}{60}=1.216666\ldots$$
Für \(q=8\) gilt dagegen \(a=10-8=2\). Dann ist
$$C_8(2)=1,\qquad A_8(2)=1,\qquad A_8(7)-A_8(2)=\frac13+\frac15+\frac17=\frac{71}{105},$$
also
$$T_{10}(8)=\frac{1+1+3\cdot \frac{71}{105}}{8}=\frac{141}{280}=0.503571\ldots$$
Die Summe aller Beiträge für \(q=2,\dots,10\) ist
$$S(10)=\frac{3485}{504}=6.914682539682\ldots,$$
genau der angegebene Kontrollwert.
Die C++-, Python- und Java-Implementierungen benutzen dieselbe mathematische Pipeline. Zuerst wird bis \(N\) eine Tabelle des kleinsten Primteilers aufgebaut; damit lässt sich jedes \(q\) schnell faktorisieren, und zugleich kann \(\varphi(q)\) in linearer Zeit mitberechnet werden. Danach wird das harmonische Präfix \(H_0,H_1,\dots,H_N\) erzeugt. Anschließend wird jedes \(q\) unabhängig verarbeitet: \(q\) faktorisieren, aus den verschiedenen Primteilern alle quadratfreien Teiler mit Möbius-Vorzeichen erzeugen, daraus \(A_q\) und \(C_q\) auswerten, den passenden Fall anwenden und den Beitrag zur Gesamtsumme addieren.
Die C++- und Java-Implementierungen verwenden zusätzlich kompensierte Summation, um Rundungsfehler klein zu halten; die C++-Variante kann die äußere Schleife außerdem in parallele Blöcke aufteilen. Diese technischen Details ändern die Formel nicht, sondern verbessern nur Stabilität und Laufzeit für die Zielgröße.
Der Aufbau der Tabelle des kleinsten Primteilers, der Totiententabelle und des harmonischen Präfixarrays kostet \(O(N)\) Zeit und \(O(N)\) Speicher. Für jedes \(q\) wird die Faktorisierung aus der vorberechneten Tabelle gewonnen, und die Möbius-Expansion durchläuft genau \(2^{\omega(q)}\) quadratfreie Teiler, wobei \(\omega(q)\) die Anzahl verschiedener Primteiler von \(q\) bezeichnet. Insgesamt ergibt sich damit
$$O\!\left(N+\sum_{q=2}^{N}2^{\omega(q)}\right)$$
Zeit bei \(O(N)\) Speicher. Für das konkrete Ziel \(N=10^7\) besitzt jede Zahl höchstens acht verschiedene Primteiler, daher bleibt die Teilermengenerzeugung in der Praxis sehr klein. Parallelisierung verkürzt nur die Laufzeit an der Uhr, nicht aber die asymptotische Arbeit.
Her \(M \ge 2\) için
$$R(M)=\sum_{\substack{1\le p \lt q \le M\\ p+q \ge M\\ \gcd(p,q)=1}}\frac{1}{pq}, \qquad S(N)=\sum_{M=2}^{N}R(M).$$
Kontrol değerleri \(S(2)=\tfrac12\), \(S(10)=6.9146825397\ldots\) ve \(S(100)=58.2962380621\ldots\) şeklindedir. Amaç \(S(10^7)\) değerini dört ondalık basamağa kadar hesaplamaktır. Tüm \((M,p,q)\) üçlülerini doğrudan dolaşmak çok pahalı olduğundan çözüm toplamın sırasını değiştirir ve aralarında asal olma koşulunu Möbius dahil etme-çıkarma ile yazar.
\(1 \le p \lt q \le N\) ve \(\gcd(p,q)=1\) olan sabit bir çift alalım. Bu çift \(R(M)\)'ye ancak \(q \le M\) ve \(M \le p+q\) olduğunda katkı verir. Ayrıca \(M \le N\) şartı da bulunduğu için uygun değerler
$$q \le M \le \min(N,p+q)$$
aralığındadır. Dolayısıyla bu çiftin çokluğu
$$w_N(p,q)=\min(N,p+q)-q+1$$
olur. Toplamların yerini değiştirince
$$S(N)=\sum_{\substack{1\le p \lt q \le N\\ \gcd(p,q)=1}}\frac{w_N(p,q)}{pq} =\sum_{q=2}^{N}\frac{1}{q}\sum_{\substack{1\le p \lt q\\ \gcd(p,q)=1}}\frac{w_N(p,q)}{p}$$
elde edilir. Böylece problem, büyük payda olan \(q\)'ya göre ayrılmış bağımsız katkılara dönüşür.
Şu iki matematiksel büyüklüğü tanımlayalım:
$$A_q(m)=\sum_{\substack{1\le p \le m\\ \gcd(p,q)=1}}\frac{1}{p}, \qquad C_q(m)=\sum_{\substack{1\le p \le m\\ \gcd(p,q)=1}}1.$$
Kimlik
$$\mathbf{1}_{\gcd(p,q)=1}=\sum_{d \mid \gcd(p,q)}\mu(d)$$
sayesinde bu ifadeleri \(q\)'nun bölenleri üzerinden yazabiliriz. Burada \(H_t=\sum_{j=1}^{t}\frac{1}{j}\) \(t\)'inci harmonik sayı olsun. O zaman
$$A_q(m)=\sum_{d \mid q}\frac{\mu(d)}{d}H_{\lfloor m/d \rfloor}, \qquad C_q(m)=\sum_{d \mid q}\mu(d)\left\lfloor \frac{m}{d}\right\rfloor$$
olur. \(\mu(d)=0\) yalnızca kareden bağımsız olmayan bölenlerde sıfır olduğu için, gerçekte sadece kareden bağımsız bölenler katkı verir. Bu nedenle uygulama her \(q\)'yu farklı asal çarpanlarına ayırır ve bu çarpanların tüm alt kümelerini kullanarak kareden bağımsız bölenleri ve işaretlerini üretir.
Eğer \(q \le \lfloor N/2 \rfloor\) ise her \(p \lt q\) için
$$p+q \le (q-1)+q = 2q-1 \le N-1$$
olduğundan \(\min(N,p+q)=p+q\) gelir. Böylece
$$w_N(p,q)=p+1$$
olur. Bu \(q\)'nun katkısı
$$\frac{1}{q}\sum_{\substack{1\le p \lt q\\ \gcd(p,q)=1}}\frac{p+1}{p} =\frac{1}{q}\left(\sum_{\substack{1\le p \lt q\\ \gcd(p,q)=1}}1+\sum_{\substack{1\le p \lt q\\ \gcd(p,q)=1}}\frac{1}{p}\right)$$
şeklindedir. İlk toplam \(\varphi(q)\)'dur; çünkü \(1,2,\dots,q-1\) arasında tam olarak \(\varphi(q)\) tane sayı \(q\) ile aralarında asaldır. Sonuç olarak
$$T_N(q)=\frac{\varphi(q)+A_q(q-1)}{q}\qquad \text{eğer } q \le \left\lfloor \frac{N}{2}\right\rfloor.$$
Şimdi \(q > \lfloor N/2 \rfloor\) olsun ve \(a=N-q\) tanımını yapalım. Bu durumda \(0 \le a \lt q\) ve ağırlık \(p=a\) noktasında değişir:
$$w_N(p,q)= \begin{cases} p+1, & 1 \le p \le a,\\ a+1, & a \lt p \lt q. \end{cases}$$
Buna göre katkı iki parçaya ayrılır:
$$T_N(q)=\frac{1}{q}\left(\sum_{\substack{1\le p \le a\\ \gcd(p,q)=1}}\left(1+\frac{1}{p}\right) +(a+1)\sum_{\substack{a \lt p \lt q\\ \gcd(p,q)=1}}\frac{1}{p}\right).$$
İkinci ters toplamı \(A_q(q-1)-A_q(a)\) diye yazarsak
$$T_N(q)=\frac{C_q(a)+A_q(a)+(a+1)\bigl(A_q(q-1)-A_q(a)\bigr)}{q}$$
elde edilir. Uygulamanın orta noktadan sonra kullandığı formül tam olarak budur.
İki bölgeyi birleştirince
$$\boxed{ S(N)= \sum_{q=2}^{\lfloor N/2 \rfloor}\frac{\varphi(q)+A_q(q-1)}{q} + \sum_{q=\lfloor N/2 \rfloor+1}^{N} \frac{C_q(N-q)+A_q(N-q)+(N-q+1)\bigl(A_q(q-1)-A_q(N-q)\bigr)}{q} .}$$
Böylece başlangıçtaki üç katmanlı toplam, önceden hazırlanabilen üç bileşene indirgenmiş olur: totient değerleri, harmonik önekler ve \(2,\dots,N\) sayılarının çarpanlara ayrılması.
\(q=5\) için ilk daldayız. \(\varphi(5)=4\) ve
$$A_5(4)=1+\frac12+\frac13+\frac14=\frac{25}{12}$$
olduğundan
$$T_{10}(5)=\frac{4+\frac{25}{12}}{5}=\frac{73}{60}=1.216666\ldots$$
\(q=8\) için ise ikinci dala geçilir ve \(a=10-8=2\) olur. Burada
$$C_8(2)=1,\qquad A_8(2)=1,\qquad A_8(7)-A_8(2)=\frac13+\frac15+\frac17=\frac{71}{105}$$
olduğu için
$$T_{10}(8)=\frac{1+1+3\cdot \frac{71}{105}}{8}=\frac{141}{280}=0.503571\ldots$$
Tüm \(q=2,\dots,10\) katkılarını toplarsak
$$S(10)=\frac{3485}{504}=6.914682539682\ldots$$
çıkar; bu da verilen kontrol değeriyle aynıdır.
C++, Python ve Java implementasyonları aynı matematiksel akışı izler. Önce \(N\)'ye kadar en küçük asal çarpan tablosu kurulur; bu tablo sayesinde her \(q\) hızlıca çarpanlara ayrılır ve aynı anda \(\varphi(q)\) değerleri de doğrusal sürede hesaplanır. Sonra harmonik önek dizisi \(H_0,H_1,\dots,H_N\) oluşturulur. Ardından her \(q\) bağımsız işlenir: \(q\)'nun farklı asal çarpanlarını bul, bunlardan kareden bağımsız bölenleri ve Möbius işaretlerini üret, \(A_q\) ve \(C_q\) formüllerini değerlendir, uygun dalı uygula ve sonucu toplam cevaba ekle.
C++ ve Java implementasyonları kayan nokta hatasını azaltmak için telafili toplama da kullanır; C++ sürümü dış döngüyü paralel parçalara da bölebilir. Bunlar matematiği değiştirmez, sadece hedef boyutta daha sağlam ve hızlı çalışmayı sağlar.
En küçük asal çarpan tablosu, totient dizisi ve harmonik önek dizisinin hazırlanması \(O(N)\) zaman ve \(O(N)\) bellek gerektirir. Her \(q\) için çarpanlara ayırma işlemi önceden kurulan tablodan alınır ve Möbius açılımı tam olarak \(2^{\omega(q)}\) kareden bağımsız bölen dolaşır; burada \(\omega(q)\), \(q\)'nun farklı asal çarpan sayısıdır. Bu yüzden toplam çalışma süresi
$$O\!\left(N+\sum_{q=2}^{N}2^{\omega(q)}\right)$$
şeklindedir ve bellek kullanımı \(O(N)\)'dir. Hedef \(N=10^7\) için herhangi bir sayının en fazla sekiz farklı asal çarpanı vardır; dolayısıyla altküme dolaşımı pratikte küçüktür. Paralelleştirme sadece gerçek çalışma süresini azaltır, asimptotik işi değiştirmez.
Para cada entero \(M \ge 2\), se define
$$R(M)=\sum_{\substack{1\le p \lt q \le M\\ p+q \ge M\\ \gcd(p,q)=1}}\frac{1}{pq}, \qquad S(N)=\sum_{M=2}^{N}R(M).$$
Los valores de control son \(S(2)=\tfrac12\), \(S(10)=6.9146825397\ldots\) y \(S(100)=58.2962380621\ldots\). El objetivo es calcular \(S(10^7)\) con cuatro decimales. Evaluar directamente todos los ternas \((M,p,q)\) es inviable, así que la solución reorganiza la suma y sustituye la condición de coprimalidad por inclusión-exclusión de Möbius.
Fijemos un par coprimo \(1 \le p \lt q \le N\). Ese par contribuye a \(R(M)\) exactamente cuando \(q \le M\) y \(M \le p+q\). Además, \(M\) nunca supera a \(N\), así que los valores admisibles son
$$q \le M \le \min(N,p+q).$$
Por tanto, la multiplicidad del par es
$$w_N(p,q)=\min(N,p+q)-q+1.$$
Al intercambiar el orden de sumación obtenemos
$$S(N)=\sum_{\substack{1\le p \lt q \le N\\ \gcd(p,q)=1}}\frac{w_N(p,q)}{pq} =\sum_{q=2}^{N}\frac{1}{q}\sum_{\substack{1\le p \lt q\\ \gcd(p,q)=1}}\frac{w_N(p,q)}{p}.$$
Esto reduce el problema a sumar contribuciones independientes indexadas por el denominador mayor \(q\).
Definimos la suma armónica coprima y la función de conteo
$$A_q(m)=\sum_{\substack{1\le p \le m\\ \gcd(p,q)=1}}\frac{1}{p}, \qquad C_q(m)=\sum_{\substack{1\le p \le m\\ \gcd(p,q)=1}}1.$$
Usando la identidad
$$\mathbf{1}_{\gcd(p,q)=1}=\sum_{d \mid \gcd(p,q)}\mu(d),$$
ambas cantidades pueden escribirse como sumas sobre divisores de \(q\). Si \(H_t=\sum_{j=1}^{t}\frac{1}{j}\) es el número armónico \(t\)-ésimo, entonces
$$A_q(m)=\sum_{d \mid q}\frac{\mu(d)}{d}H_{\lfloor m/d \rfloor}, \qquad C_q(m)=\sum_{d \mid q}\mu(d)\left\lfloor \frac{m}{d}\right\rfloor.$$
Solo importan los divisores libres de cuadrados, porque \(\mu(d)=0\) cuando \(d\) tiene un factor primo repetido. Por eso la implementación factoriza \(q\) en primos distintos y enumera todos los subconjuntos cuadrado-libres con su signo de Möbius.
Si \(q \le \lfloor N/2 \rfloor\), entonces para todo \(p \lt q\) se cumple
$$p+q \le (q-1)+q = 2q-1 \le N-1,$$
así que \(\min(N,p+q)=p+q\). En consecuencia,
$$w_N(p,q)=p+1.$$
La contribución de ese \(q\) es
$$\frac{1}{q}\sum_{\substack{1\le p \lt q\\ \gcd(p,q)=1}}\frac{p+1}{p} =\frac{1}{q}\left(\sum_{\substack{1\le p \lt q\\ \gcd(p,q)=1}}1+\sum_{\substack{1\le p \lt q\\ \gcd(p,q)=1}}\frac{1}{p}\right).$$
La primera suma vale \(\varphi(q)\), porque en \(1,2,\dots,q-1\) hay exactamente \(\varphi(q)\) enteros coprimos con \(q\). Por tanto,
$$T_N(q)=\frac{\varphi(q)+A_q(q-1)}{q}\qquad \text{para } q \le \left\lfloor \frac{N}{2}\right\rfloor.$$
Ahora supongamos \(q > \lfloor N/2 \rfloor\) y definamos \(a=N-q\). Entonces \(0 \le a \lt q\), y el peso cambia al cruzar \(p=a\):
$$w_N(p,q)= \begin{cases} p+1, & 1 \le p \le a,\\ a+1, & a \lt p \lt q. \end{cases}$$
Así, la contribución se separa en dos partes:
$$T_N(q)=\frac{1}{q}\left(\sum_{\substack{1\le p \le a\\ \gcd(p,q)=1}}\left(1+\frac{1}{p}\right) +(a+1)\sum_{\substack{a \lt p \lt q\\ \gcd(p,q)=1}}\frac{1}{p}\right).$$
Si escribimos la segunda suma de recíprocos como \(A_q(q-1)-A_q(a)\), resulta
$$T_N(q)=\frac{C_q(a)+A_q(a)+(a+1)\bigl(A_q(q-1)-A_q(a)\bigr)}{q}.$$
Esta es exactamente la rama que usan las implementaciones cuando \(q\) supera la mitad de \(N\).
Al combinar ambos rangos se obtiene
$$\boxed{ S(N)= \sum_{q=2}^{\lfloor N/2 \rfloor}\frac{\varphi(q)+A_q(q-1)}{q} + \sum_{q=\lfloor N/2 \rfloor+1}^{N} \frac{C_q(N-q)+A_q(N-q)+(N-q+1)\bigl(A_q(q-1)-A_q(N-q)\bigr)}{q} .}$$
Con esto, la suma triple original queda reducida a tres ingredientes precalculables: los valores de la función totiente, los prefijos armónicos y las factorizaciones de \(2,\dots,N\).
Para \(q=5\) seguimos en la primera rama. Como \(\varphi(5)=4\) y
$$A_5(4)=1+\frac12+\frac13+\frac14=\frac{25}{12},$$
tenemos
$$T_{10}(5)=\frac{4+\frac{25}{12}}{5}=\frac{73}{60}=1.216666\ldots$$
Para \(q=8\) entramos en la segunda rama, con \(a=10-8=2\). En este caso
$$C_8(2)=1,\qquad A_8(2)=1,\qquad A_8(7)-A_8(2)=\frac13+\frac15+\frac17=\frac{71}{105},$$
de modo que
$$T_{10}(8)=\frac{1+1+3\cdot \frac{71}{105}}{8}=\frac{141}{280}=0.503571\ldots$$
Al sumar todos los términos desde \(q=2\) hasta \(q=10\) obtenemos
$$S(10)=\frac{3485}{504}=6.914682539682\ldots,$$
exactamente el valor de control esperado.
Las implementaciones en C++, Python y Java siguen el mismo esquema matemático. Primero construyen una tabla del menor factor primo hasta \(N\); con ella se puede factorizar cada \(q\) rápidamente y, al mismo tiempo, calcular \(\varphi(q)\) en tiempo lineal. Después precalculan el arreglo de prefijos armónicos \(H_0,H_1,\dots,H_N\). A continuación, cada \(q\) se procesa de forma independiente: se factoriza \(q\), se enumeran los divisores cuadrado-libres inducidos por sus primos distintos, se evalúan las fórmulas de Möbius para \(A_q\) y \(C_q\), se aplica la rama correcta y se añade la contribución al total.
Las versiones en C++ y Java además usan suma compensada para reducir el error de redondeo, y la versión en C++ puede dividir el bucle exterior en bloques paralelos. Esas decisiones mejoran estabilidad numérica y tiempo de ejecución, pero no alteran la fórmula matemática.
Construir la tabla del menor factor primo, la tabla de totientes y el prefijo armónico cuesta \(O(N)\) tiempo y \(O(N)\) memoria. Para cada \(q\), la factorización se recupera desde la tabla precalculada y la expansión de Möbius recorre exactamente \(2^{\omega(q)}\) divisores cuadrado-libres, donde \(\omega(q)\) es el número de factores primos distintos de \(q\). Por ello, el tiempo total es
$$O\!\left(N+\sum_{q=2}^{N}2^{\omega(q)}\right),$$
con memoria \(O(N)\). Para el objetivo concreto \(N=10^7\), cualquier entero tiene a lo sumo ocho factores primos distintos, así que la enumeración de subconjuntos sigue siendo pequeña en la práctica. El paralelismo reduce solo el tiempo de reloj, no el trabajo asintótico.
对每个整数 \(M \ge 2\)，定义
$$R(M)=\sum_{\substack{1\le p \lt q \le M\\ p+q \ge M\\ \gcd(p,q)=1}}\frac{1}{pq}, \qquad S(N)=\sum_{M=2}^{N}R(M).$$
已知校验值为 \(S(2)=\tfrac12\)、\(S(10)=6.9146825397\ldots\)、\(S(100)=58.2962380621\ldots\)。目标是把 \(S(10^7)\) 计算到小数点后四位。若直接枚举所有 \((M,p,q)\) 三元组，代价会高得不可接受，因此实现采用换序求和，并用 Möbius 反演来处理互素条件。
固定一个互素对 \(1 \le p \lt q \le N\)。它对 \(R(M)\) 有贡献，当且仅当同时满足 \(q \le M\) 与 \(M \le p+q\)。再加上 \(M \le N\) 的限制，可行范围是
$$q \le M \le \min(N,p+q).$$
因此该对 \((p,q)\) 出现的次数为
$$w_N(p,q)=\min(N,p+q)-q+1.$$
把求和顺序改写后，有
$$S(N)=\sum_{\substack{1\le p \lt q \le N\\ \gcd(p,q)=1}}\frac{w_N(p,q)}{pq} =\sum_{q=2}^{N}\frac{1}{q}\sum_{\substack{1\le p \lt q\\ \gcd(p,q)=1}}\frac{w_N(p,q)}{p}.$$
这样问题就被拆成按较大分母 \(q\) 索引的一系列独立贡献。
定义互素倒数和与互素计数函数
$$A_q(m)=\sum_{\substack{1\le p \le m\\ \gcd(p,q)=1}}\frac{1}{p}, \qquad C_q(m)=\sum_{\substack{1\le p \le m\\ \gcd(p,q)=1}}1.$$
利用恒等式
$$\mathbf{1}_{\gcd(p,q)=1}=\sum_{d \mid \gcd(p,q)}\mu(d),$$
这两个量都能改写成关于 \(q\) 的因子和。若 \(H_t=\sum_{j=1}^{t}\frac{1}{j}\) 是第 \(t\) 个调和数，则
$$A_q(m)=\sum_{d \mid q}\frac{\mu(d)}{d}H_{\lfloor m/d \rfloor}, \qquad C_q(m)=\sum_{d \mid q}\mu(d)\left\lfloor \frac{m}{d}\right\rfloor.$$
由于非平方自由因子满足 \(\mu(d)=0\)，真正需要枚举的只有平方自由因子。所以实现先分解出 \(q\) 的不同质因子，再枚举这些质因子的所有子集，并附上相应的 Möbius 符号。
当 \(q \le \lfloor N/2 \rfloor\) 时，对任意 \(p \lt q\) 都有
$$p+q \le (q-1)+q = 2q-1 \le N-1,$$
因此 \(\min(N,p+q)=p+q\)。于是
$$w_N(p,q)=p+1.$$
对应的贡献为
$$\frac{1}{q}\sum_{\substack{1\le p \lt q\\ \gcd(p,q)=1}}\frac{p+1}{p} =\frac{1}{q}\left(\sum_{\substack{1\le p \lt q\\ \gcd(p,q)=1}}1+\sum_{\substack{1\le p \lt q\\ \gcd(p,q)=1}}\frac{1}{p}\right).$$
第一项正是 \(\varphi(q)\)，因为在 \(1,2,\dots,q-1\) 中恰有 \(\varphi(q)\) 个整数与 \(q\) 互素。因此
$$T_N(q)=\frac{\varphi(q)+A_q(q-1)}{q}\qquad \text{当 } q \le \left\lfloor \frac{N}{2}\right\rfloor.$$
现在设 \(q > \lfloor N/2 \rfloor\)，并记 \(a=N-q\)。此时 \(0 \le a \lt q\)，权重在 \(p=a\) 处分段：
$$w_N(p,q)= \begin{cases} p+1, & 1 \le p \le a,\\ a+1, & a \lt p \lt q. \end{cases}$$
所以贡献拆成两部分：
$$T_N(q)=\frac{1}{q}\left(\sum_{\substack{1\le p \le a\\ \gcd(p,q)=1}}\left(1+\frac{1}{p}\right) +(a+1)\sum_{\substack{a \lt p \lt q\\ \gcd(p,q)=1}}\frac{1}{p}\right).$$
把第二段写成 \(A_q(q-1)-A_q(a)\)，得到
$$T_N(q)=\frac{C_q(a)+A_q(a)+(a+1)\bigl(A_q(q-1)-A_q(a)\bigr)}{q}.$$
这就是实现中跨过中点后所使用的精确公式。
合并两个区间，可得
$$\boxed{ S(N)= \sum_{q=2}^{\lfloor N/2 \rfloor}\frac{\varphi(q)+A_q(q-1)}{q} + \sum_{q=\lfloor N/2 \rfloor+1}^{N} \frac{C_q(N-q)+A_q(N-q)+(N-q+1)\bigl(A_q(q-1)-A_q(N-q)\bigr)}{q} .}$$
原来的三层求和因此被压缩成三个可预处理的对象：欧拉函数值、调和前缀和、以及 \(2,\dots,N\) 的质因数分解信息。
当 \(q=5\) 时，还处于第一种情况。因为 \(\varphi(5)=4\)，并且
$$A_5(4)=1+\frac12+\frac13+\frac14=\frac{25}{12},$$
所以
$$T_{10}(5)=\frac{4+\frac{25}{12}}{5}=\frac{73}{60}=1.216666\ldots$$
当 \(q=8\) 时，进入第二种情况，此时 \(a=10-8=2\)。有
$$C_8(2)=1,\qquad A_8(2)=1,\qquad A_8(7)-A_8(2)=\frac13+\frac15+\frac17=\frac{71}{105},$$
因此
$$T_{10}(8)=\frac{1+1+3\cdot \frac{71}{105}}{8}=\frac{141}{280}=0.503571\ldots$$
把 \(q=2\) 到 \(10\) 的所有贡献相加，得到
$$S(10)=\frac{3485}{504}=6.914682539682\ldots,$$
与校验值完全一致。
C++、Python 和 Java 实现遵循同一条数学流程。首先预计算到 \(N\) 为止的最小质因子表，这样每个 \(q\) 都能被快速分解，同时也能在线性时间内得到 \(\varphi(q)\)。然后构造调和前缀数组 \(H_0,H_1,\dots,H_N\)。之后对每个 \(q\) 独立处理：分解 \(q\)，枚举由不同质因子生成的所有平方自由因子及其 Möbius 符号，计算 \(A_q\) 与 \(C_q\)，套用对应分支，把结果加入总和。
C++ 与 Java 版本还使用补偿求和来减小浮点累积误差，C++ 版本还可以把外层循环切成并行区块。这些工程细节不会改变数学公式，只是为了在目标规模上获得更稳定、更快的运行效果。
建立最小质因子表、欧拉函数表和调和前缀数组需要 \(O(N)\) 时间与 \(O(N)\) 空间。对每个 \(q\)，分解可以直接从预处理表恢复，而 Möbius 展开恰好枚举 \(2^{\omega(q)}\) 个平方自由因子，其中 \(\omega(q)\) 表示 \(q\) 的不同质因子个数。因此总时间为
$$O\!\left(N+\sum_{q=2}^{N}2^{\omega(q)}\right),$$
空间复杂度为 \(O(N)\)。对于目标 \(N=10^7\)，任意整数最多只有八个不同质因子，所以子集枚举在实际中非常小。并行化只会降低墙钟时间，不会改变渐近工作量。
Для каждого целого \(M \ge 2\) задаются
$$R(M)=\sum_{\substack{1\le p \lt q \le M\\ p+q \ge M\\ \gcd(p,q)=1}}\frac{1}{pq}, \qquad S(N)=\sum_{M=2}^{N}R(M).$$
Контрольные значения равны \(S(2)=\tfrac12\), \(S(10)=6.9146825397\ldots\) и \(S(100)=58.2962380621\ldots\). Требуется найти \(S(10^7)\) с округлением до четырёх знаков после запятой. Прямой перебор всех троек \((M,p,q)\) слишком дорог, поэтому решение меняет порядок суммирования и заменяет условие взаимной простоты формулами Мёбиуса.
Зафиксируем взаимно простую пару \(1 \le p \lt q \le N\). Она входит в \(R(M)\) тогда и только тогда, когда выполнены неравенства \(q \le M\) и \(M \le p+q\). С учётом ограничения \(M \le N\) допустимые значения имеют вид
$$q \le M \le \min(N,p+q).$$
Значит, кратность этой пары равна
$$w_N(p,q)=\min(N,p+q)-q+1.$$
После перестановки порядка суммирования получаем
$$S(N)=\sum_{\substack{1\le p \lt q \le N\\ \gcd(p,q)=1}}\frac{w_N(p,q)}{pq} =\sum_{q=2}^{N}\frac{1}{q}\sum_{\substack{1\le p \lt q\\ \gcd(p,q)=1}}\frac{w_N(p,q)}{p}.$$
Тем самым задача распадается на независимые вклады, зависящие только от большего знаменателя \(q\).
Введём сумму обратных по взаимно простым числам и функцию подсчёта
$$A_q(m)=\sum_{\substack{1\le p \le m\\ \gcd(p,q)=1}}\frac{1}{p}, \qquad C_q(m)=\sum_{\substack{1\le p \le m\\ \gcd(p,q)=1}}1.$$
Используя тождество
$$\mathbf{1}_{\gcd(p,q)=1}=\sum_{d \mid \gcd(p,q)}\mu(d),$$
обе величины можно записать через делители числа \(q\). Если \(H_t=\sum_{j=1}^{t}\frac{1}{j}\) обозначает \(t\)-е гармоническое число, то
$$A_q(m)=\sum_{d \mid q}\frac{\mu(d)}{d}H_{\lfloor m/d \rfloor}, \qquad C_q(m)=\sum_{d \mid q}\mu(d)\left\lfloor \frac{m}{d}\right\rfloor.$$
Учитываются только квадратсвободные делители, потому что для делителей с повторяющимся простым множителем \(\mu(d)=0\). Поэтому реализация разлагает \(q\) на различные простые множители и перебирает все квадратсвободные подмножества с нужными знаками Мёбиуса.
Если \(q \le \lfloor N/2 \rfloor\), то для любого \(p \lt q\)
$$p+q \le (q-1)+q = 2q-1 \le N-1,$$
поэтому \(\min(N,p+q)=p+q\). Следовательно,
$$w_N(p,q)=p+1.$$
Вклад такого \(q\) равен
$$\frac{1}{q}\sum_{\substack{1\le p \lt q\\ \gcd(p,q)=1}}\frac{p+1}{p} =\frac{1}{q}\left(\sum_{\substack{1\le p \lt q\\ \gcd(p,q)=1}}1+\sum_{\substack{1\le p \lt q\\ \gcd(p,q)=1}}\frac{1}{p}\right).$$
Первая сумма есть \(\varphi(q)\), потому что среди чисел \(1,2,\dots,q-1\) ровно \(\varphi(q)\) взаимно просты с \(q\). Значит,
$$T_N(q)=\frac{\varphi(q)+A_q(q-1)}{q}\qquad \text{при } q \le \left\lfloor \frac{N}{2}\right\rfloor.$$
Теперь предположим, что \(q > \lfloor N/2 \rfloor\), и положим \(a=N-q\). Тогда \(0 \le a \lt q\), а вес меняется в точке \(p=a\):
$$w_N(p,q)= \begin{cases} p+1, & 1 \le p \le a,\\ a+1, & a \lt p \lt q. \end{cases}$$
Поэтому вклад раскладывается так:
$$T_N(q)=\frac{1}{q}\left(\sum_{\substack{1\le p \le a\\ \gcd(p,q)=1}}\left(1+\frac{1}{p}\right) +(a+1)\sum_{\substack{a \lt p \lt q\\ \gcd(p,q)=1}}\frac{1}{p}\right).$$
Если вторую сумму записать как \(A_q(q-1)-A_q(a)\), получаем
$$T_N(q)=\frac{C_q(a)+A_q(a)+(a+1)\bigl(A_q(q-1)-A_q(a)\bigr)}{q}.$$
Именно эту формулу используют реализации после середины диапазона.
Объединяя оба случая, получаем
$$\boxed{ S(N)= \sum_{q=2}^{\lfloor N/2 \rfloor}\frac{\varphi(q)+A_q(q-1)}{q} + \sum_{q=\lfloor N/2 \rfloor+1}^{N} \frac{C_q(N-q)+A_q(N-q)+(N-q+1)\bigl(A_q(q-1)-A_q(N-q)\bigr)}{q} .}$$
Исходная тройная сумма тем самым сведена к трём заранее вычислимым объектам: значениям функции Эйлера, гармоническим префиксам и факторизациям чисел \(2,\dots,N\).
Для \(q=5\) мы ещё находимся в первом случае. Поскольку \(\varphi(5)=4\) и
$$A_5(4)=1+\frac12+\frac13+\frac14=\frac{25}{12},$$
то
$$T_{10}(5)=\frac{4+\frac{25}{12}}{5}=\frac{73}{60}=1.216666\ldots$$
Для \(q=8\) действует второй случай, где \(a=10-8=2\). Здесь
$$C_8(2)=1,\qquad A_8(2)=1,\qquad A_8(7)-A_8(2)=\frac13+\frac15+\frac17=\frac{71}{105},$$
поэтому
$$T_{10}(8)=\frac{1+1+3\cdot \frac{71}{105}}{8}=\frac{141}{280}=0.503571\ldots$$
Суммирование вкладов при \(q=2,\dots,10\) даёт
$$S(10)=\frac{3485}{504}=6.914682539682\ldots,$$
что точно совпадает с контрольным значением.
Реализации на C++, Python и Java следуют одной и той же математической схеме. Сначала строится таблица минимального простого делителя до \(N\); благодаря ей каждое \(q\) быстро факторизуется, и одновременно за линейное время вычисляется \(\varphi(q)\). Затем заранее заполняется массив гармонических префиксов \(H_0,H_1,\dots,H_N\). После этого каждое значение \(q\) обрабатывается независимо: разложить \(q\), перечислить все квадратсвободные делители, полученные из различных простых факторов, вычислить по формулам Мёбиуса величины \(A_q\) и \(C_q\), применить нужную ветвь и добавить вклад к общей сумме.
Версии на C++ и Java дополнительно используют компенсированное суммирование, чтобы уменьшить накопление ошибок округления, а версия на C++ может делить внешний цикл на параллельные блоки. Эти инженерные детали не меняют математику, а лишь улучшают устойчивость и время работы на целевом размере.
Построение таблицы минимального простого делителя, таблицы значений \(\varphi\) и массива гармонических префиксов требует \(O(N)\) времени и \(O(N)\) памяти. Для каждого \(q\) разложение берётся из предвычисленной таблицы, а разложение по Мёбиусу перечисляет ровно \(2^{\omega(q)}\) квадратсвободных делителей, где \(\omega(q)\) означает число различных простых делителей \(q\). Следовательно, общее время равно
$$O\!\left(N+\sum_{q=2}^{N}2^{\omega(q)}\right),$$
а память остаётся \(O(N)\). Для целевого \(N=10^7\) любое число имеет не более восьми различных простых делителей, поэтому перебор подмножеств остаётся небольшим на практике. Параллелизм уменьшает только реальное время исполнения, но не асимптотическую трудоёмкость.
لكل عدد صحيح \(M \ge 2\) نعرّف
$$R(M)=\sum_{\substack{1\le p \lt q \le M\\ p+q \ge M\\ \gcd(p,q)=1}}\frac{1}{pq}, \qquad S(N)=\sum_{M=2}^{N}R(M).$$
قيم التحقق هي \(S(2)=\tfrac12\)، و\(S(10)=6.9146825397\ldots\)، و\(S(100)=58.2962380621\ldots\). المطلوب هو حساب \(S(10^7)\) مع التقريب إلى أربع منازل عشرية. الفحص المباشر لكل الثلاثيات \((M,p,q)\) مكلف جداً، لذلك يعيد الحل ترتيب الجمع ويستبدل شرط التباين الأولي بصياغة تعتمد على دالة Möbius.
لنثبت زوجاً أولياً فيما بينه \(1 \le p \lt q \le N\). هذا الزوج يساهم في \(R(M)\) بالضبط عندما يتحقق الشرطان \(q \le M\) و\(M \le p+q\). ومع القيد الإضافي \(M \le N\) تصبح القيم الممكنة
$$q \le M \le \min(N,p+q).$$
إذن عدد مرات ظهور الزوج هو
$$w_N(p,q)=\min(N,p+q)-q+1.$$
وبعد تبديل ترتيب الجمع نحصل على
$$S(N)=\sum_{\substack{1\le p \lt q \le N\\ \gcd(p,q)=1}}\frac{w_N(p,q)}{pq} =\sum_{q=2}^{N}\frac{1}{q}\sum_{\substack{1\le p \lt q\\ \gcd(p,q)=1}}\frac{w_N(p,q)}{p}.$$
وهكذا تتحول المسألة إلى مجموع مساهمات مستقلة مفهرسة بالقيمة الأكبر \(q\).
نعرّف مجموع المقلوبات على الأعداد الأولية نسبياً مع \(q\)، وكذلك دالة العد:
$$A_q(m)=\sum_{\substack{1\le p \le m\\ \gcd(p,q)=1}}\frac{1}{p}, \qquad C_q(m)=\sum_{\substack{1\le p \le m\\ \gcd(p,q)=1}}1.$$
وباستخدام الهوية
$$\mathbf{1}_{\gcd(p,q)=1}=\sum_{d \mid \gcd(p,q)}\mu(d),$$
يمكن إعادة كتابة الكميتين على صورة مجموعات على قواسم \(q\). إذا كتبنا \(H_t=\sum_{j=1}^{t}\frac{1}{j}\) للعدد التوافقي رقم \(t\)، فإن
$$A_q(m)=\sum_{d \mid q}\frac{\mu(d)}{d}H_{\lfloor m/d \rfloor}, \qquad C_q(m)=\sum_{d \mid q}\mu(d)\left\lfloor \frac{m}{d}\right\rfloor.$$
لا تسهم إلا القواسم الخالية من المربعات، لأن \(\mu(d)=0\) إذا احتوى \(d\) على مربع أولي. لهذا السبب تقوم الخوارزمية بتحليل \(q\) إلى عوامله الأولية المميزة ثم تعدّد جميع القواسم الخالية من المربعات الناتجة من هذه العوامل مع إشارات Möbius الموافقة.
إذا كان \(q \le \lfloor N/2 \rfloor\)، فعند كل \(p \lt q\) يكون
$$p+q \le (q-1)+q = 2q-1 \le N-1,$$
ومن ثم \(\min(N,p+q)=p+q\). وبالتالي
$$w_N(p,q)=p+1.$$
وعندئذ تصبح مساهمة هذا \(q\)
$$\frac{1}{q}\sum_{\substack{1\le p \lt q\\ \gcd(p,q)=1}}\frac{p+1}{p} =\frac{1}{q}\left(\sum_{\substack{1\le p \lt q\\ \gcd(p,q)=1}}1+\sum_{\substack{1\le p \lt q\\ \gcd(p,q)=1}}\frac{1}{p}\right).$$
المجموع الأول يساوي \(\varphi(q)\)، لأن بين الأعداد \(1,2,\dots,q-1\) يوجد بالضبط \(\varphi(q)\) عدداً أولياً نسبياً مع \(q\). إذن
$$T_N(q)=\frac{\varphi(q)+A_q(q-1)}{q}$$
وذلك عندما \(q \le \left\lfloor \frac{N}{2}\right\rfloor\).
الآن افترض أن \(q > \lfloor N/2 \rfloor\)، ولنضع \(a=N-q\). عندئذ \(0 \le a \lt q\)، ويتغير الوزن عند النقطة \(p=a\):
$$w_N(p,q)= \begin{cases} p+1, & 1 \le p \le a,\\ a+1, & a \lt p \lt q. \end{cases}$$
ومن ثم تنقسم المساهمة إلى جزأين:
$$T_N(q)=\frac{1}{q}\left(\sum_{\substack{1\le p \le a\\ \gcd(p,q)=1}}\left(1+\frac{1}{p}\right) +(a+1)\sum_{\substack{a \lt p \lt q\\ \gcd(p,q)=1}}\frac{1}{p}\right).$$
وبكتابة المجموع الثاني على الصورة \(A_q(q-1)-A_q(a)\) نحصل على
$$T_N(q)=\frac{C_q(a)+A_q(a)+(a+1)\bigl(A_q(q-1)-A_q(a)\bigr)}{q}.$$
وهذا هو الفرع نفسه الذي تستخدمه الخوارزمية عندما يتجاوز \(q\) منتصف المجال.
بجمع الحالتين نحصل على
$$\boxed{ S(N)= \sum_{q=2}^{\lfloor N/2 \rfloor}\frac{\varphi(q)+A_q(q-1)}{q} + \sum_{q=\lfloor N/2 \rfloor+1}^{N} \frac{C_q(N-q)+A_q(N-q)+(N-q+1)\bigl(A_q(q-1)-A_q(N-q)\bigr)}{q} .}$$
وبذلك ينخفض التعقيد من جمع ثلاثي مباشر إلى ثلاثة مكونات قابلة للتهيئة المسبقة: قيم دالة أويلر، والبادئات التوافقية، وتحليل الأعداد \(2,\dots,N\) إلى عواملها الأولية.
عند \(q=5\) ما زلنا في الفرع الأول. بما أن \(\varphi(5)=4\) و
$$A_5(4)=1+\frac12+\frac13+\frac14=\frac{25}{12},$$
فإن
$$T_{10}(5)=\frac{4+\frac{25}{12}}{5}=\frac{73}{60}=1.216666\ldots$$
أما عند \(q=8\) فننتقل إلى الفرع الثاني مع \(a=10-8=2\). وهنا
$$C_8(2)=1,\qquad A_8(2)=1,\qquad A_8(7)-A_8(2)=\frac13+\frac15+\frac17=\frac{71}{105},$$
ولذلك
$$T_{10}(8)=\frac{1+1+3\cdot \frac{71}{105}}{8}=\frac{141}{280}=0.503571\ldots$$
ومجموع جميع المساهمات من \(q=2\) إلى \(q=10\) يساوي
$$S(10)=\frac{3485}{504}=6.914682539682\ldots,$$
وهو مطابق تماماً لقيمة التحقق.
تتبع تطبيقات C++ وPython وJava المسار الرياضي نفسه. أولاً تُبنى حتى \(N\) جدولٌ لأصغر عامل أولي، مما يسمح بتحليل كل \(q\) بسرعة، كما يتيح حساب \(\varphi(q)\) في زمن خطي. بعد ذلك تُحسب البادئة التوافقية \(H_0,H_1,\dots,H_N\). ثم يُعالج كل \(q\) بصورة مستقلة: تحليل \(q\)، وتوليد جميع القواسم الخالية من المربعات الناتجة من عوامله الأولية المميزة، وتقييم صيغ Möbius الخاصة بـ \(A_q\) و\(C_q\)، واختيار الفرع المناسب، ثم إضافة الناتج إلى المجموع الكلي.
تستخدم نسختا C++ وJava أيضاً جمعاً تعويضياً لتقليل تراكم الخطأ العددي، ويمكن لنسخة C++ تقسيم الحلقة الخارجية إلى كتل متوازية. هذه الخيارات الهندسية لا تغير الرياضيات، لكنها تحسن الاستقرار العددي وزمن التنفيذ عند الحجم المطلوب.
بناء جدول أصغر عامل أولي، وجدول \(\varphi\)، ومصفوفة البادئات التوافقية يتطلب \(O(N)\) زمناً و\(O(N)\) ذاكرة. ولكل \(q\)، يُستخرج التحليل من الجدول المحسوب مسبقاً، بينما يمر توسيع Möbius على عدد يساوي تماماً \(2^{\omega(q)}\) من القواسم الخالية من المربعات، حيث تمثل \(\omega(q)\) عدد العوامل الأولية المميزة لـ \(q\). لذلك يكون الزمن الكلي
$$O\!\left(N+\sum_{q=2}^{N}2^{\omega(q)}\right),$$
مع ذاكرة \(O(N)\). وبالنسبة للهدف \(N=10^7\)، فإن أي عدد لا يملك أكثر من ثمانية عوامل أولية مميزة، لذا يبقى تعداد المجموعات الفرعية صغيراً عملياً. التنفيذ المتوازي يخفض زمن الساعة فقط ولا يغير العمل الأسيمبتوطي.