For coprime positive integers \(a\le b\), the quantity studied here is built from Bézout relations between \(a\) and \(b\). In the Project Euler instance, the inputs are Mersenne numbers \(M_u=2^u-1\) and \(M_v=2^v-1\), where \(u=r^5\) and \(v=s^5\) for distinct primes \(r\lt s\lt 1000\). The goal is to evaluate the contribution for every such prime pair and sum all contributions modulo \(10^9+7\).
Directly constructing \(2^u-1\) and \(2^v-1\) is hopeless because the exponents are enormous. The key observation is that the Euclidean algorithm on Mersenne numbers can be simulated on the exponents alone, while every coefficient that matters is tracked only modulo the final modulus.
Write \(M_n=2^n-1\). For two coprime positive integers \(a\le b\), choose the least positive inverse \(x\) of \(a\) modulo \(b\):
$$a x \equiv 1 \pmod{b},\qquad 1\le x \lt b.$$
Then
$$y=\frac{a x - 1}{b}$$
is a positive integer and satisfies \(a x - b y = 1\). The complementary positive solution to the opposite sign equation is \((b-x,\ a-y)\), so the implementations evaluate
$$P(a,b)=2\min\left(x+y-1,\ (b-x)+(a-y)-1\right).$$
The problem is therefore to compute \(P(M_u,M_v)\) quickly for very large exponents \(u\) and \(v\).
Let \(a=M_u\) and \(b=M_v\) with \(u\lt v\). If \(u\) and \(v\) are coprime, then the classical identity
$$\gcd(M_u,M_v)=M_{\gcd(u,v)}$$
gives
$$\gcd(M_u,M_v)=M_1=1.$$
This is exactly the situation here, because distinct prime fifth powers remain coprime:
$$\gcd(r^5,s^5)=1\qquad (r\ne s).$$
So every target pair is valid for the Bézout setup, even though the numbers themselves are astronomically large.
If \(A=qB+r\) with \(0\le r \lt B\), then
$$2^A-1=2^r\left(2^{qB}-1\right)+\left(2^r-1\right).$$
After factoring \(2^{qB}-1=(2^B-1)\sum_{j=0}^{q-1}2^{jB}\), this becomes
$$M_A = Q\,M_B + M_r,$$
with quotient
$$Q = 2^r\sum_{j=0}^{q-1}2^{jB}.$$
Therefore the remainder when dividing \(M_A\) by \(M_B\) is exactly \(M_r\). The full remainder sequence for \((M_A,M_B)\) is obtained by applying ordinary Euclid to the exponent pair \((A,B)\). This is the compression that makes the problem tractable.
Each Euclidean step needs the quotient \(Q\) only modulo \(10^9+7\). Using the formula above, we only need
$$Q \equiv 2^r\sum_{j=0}^{q-1}\left(2^B\right)^j \pmod{10^9+7}.$$
The implementation evaluates this geometric sum by divide and conquer, simultaneously tracking a power and the corresponding partial sum. As a result, every quotient update is handled with modular arithmetic and logarithmic recursion depth, not with gigantic integers.
Feeding these modular quotients into the extended Euclidean recurrence yields the inverse of a Mersenne remainder modulo another Mersenne number, but represented only modulo \(10^9+7\).
Let
$$w = v \bmod u,$$
and let \(t\) be the inverse of \(w\) modulo \(u\):
$$w t \equiv 1 \pmod{u},\qquad 1\le t \lt u.$$
Now consider
$$R = 1 + 2^w + 2^{2w} + \cdots + 2^{(t-1)w}.$$
Then
$$\begin{aligned} (2^w-1)R &= 2^{tw}-1 \\ &\equiv 2^1-1 \equiv 1 \pmod{2^u-1}, \end{aligned}$$
so \(R\) is an inverse of \(M_w\) modulo \(M_u\). Because \(M_v\equiv M_w \pmod{M_u}\), the same value also acts as an inverse of \(M_v\) modulo \(M_u\).
The opposite Bézout sign branch is obtained by replacing \(R\) with \(M_u-R\). The implementation chooses between these two branches with the criterion
$$2t\le u,$$
which is equivalent to choosing between the \(t\)-term geometric sum and its \((u-t)\)-term complement inside \(M_u=1+2+\cdots+2^{u-1}\).
Let \(d\) be the branch selected in the previous step. Then one of the two relations
$$M_v d = M_u x \pm 1.$$
holds, depending on the chosen sign. Once \(x\) is recovered, the contribution is
$$P(M_u,M_v)=2(x+d-1).$$
This is the same minimum as in the generic formula for \(P(a,b)\); the only difference is that the implementation works with the coefficient attached to \(M_v\), because that is the quantity naturally produced by the Mersenne inverse computation.
Take \(u=3\) and \(v=5\). Then
$$M_3=7,\qquad M_5=31,\qquad w=5\bmod 3=2.$$
The inverse of \(w=2\) modulo \(u=3\) is \(t=2\), since
$$2\cdot 2 \equiv 1 \pmod{3}.$$
Hence
$$R = 1 + 2^2 = 5,$$
and indeed
$$M_2 R = 3\cdot 5 = 15 \equiv 1 \pmod{7}.$$
Because \(2t=4>3\), the smaller branch uses the complement
$$d = M_3 - R = 7-5=2.$$
Now
$$x=\frac{M_5 d + 1}{M_3}=\frac{31\cdot 2 + 1}{7}=9,$$
so
$$P(M_3,M_5)=2(x+d-1)=2(9+2-1)=20.$$
This matches the small verification case built into the implementation.
If \(p_1,p_2,\dots,p_k\) are the primes below \(1000\), the final value is
$$\sum_{1\le i \lt j\le k} P\left(M_{p_i^5},M_{p_j^5}\right)\pmod{10^9+7}.$$
The entire mathematical task is therefore to evaluate each pair contribution via exponent-level Euclid and accumulate the results modulo the final prime modulus.
The C++, Python, and Java implementations all expose the same computation, but the actual number-theoretic work is performed by the C++ implementation. It begins with several small checkpoints: first on ordinary coprime integer pairs, and then on small Mersenne exponents where a direct comparison against the generic Bézout formula is still possible.
After the checkpoints, the implementation generates all primes below \(1000\), raises each of them to the fifth power, and iterates over every unordered pair. For each pair \((u,v)\), it computes \(M_u \bmod 10^9+7\) and \(M_v \bmod 10^9+7\), runs the Euclidean algorithm on the exponent pair, and reconstructs the required quotient information only modulo \(10^9+7\).
The key accelerator is the divide-and-conquer evaluation of geometric sums. Instead of forming the true quotient between two huge Mersenne numbers, the implementation directly produces that quotient modulo \(10^9+7\), which is exactly what the extended Euclidean recurrence needs. The branch decision based on the inverse of \(w\) modulo \(u\) then determines which Bézout sign gives the smaller contribution.
Finally, the contribution of every prime pair is added into a running sum modulo \(10^9+7\). The Python and Java implementations are thin launchers around that same compiled computation, so all three language entries stay synchronized on both the mathematics and the final output.
There are \(168\) primes below \(1000\), so the outer loop evaluates \(\binom{168}{2}=14028\) unordered pairs. For one pair, the exponent-level Euclidean algorithm takes \(O(\log \max(u,v))\) divisions, just as ordinary Euclid does on integers. Each division performs a constant number of modular exponentiations and one divide-and-conquer geometric-sum evaluation, both logarithmic in the relevant exponent sizes. Therefore the total running time is \(O(P^2\operatorname{polylog} U)\) with \(P=168\) and \(U=\max p^5\), while memory usage is \(O(P)\) for the prime list plus \(O(\log U)\) recursion depth.
Für teilerfremde positive ganze Zahlen \(a\le b\) wird hier eine Größe betrachtet, die aus Bézout-Beziehungen zwischen \(a\) und \(b\) entsteht. Im Project-Euler-Fall sind die Eingaben Mersenne-Zahlen \(M_u=2^u-1\) und \(M_v=2^v-1\), wobei \(u=r^5\) und \(v=s^5\) für verschiedene Primzahlen \(r\lt s\lt 1000\) gelten. Gesucht ist die Summe aller Beiträge über diese Primzahlpaare modulo \(10^9+7\).
Ein direktes Konstruieren von \(2^u-1\) und \(2^v-1\) ist ausgeschlossen, weil die Exponenten riesig sind. Die entscheidende Beobachtung ist, dass man den euklidischen Algorithmus für Mersenne-Zahlen vollständig auf Exponentenebene nachbilden kann, während alle relevanten Koeffizienten nur modulo des Endmoduls verfolgt werden.
Schreibe \(M_n=2^n-1\). Für zwei teilerfremde positive ganze Zahlen \(a\le b\) sei \(x\) das kleinste positive Inverse von \(a\) modulo \(b\):
$$a x \equiv 1 \pmod{b},\qquad 1\le x \lt b.$$
Dann ist
$$y=\frac{a x - 1}{b}$$
eine positive ganze Zahl mit \(a x - b y = 1\). Die komplementäre positive Lösung zur Gleichung mit umgekehrtem Vorzeichen ist \((b-x,\ a-y)\). Daher berechnen die Implementierungen
$$P(a,b)=2\min\left(x+y-1,\ (b-x)+(a-y)-1\right).$$
Die eigentliche Aufgabe besteht also darin, \(P(M_u,M_v)\) für sehr große Exponenten \(u\) und \(v\) effizient zu bestimmen.
Setze \(a=M_u\) und \(b=M_v\) mit \(u\lt v\). Sind \(u\) und \(v\) teilerfremd, dann liefert die klassische Identität
$$\gcd(M_u,M_v)=M_{\gcd(u,v)}$$
sofort
$$\gcd(M_u,M_v)=M_1=1.$$
Genau das passiert hier, denn verschiedene fünfte Primzahlpotenzen bleiben teilerfremd:
$$\gcd(r^5,s^5)=1\qquad (r\ne s).$$
Jedes Zielpaar passt also zur Bézout-Struktur, obwohl die Zahlen selbst astronomisch groß sind.
Gilt \(A=qB+r\) mit \(0\le r \lt B\), dann
$$2^A-1=2^r\left(2^{qB}-1\right)+\left(2^r-1\right).$$
Nach dem Ausklammern von \(2^{qB}-1=(2^B-1)\sum_{j=0}^{q-1}2^{jB}\) folgt
$$M_A = Q\,M_B + M_r,$$
wobei der Quotient
$$Q = 2^r\sum_{j=0}^{q-1}2^{jB}$$
ist. Der Rest bei der Division von \(M_A\) durch \(M_B\) ist also genau \(M_r\). Damit erhält man die gesamte Restfolge von \((M_A,M_B)\), indem man den gewöhnlichen euklidischen Algorithmus nur auf das Exponentenpaar \((A,B)\) anwendet.
In jedem Schritt wird der Quotient \(Q\) nur modulo \(10^9+7\) benötigt. Aus der Formel oben genügt also
$$Q \equiv 2^r\sum_{j=0}^{q-1}\left(2^B\right)^j \pmod{10^9+7}.$$
Die Implementierung berechnet diese geometrische Summe per Divide-and-Conquer und verfolgt dabei gleichzeitig eine Potenz sowie die zugehörige Teilsumme. Dadurch bleibt jeder Quotienten-Schritt logarithmisch und arbeitet nur mit modularer Arithmetik statt mit gewaltigen Ganzzahlen.
Setzt man diese modularen Quotienten in die Rekursion des erweiterten Euklid ein, erhält man das Inverse eines Mersenne-Rests modulo einer anderen Mersenne-Zahl, jedoch nur noch in der Endmodulo-Darstellung.
Sei
$$w = v \bmod u,$$
und sei \(t\) das Inverse von \(w\) modulo \(u\):
$$w t \equiv 1 \pmod{u},\qquad 1\le t \lt u.$$
Betrachte nun
$$R = 1 + 2^w + 2^{2w} + \cdots + 2^{(t-1)w}.$$
Dann gilt
$$\begin{aligned} (2^w-1)R &= 2^{tw}-1 \\ &\equiv 2^1-1 \equiv 1 \pmod{2^u-1}, \end{aligned}$$
also ist \(R\) ein Inverses von \(M_w\) modulo \(M_u\). Da außerdem \(M_v\equiv M_w \pmod{M_u}\) gilt, fungiert derselbe Wert auch als Inverses von \(M_v\) modulo \(M_u\).
Der Bézout-Zweig mit dem entgegengesetzten Vorzeichen entsteht durch Ersetzen von \(R\) durch \(M_u-R\). Zwischen diesen beiden Zweigen wählt die Implementierung mit dem Kriterium
$$2t\le u,$$
also zwischen einer geometrischen Summe mit \(t\) Termen und ihrem Komplement mit \(u-t\) Termen in \(M_u=1+2+\cdots+2^{u-1}\).
Sei \(d\) der im vorigen Schritt gewählte Koeffizient. Dann gilt je nach Vorzeichen genau eine der beiden Relationen
$$M_v d = M_u x \pm 1.$$
Sobald \(x\) bestimmt ist, lautet der Beitrag
$$P(M_u,M_v)=2(x+d-1).$$
Das ist dasselbe Minimum wie in der allgemeinen Formel für \(P(a,b)\); lediglich der natürliche Koeffizient aus der Mersenne-Inversenberechnung sitzt hier an \(M_v\) statt an \(M_u\).
Nimm \(u=3\) und \(v=5\). Dann
$$M_3=7,\qquad M_5=31,\qquad w=5\bmod 3=2.$$
Das Inverse von \(w=2\) modulo \(u=3\) ist \(t=2\), denn
$$2\cdot 2 \equiv 1 \pmod{3}.$$
Also
$$R = 1 + 2^2 = 5,$$
und tatsächlich
$$M_2 R = 3\cdot 5 = 15 \equiv 1 \pmod{7}.$$
Weil \(2t=4>3\), nimmt der kleinere Zweig das Komplement
$$d = M_3 - R = 7-5=2.$$
Nun ist
$$x=\frac{M_5 d + 1}{M_3}=\frac{31\cdot 2 + 1}{7}=9,$$
also
$$P(M_3,M_5)=2(x+d-1)=2(9+2-1)=20.$$
Genau dieser Wert erscheint als kleiner Kontrollfall in der Implementierung.
Seien \(p_1,p_2,\dots,p_k\) die Primzahlen unter \(1000\). Dann ist der gesuchte Endwert
$$\sum_{1\le i \lt j\le k} P\left(M_{p_i^5},M_{p_j^5}\right)\pmod{10^9+7}.$$
Die mathematische Hauptaufgabe des Programms besteht also darin, jeden Paarbeitrag per Exponenten-Euklid auszurechnen und alles modulo des Primmoduls zu akkumulieren.
Die C++-, Python- und Java-Implementierungen liefern dasselbe Ergebnis, aber die eigentliche Zahlentheorie steckt in der C++-Implementierung. Zu Beginn werden mehrere kleine Kontrollfälle geprüft: zuerst für gewöhnliche teilerfremde Ganzzahlen und danach für kleine Mersenne-Exponenten, bei denen man das Ergebnis noch direkt mit der allgemeinen Bézout-Formel vergleichen kann.
Danach erzeugt die Implementierung alle Primzahlen unter \(1000\), hebt jede davon in die fünfte Potenz und durchläuft jedes ungeordnete Paar. Für jedes Paar \((u,v)\) werden \(M_u \bmod 10^9+7\) und \(M_v \bmod 10^9+7\) berechnet, der euklidische Algorithmus auf dem Exponentenpaar ausgeführt und die nötige Quotienteninformation nur modulo \(10^9+7\) rekonstruiert.
Der entscheidende Beschleuniger ist die Divide-and-Conquer-Auswertung der geometrischen Summen. Statt den echten Quotienten zweier riesiger Mersenne-Zahlen auszurechnen, produziert die Implementierung diesen Quotienten direkt modulo \(10^9+7\), und genau das benötigt die Rekursion des erweiterten Euklid. Anschließend bestimmt der auf \(w^{-1}\bmod u\) beruhende Zweigtest, welches Bézout-Vorzeichen den kleineren Beitrag liefert.
Am Ende wird jeder Paarbeitrag zu einer laufenden Summe modulo \(10^9+7\) addiert. Die Python- und Java-Implementierungen sind lediglich schlanke Starter um dieselbe kompilierte Berechnung herum, sodass Mathematik und Ausgabe in allen drei Sprachen deckungsgleich bleiben.
Es gibt \(168\) Primzahlen unter \(1000\), also \(\binom{168}{2}=14028\) ungeordnete Paare. Für ein Paar benötigt der Exponenten-Euklid \(O(\log \max(u,v))\) Divisionen, genau wie der gewöhnliche euklidische Algorithmus auf ganzen Zahlen. Jede Division führt eine konstante Anzahl modularer Potenzierungen und eine Divide-and-Conquer-Auswertung einer geometrischen Summe aus, beide logarithmisch in den beteiligten Exponenten. Insgesamt ergibt sich daher \(O(P^2\operatorname{polylog} U)\) mit \(P=168\) und \(U=\max p^5\), bei \(O(P)\) Speicher plus \(O(\log U)\) Rekursionstiefe.
Burada incelenen nicelik, aralarında asal pozitif tamsayılar \(a\le b\) için Bézout ilişkilerinden türetilir. Project Euler örneğinde girdiler Mersenne sayılarıdır: \(M_u=2^u-1\) ve \(M_v=2^v-1\). Burada \(u=r^5\) ve \(v=s^5\) olacak şekilde \(r\lt s\lt 1000\) koşulunu sağlayan farklı asallar kullanılır. Amaç, bu tüm asal çiftleri için katkıyı hesaplayıp sonucu \(10^9+7\) modunda toplamaktır.
\(2^u-1\) ve \(2^v-1\) sayılarını doğrudan kurmak mümkün değildir; çünkü üsler aşırı büyüktür. Asıl fikir, Mersenne sayıları üzerindeki Euklid algoritmasının yalnızca üsler üzerinde simüle edilebilmesidir. Gerekli katsayılar da sadece nihai modda izlendiği için dev tamsayılara hiç girilmez.
\(M_n=2^n-1\) yazalım. Aralarında asal iki pozitif tamsayı \(a\le b\) için, \(a\)'nın \(b\) modundaki en küçük pozitif tersini \(x\) ile gösterelim:
$$a x \equiv 1 \pmod{b},\qquad 1\le x \lt b.$$
O zaman
$$y=\frac{a x - 1}{b}$$
pozitif bir tamsayıdır ve \(a x - b y = 1\) eşitliğini sağlar. Ters işaretli denklemin pozitif tamamlayıcı çözümü \((b-x,\ a-y)\) olduğundan, uygulamalar şu niceliği hesaplar:
$$P(a,b)=2\min\left(x+y-1,\ (b-x)+(a-y)-1\right).$$
Dolayısıyla asıl iş, çok büyük \(u\) ve \(v\) için \(P(M_u,M_v)\) değerini hızlı bulmaktır.
\(a=M_u\) ve \(b=M_v\) olsun, ayrıca \(u\lt v\) kabul edelim. Eğer \(u\) ile \(v\) aralarında asalsa, klasik özdeşlik
$$\gcd(M_u,M_v)=M_{\gcd(u,v)}$$
bize hemen
$$\gcd(M_u,M_v)=M_1=1$$
sonucunu verir. Burada tam olarak bu durum vardır; çünkü farklı asal beşinci kuvvetleri yine aralarında asaldır:
$$\gcd(r^5,s^5)=1\qquad (r\ne s).$$
Yani hedefteki her giriş çifti Bézout düzeni için uygundur, fakat sayılar yine de doğrudan oluşturulamayacak kadar büyüktür.
\(A=qB+r\) ve \(0\le r \lt B\) olsun. O zaman
$$2^A-1=2^r\left(2^{qB}-1\right)+\left(2^r-1\right).$$
Ayrıca \(2^{qB}-1=(2^B-1)\sum_{j=0}^{q-1}2^{jB}\) olduğundan
$$M_A = Q\,M_B + M_r,$$
burada bölüm katsayısı
$$Q = 2^r\sum_{j=0}^{q-1}2^{jB}$$
şeklindedir. Demek ki \(M_A\)'yı \(M_B\)'ye böldüğümüzde kalan tam olarak \(M_r\) olur. Yani \((M_A,M_B)\) için kalan dizisi, \((A,B)\) üs çifti üzerinde sıradan Euklid algoritmasını çalıştırarak elde edilir. Problemi yapılabilir kılan sıkıştırma tam olarak budur.
Her Euklid adımında gerekli olan şey, \(Q\) bölümünün yalnızca \(10^9+7\) modundaki değeridir. Yukarıdaki formüle göre ihtiyaç duyulan miktar
$$Q \equiv 2^r\sum_{j=0}^{q-1}\left(2^B\right)^j \pmod{10^9+7}$$
olur. Uygulama bu geometrik toplamı böl-ve-yönet yaklaşımıyla hesaplar; aynı anda bir kuvvet ve ona karşılık gelen kısmi toplam izlenir. Böylece her güncelleme logaritmik derinlikte ve tamamen modüler aritmetikle yapılır.
Bu modüler bölümler, genişletilmiş Euklid yinelemesine beslenerek bir Mersenne kalanın başka bir Mersenne sayısı modundaki tersini, fakat sadece nihai mod içinde, üretir.
Şimdi
$$w = v \bmod u$$
olsun. Ayrıca \(t\), \(w\)'nun \(u\) modundaki tersi olsun:
$$w t \equiv 1 \pmod{u},\qquad 1\le t \lt u.$$
Şu geometrik toplamı ele alalım:
$$R = 1 + 2^w + 2^{2w} + \cdots + 2^{(t-1)w}.$$
O zaman
$$\begin{aligned} (2^w-1)R &= 2^{tw}-1 \\ &\equiv 2^1-1 \equiv 1 \pmod{2^u-1}, \end{aligned}$$
dolayısıyla \(R\), \(M_w\)'nun \(M_u\) modundaki tersidir. Üstelik \(M_v\equiv M_w \pmod{M_u}\) olduğundan aynı değer \(M_v\)'nin \(M_u\) modundaki tersi gibi de davranır.
Ters işaretli Bézout dalı, \(R\) yerine \(M_u-R\) alınarak elde edilir. Uygulama bu iki dal arasında
$$2t\le u$$
koşuluyla seçim yapar. Başka bir deyişle, \(t\) terimli geometrik toplam ile onun \(u-t\) terimli tümleyeni karşılaştırılır.
Bir önceki adımda seçilen katsayıya \(d\) diyelim. O zaman işarete bağlı olarak şu iki bağıntıdan biri geçerlidir:
$$M_v d = M_u x \pm 1.$$
\(x\) bulunduğunda katkı
$$P(M_u,M_v)=2(x+d-1)$$
şeklinde yazılır. Bu, genel \(P(a,b)\) formülündeki aynı minimumdur; fark sadece Mersenne tersi hesabının doğal olarak \(M_v\)'ye bağlı katsayıyı üretmesidir.
\(u=3\) ve \(v=5\) alınsın. O halde
$$M_3=7,\qquad M_5=31,\qquad w=5\bmod 3=2.$$
\(w=2\)'nin \(u=3\) modundaki tersi \(t=2\)'dir; çünkü
$$2\cdot 2 \equiv 1 \pmod{3}.$$
Bu nedenle
$$R = 1 + 2^2 = 5,$$
ve gerçekten
$$M_2 R = 3\cdot 5 = 15 \equiv 1 \pmod{7}.$$
\(2t=4>3\) olduğu için daha küçük dal tümleyeni kullanır:
$$d = M_3 - R = 7-5=2.$$
Buradan
$$x=\frac{M_5 d + 1}{M_3}=\frac{31\cdot 2 + 1}{7}=9$$
ve dolayısıyla
$$P(M_3,M_5)=2(x+d-1)=2(9+2-1)=20$$
elde edilir. Bu değer, uygulamadaki küçük doğrulama örneklerinden biriyle aynıdır.
\(1000\)'den küçük asallar \(p_1,p_2,\dots,p_k\) olsun. Nihai değer
$$\sum_{1\le i \lt j\le k} P\left(M_{p_i^5},M_{p_j^5}\right)\pmod{10^9+7}$$
biçimindedir. Yani programın matematiksel görevi, her asal çiftinin katkısını üs-düzeyli Euklid ile bulup bunları nihai mod altında toplamaktır.
C++, Python ve Java uygulamalarının ürettiği sonuç aynıdır; fakat gerçek sayısal işin tamamı C++ uygulamasında yapılır. Başlangıçta birkaç küçük kontrol çalıştırılır: önce sıradan aralarında asal tam sayılar için, sonra da genel Bézout formülüyle doğrudan karşılaştırma yapılabilen küçük Mersenne üsleri için.
Daha sonra uygulama \(1000\)'den küçük tüm asalları üretir, her birini beşinci kuvvete yükseltir ve her sırasız çifti dolaşır. Her \((u,v)\) çifti için \(M_u \bmod 10^9+7\) ve \(M_v \bmod 10^9+7\) hesaplanır, üs çifti üzerinde Euklid algoritması çalıştırılır ve gerekli bölüm katsayıları yalnızca \(10^9+7\) modunda geri kurulur.
Asıl hızlandırıcı, geometrik toplamların böl-ve-yönet yöntemiyle hesaplanmasıdır. Uygulama iki dev Mersenne sayısının gerçek bölümünü oluşturmaz; bunun yerine bu bölümün \(10^9+7\) modundaki değerini doğrudan üretir. Genişletilmiş Euklid yinelemesinin ihtiyacı da tam olarak budur. Ardından \(w^{-1}\bmod u\) bilgisinden yararlanan dal seçimi, hangi Bézout işaretinin daha küçük katkı vereceğini belirler.
Son aşamada her asal çiftinin katkısı akan toplama \(10^9+7\) modunda eklenir. Python ve Java uygulamaları ise aynı derlenmiş hesabı başlatan ince sarmalayıcılardır; bu nedenle üç dilde de hem matematik hem çıktı birebir uyumludur.
\(1000\)'den küçük \(168\) asal vardır; dolayısıyla dış döngü \(\binom{168}{2}=14028\) sırasız çifti işler. Tek bir çift için üs düzeyindeki Euklid algoritması, sıradan Euklid gibi \(O(\log \max(u,v))\) bölme adımı ister. Her adım sabit sayıda modüler üs alma ve bir böl-ve-yönet geometrik toplam hesabı yapar; bunların ikisi de ilgili üs büyüklüklerinde logaritmiktir. Sonuç olarak toplam zaman \(P=168\) ve \(U=\max p^5\) olmak üzere \(O(P^2\operatorname{polylog} U)\), bellek ise asal listesi için \(O(P)\) ve özyineleme derinliği için \(O(\log U)\) mertebesindedir.
Para enteros positivos coprimos \(a\le b\), la cantidad estudiada aquí se construye a partir de relaciones de Bézout entre \(a\) y \(b\). En el caso de Project Euler, las entradas son números de Mersenne \(M_u=2^u-1\) y \(M_v=2^v-1\), donde \(u=r^5\) y \(v=s^5\) para primos distintos \(r\lt s\lt 1000\). El objetivo es evaluar la contribución de cada uno de esos pares de primos y sumar todo módulo \(10^9+7\).
Construir directamente \(2^u-1\) y \(2^v-1\) es inviable porque los exponentes son enormes. La observación decisiva es que el algoritmo de Euclides sobre números de Mersenne puede simularse por completo sobre los exponentes, mientras todos los coeficientes relevantes se mantienen solo módulo del valor final.
Escribamos \(M_n=2^n-1\). Para dos enteros positivos coprimos \(a\le b\), sea \(x\) el menor inverso positivo de \(a\) módulo \(b\):
$$a x \equiv 1 \pmod{b},\qquad 1\le x \lt b.$$
Entonces
$$y=\frac{a x - 1}{b}$$
es un entero positivo y satisface \(a x - b y = 1\). La solución positiva complementaria para la ecuación con el signo opuesto es \((b-x,\ a-y)\), de modo que las implementaciones evalúan
$$P(a,b)=2\min\left(x+y-1,\ (b-x)+(a-y)-1\right).$$
Así, el problema pasa a ser calcular \(P(M_u,M_v)\) rápidamente para exponentes muy grandes.
Tomemos \(a=M_u\) y \(b=M_v\) con \(u\lt v\). Si \(u\) y \(v\) son coprimos, la identidad clásica
$$\gcd(M_u,M_v)=M_{\gcd(u,v)}$$
implica
$$\gcd(M_u,M_v)=M_1=1.$$
Eso es exactamente lo que ocurre aquí, porque quintas potencias de primos distintos siguen siendo coprimas:
$$\gcd(r^5,s^5)=1\qquad (r\ne s).$$
Por tanto, cada par objetivo es válido para la construcción de Bézout, aunque los propios números sean enormes.
Si \(A=qB+r\) con \(0\le r \lt B\), entonces
$$2^A-1=2^r\left(2^{qB}-1\right)+\left(2^r-1\right).$$
Usando \(2^{qB}-1=(2^B-1)\sum_{j=0}^{q-1}2^{jB}\), esto se reescribe como
$$M_A = Q\,M_B + M_r,$$
con cociente
$$Q = 2^r\sum_{j=0}^{q-1}2^{jB}.$$
Así, el resto al dividir \(M_A\) por \(M_B\) es exactamente \(M_r\). En consecuencia, toda la sucesión de restos para \((M_A,M_B)\) se obtiene ejecutando el algoritmo de Euclides ordinario sobre el par de exponentes \((A,B)\). Esa es la compresión esencial del método.
En cada paso solo hace falta el cociente \(Q\) módulo \(10^9+7\). Por la fórmula anterior, basta con calcular
$$Q \equiv 2^r\sum_{j=0}^{q-1}\left(2^B\right)^j \pmod{10^9+7}.$$
La implementación evalúa esta suma geométrica mediante divide y vencerás, llevando al mismo tiempo una potencia y la suma parcial asociada. De este modo, cada actualización del cociente se hace con aritmética modular y profundidad logarítmica, sin construir enteros monstruosos.
Al introducir esos cocientes modulares en la recurrencia del algoritmo extendido de Euclides, se obtiene el inverso de un resto de Mersenne módulo otro número de Mersenne, pero expresado solo módulo \(10^9+7\).
Sea
$$w = v \bmod u,$$
y sea \(t\) el inverso de \(w\) módulo \(u\):
$$w t \equiv 1 \pmod{u},\qquad 1\le t \lt u.$$
Consideremos ahora
$$R = 1 + 2^w + 2^{2w} + \cdots + 2^{(t-1)w}.$$
Entonces
$$\begin{aligned} (2^w-1)R &= 2^{tw}-1 \\ &\equiv 2^1-1 \equiv 1 \pmod{2^u-1}, \end{aligned}$$
de modo que \(R\) es un inverso de \(M_w\) módulo \(M_u\). Como además \(M_v\equiv M_w \pmod{M_u}\), el mismo valor también actúa como inverso de \(M_v\) módulo \(M_u\).
La rama de Bézout con signo opuesto se obtiene sustituyendo \(R\) por \(M_u-R\). La implementación elige entre ambas ramas mediante el criterio
$$2t\le u,$$
es decir, comparando la suma geométrica de \(t\) términos con su complemento de \(u-t\) términos dentro de \(M_u=1+2+\cdots+2^{u-1}\).
Sea \(d\) el coeficiente elegido en el paso anterior. Entonces, según el signo seleccionado, vale una de las dos relaciones
$$M_v d = M_u x \pm 1.$$
Una vez recuperado \(x\), la contribución queda
$$P(M_u,M_v)=2(x+d-1).$$
Este es exactamente el mismo mínimo que aparecía en la fórmula general de \(P(a,b)\); la diferencia es solo que aquí se trabaja con el coeficiente asociado a \(M_v\), porque es el que sale de manera natural al calcular el inverso de Mersenne.
Tomemos \(u=3\) y \(v=5\). Entonces
$$M_3=7,\qquad M_5=31,\qquad w=5\bmod 3=2.$$
El inverso de \(w=2\) módulo \(u=3\) es \(t=2\), porque
$$2\cdot 2 \equiv 1 \pmod{3}.$$
Por tanto,
$$R = 1 + 2^2 = 5,$$
y efectivamente
$$M_2 R = 3\cdot 5 = 15 \equiv 1 \pmod{7}.$$
Como \(2t=4>3\), la rama menor usa el complemento
$$d = M_3 - R = 7-5=2.$$
Ahora
$$x=\frac{M_5 d + 1}{M_3}=\frac{31\cdot 2 + 1}{7}=9,$$
de donde
$$P(M_3,M_5)=2(x+d-1)=2(9+2-1)=20.$$
Ese valor coincide con una de las verificaciones pequeñas incluidas en la implementación.
Si \(p_1,p_2,\dots,p_k\) son los primos menores que \(1000\), el resultado final es
$$\sum_{1\le i \lt j\le k} P\left(M_{p_i^5},M_{p_j^5}\right)\pmod{10^9+7}.$$
Por lo tanto, la tarea matemática del programa consiste en evaluar cada contribución de par mediante Euclides sobre exponentes y acumular todo bajo el módulo primo final.
Las implementaciones en C++, Python y Java producen el mismo resultado, pero el trabajo aritmético real se realiza en la implementación de C++. Primero se ejecutan varios controles pequeños: al principio con pares coprimos ordinarios y luego con exponentes de Mersenne pequeños, donde todavía es posible comparar directamente contra la fórmula general de Bézout.
Después, la implementación genera todos los primos menores que \(1000\), eleva cada uno a la quinta potencia y recorre cada par no ordenado. Para cada \((u,v)\), calcula \(M_u \bmod 10^9+7\) y \(M_v \bmod 10^9+7\), ejecuta el algoritmo de Euclides sobre el par de exponentes y reconstruye la información necesaria de los cocientes solo módulo \(10^9+7\).
El acelerador clave es la evaluación divide y vencerás de las sumas geométricas. En lugar de formar el cociente real entre dos enormes números de Mersenne, la implementación produce directamente ese cociente módulo \(10^9+7\), que es justo lo que necesita la recurrencia del Euclides extendido. A continuación, la decisión de rama basada en el inverso de \(w\) módulo \(u\) determina qué signo de Bézout da la contribución menor.
Por último, cada contribución de par primo se añade a una suma acumulada módulo \(10^9+7\). Las implementaciones en Python y Java son lanzadores ligeros alrededor de ese mismo cálculo compilado, así que tanto la matemática como la salida final permanecen sincronizadas en los tres lenguajes.
Hay \(168\) primos por debajo de \(1000\), de modo que el bucle exterior procesa \(\binom{168}{2}=14028\) pares no ordenados. Para un par, el algoritmo de Euclides sobre exponentes realiza \(O(\log \max(u,v))\) divisiones, igual que Euclides ordinario sobre enteros. Cada división hace un número constante de exponenciaciones modulares y una evaluación divide y vencerás de una suma geométrica, ambas logarítmicas en los tamaños de exponente relevantes. En conjunto, el tiempo total es \(O(P^2\operatorname{polylog} U)\), con \(P=168\) y \(U=\max p^5\), mientras que la memoria es \(O(P)\) para la lista de primos más una profundidad recursiva de \(O(\log U)\).
这里要计算的量，来自互素正整数 \(a\le b\) 的 Bézout 关系。在 Project Euler 的具体版本中，输入不是普通整数，而是 Mersenne 数 \(M_u=2^u-1\) 与 \(M_v=2^v-1\)。其中 \(u=r^5\)、\(v=s^5\)，而 \(r\lt s\lt 1000\) 是不同的素数。目标是对所有这样的素数对分别求出贡献，并把总和对 \(10^9+7\) 取模。
由于指数极大，直接构造 \(2^u-1\) 与 \(2^v-1\) 完全不可行。真正的突破点是：Mersenne 数上的欧几里得算法，可以在指数层面完整模拟；与此同时，所有真正需要保留的系数，只要在最终模数下维护即可。
记 \(M_n=2^n-1\)。对两个互素正整数 \(a\le b\)，令 \(x\) 为 \(a\) 在模 \(b\) 意义下最小的正逆元：
$$a x \equiv 1 \pmod{b},\qquad 1\le x \lt b.$$
则
$$y=\frac{a x - 1}{b}$$
是一个正整数，并满足 \(a x - b y = 1\)。而符号相反的那条 Bézout 分支，其对应的正解正好是 \((b-x,\ a-y)\)。因此，程序实际计算的是
$$P(a,b)=2\min\left(x+y-1,\ (b-x)+(a-y)-1\right).$$
所以核心问题就变成：如何在 \(u\)、\(v\) 极大的情况下，高效计算 \(P(M_u,M_v)\)。
设 \(a=M_u\)、\(b=M_v\)，并且 \(u\lt v\)。如果 \(u\) 与 \(v\) 互素，那么由经典恒等式
$$\gcd(M_u,M_v)=M_{\gcd(u,v)}$$
立刻得到
$$\gcd(M_u,M_v)=M_1=1.$$
本题正好属于这种情形，因为不同素数的五次幂仍然互素：
$$\gcd(r^5,s^5)=1\qquad (r\ne s).$$
这说明所有目标输入对都满足 Bézout 结构的前提，只是对应的 Mersenne 数本身大到无法显式写出。
若 \(A=qB+r\)，且 \(0\le r \lt B\)，那么
$$2^A-1=2^r\left(2^{qB}-1\right)+\left(2^r-1\right).$$
再把 \(2^{qB}-1\) 写成 \((2^B-1)\sum_{j=0}^{q-1}2^{jB}\)，就得到
$$M_A = Q\,M_B + M_r,$$
其中商为
$$Q = 2^r\sum_{j=0}^{q-1}2^{jB}.$$
因此，用 \(M_B\) 去除 \(M_A\) 时，余数恰好就是 \(M_r\)。也就是说，\((M_A,M_B)\) 的整个余数序列，和普通整数对 \((A,B)\) 的欧几里得算法余数序列完全对应。这个观察把问题从“巨大整数”压缩成了“中等大小的指数”。
在每一步欧几里得递推中，真正需要的只是商 \(Q\) 在模 \(10^9+7\) 下的值。由上式可知，只需求
$$Q \equiv 2^r\sum_{j=0}^{q-1}\left(2^B\right)^j \pmod{10^9+7}.$$
实现采用分治方式计算这类几何级数：同时维护一个幂值以及对应的部分和。这样一来，每一步都只涉及模运算与对数级递归深度，而不会真的去构造两个巨大 Mersenne 数相除得到的完整商。
把这些“只在最终模数下表示”的商送入扩展欧几里得递推后，就能得到某个 Mersenne 余数在另一个 Mersenne 数模下的逆元，不过整个过程始终停留在 \(10^9+7\) 这个模数之内。
令
$$w = v \bmod u,$$
再令 \(t\) 表示 \(w\) 在模 \(u\) 下的逆元：
$$w t \equiv 1 \pmod{u},\qquad 1\le t \lt u.$$
考虑几何和
$$R = 1 + 2^w + 2^{2w} + \cdots + 2^{(t-1)w}.$$
则有
$$\begin{aligned} (2^w-1)R &= 2^{tw}-1 \\ &\equiv 2^1-1 \equiv 1 \pmod{2^u-1}, \end{aligned}$$
因此 \(R\) 就是 \(M_w\) 在模 \(M_u\) 下的逆元。又因为 \(M_v\equiv M_w \pmod{M_u}\)，所以同一个 \(R\) 也可以看作 \(M_v\) 在模 \(M_u\) 下的逆元。
另一条符号相反的 Bézout 分支，则对应把 \(R\) 换成 \(M_u-R\)。实现用非常简单的判据
$$2t\le u$$
在这两个候选之间做选择。直观上，这等价于比较一个含有 \(t\) 项的几何和，与它在 \(M_u=1+2+\cdots+2^{u-1}\) 中的 \(u-t\) 项补集，哪一边更短。
设上一步选中的系数为 \(d\)。那么根据选择的是哪一个符号分支，以下两式必有其一成立：
$$M_v d = M_u x \pm 1.$$
一旦恢复出 \(x\)，对应贡献就是
$$P(M_u,M_v)=2(x+d-1).$$
这与开头给出的通用 \(P(a,b)\) 公式是同一个最小值，只不过在 Mersenne 情形下，程序自然先得到的是附着在 \(M_v\) 上的那个系数，因此表达式写成了 \(x+d-1\) 的样子。
取 \(u=3\)、\(v=5\)。此时
$$M_3=7,\qquad M_5=31,\qquad w=5\bmod 3=2.$$
\(w=2\) 在模 \(u=3\) 下的逆元是 \(t=2\)，因为
$$2\cdot 2 \equiv 1 \pmod{3}.$$
于是
$$R = 1 + 2^2 = 5,$$
并且确实有
$$M_2 R = 3\cdot 5 = 15 \equiv 1 \pmod{7}.$$
由于 \(2t=4>3\)，较小的那条分支选用补系数
$$d = M_3 - R = 7-5=2.$$
接着
$$x=\frac{M_5 d + 1}{M_3}=\frac{31\cdot 2 + 1}{7}=9,$$
所以
$$P(M_3,M_5)=2(x+d-1)=2(9+2-1)=20.$$
这个结果正好与实现中使用的小规模校验值一致。
设所有小于 \(1000\) 的素数为 \(p_1,p_2,\dots,p_k\)。那么最终需要计算的是
$$\sum_{1\le i \lt j\le k} P\left(M_{p_i^5},M_{p_j^5}\right)\pmod{10^9+7}.$$
因此，整道题的数学核心可以概括为：对每一对素数五次幂指数，利用指数层面的欧几里得算法求出对应的 Bézout 贡献，再在最终模数下把它们加总。
C++、Python 和 Java 三个实现对外表现一致，但真正执行数论计算的是 C++ 实现。它一开始先跑若干个小型校验：先验证普通互素整数对的通用 Bézout 公式，再验证若干小指数 Mersenne 情形，确保“指数层算法”和“直接算小整数”得到同样的结果。
随后，实现生成所有小于 \(1000\) 的素数，把每个素数提升到五次幂，并遍历每一个无序素数对。对于每个 \((u,v)\)，程序先算出 \(M_u \bmod 10^9+7\) 和 \(M_v \bmod 10^9+7\)，然后在指数对 \((u,v)\) 上运行欧几里得算法，把每一步所需的商都只以“模 \(10^9+7\) 的值”恢复出来。
实现中的关键加速点，是对几何级数进行分治求值。程序不会真正构造两个巨大 Mersenne 数相除的完整商，而是直接产出这个商在最终模数下的值；而扩展欧几里得递推所需要的恰恰就是这个量。之后，再结合 \(w\) 在模 \(u\) 下的逆元，决定两条 Bézout 符号分支中哪一条给出更小的贡献。
最后，每个素数对的贡献都被加入到一个模 \(10^9+7\) 的滚动总和中。Python 与 Java 的实现则是对同一份编译后核心计算的轻量封装，因此三种语言在数学逻辑和最终输出上保持完全一致。
小于 \(1000\) 的素数共有 \(168\) 个，所以外层总共处理 \(\binom{168}{2}=14028\) 个无序配对。对任意一个配对，指数层面的欧几里得算法需要 \(O(\log \max(u,v))\) 次除法，这与普通整数欧几里得算法相同。每一次除法会触发常数次模幂运算，以及一次分治式几何和计算，这两部分对相关指数规模都是对数级的。因此总时间可以记为 \(O(P^2\operatorname{polylog} U)\)，其中 \(P=168\)、\(U=\max p^5\)；空间复杂度为素数表的 \(O(P)\)，外加递归所需的 \(O(\log U)\) 深度。
Здесь рассматривается величина, возникающая из соотношений Безу для взаимно простых положительных целых чисел \(a\le b\). В варианте Project Euler входами являются числа Мерсенна \(M_u=2^u-1\) и \(M_v=2^v-1\), где \(u=r^5\) и \(v=s^5\) для различных простых \(r\lt s\lt 1000\). Нужно вычислить вклад каждого такого простого пары и просуммировать все вклады по модулю \(10^9+7\).
Непосредственно строить числа \(2^u-1\) и \(2^v-1\) невозможно: показатели слишком велики. Ключевая идея состоит в том, что алгоритм Евклида для чисел Мерсенна можно целиком перенести на уровень показателей, а все действительно нужные коэффициенты хранить только по модулю конечного модуля.
Обозначим \(M_n=2^n-1\). Для двух взаимно простых положительных целых \(a\le b\) пусть \(x\) — наименьший положительный обратный элемент к \(a\) по модулю \(b\):
$$a x \equiv 1 \pmod{b},\qquad 1\le x \lt b.$$
Тогда
$$y=\frac{a x - 1}{b}$$
является положительным целым числом и удовлетворяет равенству \(a x - b y = 1\). Комплементарное положительное решение для уравнения с противоположным знаком равно \((b-x,\ a-y)\), поэтому реализации вычисляют
$$P(a,b)=2\min\left(x+y-1,\ (b-x)+(a-y)-1\right).$$
Следовательно, главная задача — быстро находить \(P(M_u,M_v)\) для очень больших показателей \(u\) и \(v\).
Положим \(a=M_u\) и \(b=M_v\), где \(u\lt v\). Если \(u\) и \(v\) взаимно просты, то классическая формула
$$\gcd(M_u,M_v)=M_{\gcd(u,v)}$$
немедленно даёт
$$\gcd(M_u,M_v)=M_1=1.$$
Именно это и происходит в задаче, потому что пятые степени различных простых остаются взаимно простыми:
$$\gcd(r^5,s^5)=1\qquad (r\ne s).$$
Значит, каждая целевая пара подходит под схему Безу, хотя сами числа слишком велики для явного представления.
Если \(A=qB+r\) и \(0\le r \lt B\), то
$$2^A-1=2^r\left(2^{qB}-1\right)+\left(2^r-1\right).$$
После разложения \(2^{qB}-1=(2^B-1)\sum_{j=0}^{q-1}2^{jB}\) получаем
$$M_A = Q\,M_B + M_r,$$
где частное равно
$$Q = 2^r\sum_{j=0}^{q-1}2^{jB}.$$
Таким образом, остаток при делении \(M_A\) на \(M_B\) в точности равен \(M_r\). Следовательно, вся последовательность остатков для \((M_A,M_B)\) получается обычным алгоритмом Евклида, применённым к паре показателей \((A,B)\). Именно это сжатие делает задачу вычислимой.
На каждом шаге нужен только остаток частного \(Q\) по модулю \(10^9+7\). Из формулы выше следует, что достаточно найти
$$Q \equiv 2^r\sum_{j=0}^{q-1}\left(2^B\right)^j \pmod{10^9+7}.$$
Реализация вычисляет эту геометрическую сумму методом «разделяй и властвуй», одновременно отслеживая степень и соответствующую частичную сумму. Поэтому каждый шаг обновления частного требует только модульной арифметики и логарифмической глубины рекурсии, а не работы с чудовищно большими целыми числами.
Если подать такие модульные частные в рекуррентные формулы расширенного алгоритма Евклида, то можно получить обратный элемент для одного остатка Мерсенна по модулю другого числа Мерсенна, но уже полностью в представлении по модулю \(10^9+7\).
Пусть
$$w = v \bmod u,$$
а \(t\) — обратный к \(w\) по модулю \(u\):
$$w t \equiv 1 \pmod{u},\qquad 1\le t \lt u.$$
Рассмотрим сумму
$$R = 1 + 2^w + 2^{2w} + \cdots + 2^{(t-1)w}.$$
Тогда
$$\begin{aligned} (2^w-1)R &= 2^{tw}-1 \\ &\equiv 2^1-1 \equiv 1 \pmod{2^u-1}, \end{aligned}$$
поэтому \(R\) является обратным элементом к \(M_w\) по модулю \(M_u\). Так как \(M_v\equiv M_w \pmod{M_u}\), тот же самый \(R\) работает и как обратный элемент к \(M_v\) по модулю \(M_u\).
Вторая ветвь Безу с противоположным знаком получается заменой \(R\) на \(M_u-R\). Реализация выбирает между двумя ветвями по простому правилу
$$2t\le u,$$
то есть фактически сравнивает геометрическую сумму из \(t\) членов с её дополнением из \(u-t\) членов внутри \(M_u=1+2+\cdots+2^{u-1}\).
Пусть \(d\) — коэффициент, выбранный на предыдущем шаге. Тогда в зависимости от знака выполняется одно из двух равенств
$$M_v d = M_u x \pm 1.$$
После восстановления \(x\) вклад имеет вид
$$P(M_u,M_v)=2(x+d-1).$$
Это тот же самый минимум, что и в общей формуле для \(P(a,b)\); отличие только в том, что в случае Мерсенна естественным образом сначала вычисляется коэффициент при \(M_v\).
Возьмём \(u=3\) и \(v=5\). Тогда
$$M_3=7,\qquad M_5=31,\qquad w=5\bmod 3=2.$$
Обратный к \(w=2\) по модулю \(u=3\) равен \(t=2\), поскольку
$$2\cdot 2 \equiv 1 \pmod{3}.$$
Значит,
$$R = 1 + 2^2 = 5,$$
и действительно
$$M_2 R = 3\cdot 5 = 15 \equiv 1 \pmod{7}.$$
Так как \(2t=4>3\), меньшая ветвь использует дополнение
$$d = M_3 - R = 7-5=2.$$
Теперь
$$x=\frac{M_5 d + 1}{M_3}=\frac{31\cdot 2 + 1}{7}=9,$$
следовательно,
$$P(M_3,M_5)=2(x+d-1)=2(9+2-1)=20.$$
Именно этот результат совпадает с одним из малых контрольных примеров в реализации.
Если \(p_1,p_2,\dots,p_k\) — все простые числа меньше \(1000\), то окончательный ответ равен
$$\sum_{1\le i \lt j\le k} P\left(M_{p_i^5},M_{p_j^5}\right)\pmod{10^9+7}.$$
Значит, математическая сущность программы состоит в том, чтобы для каждой пары пятых степеней простых вычислить вклад через алгоритм Евклида на показателях, а затем сложить все эти вклады по модулю конечного простого числа.
Реализации на C++, Python и Java выдают один и тот же результат, но всё реальное числовое вычисление находится в C++-реализации. В начале она выполняет несколько небольших проверок: сначала на обычных взаимно простых парах целых чисел, затем на малых показателях Мерсенна, где ещё можно напрямую сравнить ответ с общей формулой Безу.
После этого реализация строит все простые меньше \(1000\), возводит каждое из них в пятую степень и перебирает все неупорядоченные пары. Для каждой пары \((u,v)\) вычисляются \(M_u \bmod 10^9+7\) и \(M_v \bmod 10^9+7\), запускается алгоритм Евклида на паре показателей и восстанавливается нужная информация о частных только по модулю \(10^9+7\).
Главное ускорение даёт вычисление геометрических сумм по схеме «разделяй и властвуй». Вместо того чтобы строить настоящее частное при делении двух гигантских чисел Мерсенна, реализация сразу получает это частное по модулю \(10^9+7\), а именно это и нужно для рекуррентных формул расширенного алгоритма Евклида. Затем выбор ветви, основанный на обратном элементе к \(w\) по модулю \(u\), определяет, какой знак Безу даёт меньший вклад.
Наконец, вклад каждой простой пары добавляется в текущую сумму по модулю \(10^9+7\). Реализации на Python и Java представляют собой тонкие оболочки вокруг той же самой скомпилированной вычислительной части, поэтому математика и итоговый ответ полностью совпадают во всех трёх языках.
Простых чисел меньше \(1000\) ровно \(168\), поэтому внешний цикл обрабатывает \(\binom{168}{2}=14028\) неупорядоченных пар. Для одной пары алгоритм Евклида на показателях делает \(O(\log \max(u,v))\) делений, как и обычный алгоритм Евклида на целых числах. Каждое деление вызывает постоянное число модульных возведений в степень и одно вычисление геометрической суммы по схеме «разделяй и властвуй»; обе операции логарифмичны по соответствующим размерам показателей. Итак, суммарное время равно \(O(P^2\operatorname{polylog} U)\), где \(P=168\), \(U=\max p^5\), а память составляет \(O(P)\) для списка простых плюс \(O(\log U)\) глубины рекурсии.
الكمية المدروسة هنا تُبنى من علاقات بيزو بين عددين صحيحين موجبين ومتحابين نسبيًا \(a\le b\). في نسخة Project Euler تكون المدخلات من أعداد ميرسين \(M_u=2^u-1\) و \(M_v=2^v-1\)، حيث \(u=r^5\) و \(v=s^5\) لعددين أوليين مختلفين \(r\lt s\lt 1000\). المطلوب هو حساب مساهمة كل زوج أولي من هذا النوع ثم جمع جميع المساهمات بترديد \(10^9+7\).
لا يمكن بناء \(2^u-1\) و \(2^v-1\) مباشرة لأن الأسس هائلة جدًا. الفكرة الحاسمة هي أن خوارزمية إقليدس على أعداد ميرسين يمكن محاكاتها بالكامل على مستوى الأسس فقط، بينما تُتَابَع كل المعاملات المهمة تحت الترديد النهائي لا غير.
لنكتب \(M_n=2^n-1\). إذا كان \(a\le b\) عددين صحيحين موجبين ومتحابين نسبيًا، فليكن \(x\) أصغر معكوس موجب لـ \(a\) بترديد \(b\):
$$a x \equiv 1 \pmod{b},\qquad 1\le x \lt b.$$
عندئذ يكون
$$y=\frac{a x - 1}{b}$$
عددًا صحيحًا موجبًا ويحقق \(a x - b y = 1\). أما الحل الموجب المكمّل للمعادلة ذات الإشارة المعاكسة فهو \((b-x,\ a-y)\)، ولذلك تحسب التطبيقات القيمة
$$P(a,b)=2\min\left(x+y-1,\ (b-x)+(a-y)-1\right).$$
إذن لبّ المسألة هو حساب \(P(M_u,M_v)\) بكفاءة عندما يكون \(u\) و \(v\) كبيرين جدًا.
ضع \(a=M_u\) و \(b=M_v\) مع \(u\lt v\). إذا كان \(u\) و \(v\) متحابين نسبيًا، فإن الهوية الكلاسيكية
$$\gcd(M_u,M_v)=M_{\gcd(u,v)}$$
تعطي فورًا
$$\gcd(M_u,M_v)=M_1=1.$$
وهذا هو الوضع هنا بالضبط، لأن القوى الخامسة لعددين أوليين مختلفين تبقى متحابة نسبيًا:
$$\gcd(r^5,s^5)=1\qquad (r\ne s).$$
إذًا كل زوج مستهدف يحقق شرط بيزو، حتى لو كانت الأعداد نفسها فلكية الحجم.
إذا كان \(A=qB+r\) مع \(0\le r \lt B\)، فإن
$$2^A-1=2^r\left(2^{qB}-1\right)+\left(2^r-1\right).$$
وبعد كتابة \(2^{qB}-1=(2^B-1)\sum_{j=0}^{q-1}2^{jB}\) نحصل على
$$M_A = Q\,M_B + M_r,$$
حيث خارج القسمة هو
$$Q = 2^r\sum_{j=0}^{q-1}2^{jB}.$$
إذن الباقي عند قسمة \(M_A\) على \(M_B\) هو بالضبط \(M_r\). وهذا يعني أن تسلسل البواقي الكامل للزوج \((M_A,M_B)\) يُستخرج بمجرد تشغيل خوارزمية إقليدس العادية على زوج الأسس \((A,B)\). هذا هو الاختزال الأساسي الذي يجعل المسألة ممكنة حسابيًا.
في كل خطوة لا نحتاج إلا إلى قيمة \(Q\) بترديد \(10^9+7\). ومن الصيغة السابقة يكفي حساب
$$Q \equiv 2^r\sum_{j=0}^{q-1}\left(2^B\right)^j \pmod{10^9+7}.$$
تقوم الشيفرة بحساب هذه المتسلسلة الهندسية بطريقة التقسيم والسيطرة، مع تتبع قوة واحدة ومجموعها الجزئي في الوقت نفسه. وهكذا تبقى كل عملية تحديث لوخارج القسمة داخل الحسابات المعيارية وضمن عمق لوغاريتمي، من غير تكوين أعداد ميرسين الهائلة فعليًا.
وعند تغذية هذه القيم المعيارية في علاقة إقليدس الموسع، نحصل على معكوس باقي ميرسين بترديد عدد ميرسين آخر، لكن كل ذلك ممثّل فقط تحت الترديد النهائي \(10^9+7\).
ليكن
$$w = v \bmod u,$$
ولتكن \(t\) هي معكوس \(w\) بترديد \(u\):
$$w t \equiv 1 \pmod{u},\qquad 1\le t \lt u.$$
انظر الآن إلى المجموع
$$R = 1 + 2^w + 2^{2w} + \cdots + 2^{(t-1)w}.$$
عندئذ
$$\begin{aligned} (2^w-1)R &= 2^{tw}-1 \\ &\equiv 2^1-1 \equiv 1 \pmod{2^u-1}, \end{aligned}$$
ومن ثم فإن \(R\) هو معكوس \(M_w\) بترديد \(M_u\). وبما أن \(M_v\equiv M_w \pmod{M_u}\)، فإن القيمة نفسها تعمل أيضًا بوصفها معكوسًا لـ \(M_v\) بترديد \(M_u\).
أما فرع بيزو ذي الإشارة المعاكسة فينتج من استبدال \(R\) بـ \(M_u-R\). وتختار الشيفرة بين الفرعين وفق المعيار البسيط
$$2t\le u,$$
أي أنها تقارن عمليًا بين مجموع هندسي طوله \(t\) ومتمّمه ذي الطول \(u-t\) داخل \(M_u=1+2+\cdots+2^{u-1}\).
لنسمّ المعامل المختار في الخطوة السابقة \(d\). عندها تتحقق، بحسب الإشارة المختارة، إحدى العلاقتين
$$M_v d = M_u x \pm 1.$$
وبمجرد استعادة \(x\)، تصبح المساهمة
$$P(M_u,M_v)=2(x+d-1).$$
وهذه هي القيمة الصغرى نفسها التي ظهرت في الصيغة العامة لـ \(P(a,b)\)، لكن التعبير هنا مكتوب بمعامل مرتبط بـ \(M_v\) لأن هذا هو الشكل الطبيعي الذي تنتجه عملية معكوس ميرسين.
خذ \(u=3\) و \(v=5\). عندئذ
$$M_3=7,\qquad M_5=31,\qquad w=5\bmod 3=2.$$
معكوس \(w=2\) بترديد \(u=3\) هو \(t=2\)، لأن
$$2\cdot 2 \equiv 1 \pmod{3}.$$
ومن ثم
$$R = 1 + 2^2 = 5,$$
وفعلًا
$$M_2 R = 3\cdot 5 = 15 \equiv 1 \pmod{7}.$$
وبما أن \(2t=4>3\)، فإن الفرع الأصغر يستخدم المتمّم
$$d = M_3 - R = 7-5=2.$$
والآن
$$x=\frac{M_5 d + 1}{M_3}=\frac{31\cdot 2 + 1}{7}=9,$$
فتكون
$$P(M_3,M_5)=2(x+d-1)=2(9+2-1)=20.$$
وهذا يطابق إحدى حالات التحقق الصغيرة الموجودة في التنفيذ.
إذا كانت \(p_1,p_2,\dots,p_k\) هي جميع الأعداد الأولية الأصغر من \(1000\)، فإن القيمة النهائية هي
$$\sum_{1\le i \lt j\le k} P\left(M_{p_i^5},M_{p_j^5}\right)\pmod{10^9+7}.$$
وبذلك يمكن تلخيص المهمة الرياضية للبرنامج في حساب مساهمة كل زوج من أسس أولية خامسة بواسطة إقليدس على مستوى الأسس، ثم جمع كل المساهمات تحت الترديد الأولي النهائي.
تعطي تطبيقات C++ وPython وJava النتيجة نفسها، لكن العمل العددي الحقيقي ينفَّذ داخل تطبيق C++. يبدأ التنفيذ بعدة اختبارات صغيرة: أولًا على أزواج عادية من الأعداد المتحابة نسبيًا، ثم على أسس صغيرة لأعداد ميرسين، حيث يظل من الممكن مقارنة الناتج مباشرة مع صيغة بيزو العامة.
بعد ذلك تُولِّد الشيفرة جميع الأعداد الأولية الأصغر من \(1000\)، ثم ترفع كلًّا منها إلى القوة الخامسة، وتستعرض كل زوج غير مرتب. ولكل زوج \((u,v)\) تحسب \(M_u \bmod 10^9+7\) و \(M_v \bmod 10^9+7\)، ثم تشغّل خوارزمية إقليدس على زوج الأسس، وتسترجع معلومات خارج القسمة المطلوبة فقط بترديد \(10^9+7\).
أهم تسريع في التنفيذ هو حساب المتسلسلات الهندسية بطريقة التقسيم والسيطرة. فبدلًا من تكوين خارج القسمة الحقيقي بين عددين هائلين من أعداد ميرسين، تنتج الشيفرة هذا الخارج مباشرة بترديد \(10^9+7\)، وهو بالضبط ما تحتاجه علاقة إقليدس الموسع. وبعد ذلك يحدد اختبار الفرع المعتمد على معكوس \(w\) بترديد \(u\) أي إشارة من إشارات بيزو تعطي المساهمة الأصغر.
وفي النهاية تُضاف مساهمة كل زوج أولي إلى مجموع جارٍ بترديد \(10^9+7\). أما تطبيقا Python وJava فهما مجرد غلافين خفيفين حول الحساب المترجَم نفسه، ولذلك تبقى الرياضيات والمخرج النهائي متطابقين تمامًا في اللغات الثلاث.
يوجد \(168\) عددًا أوليًا أصغر من \(1000\)، ولذلك تعالج الحلقة الخارجية \(\binom{168}{2}=14028\) زوجًا غير مرتب. ولكل زوج يحتاج إقليدس على مستوى الأسس إلى \(O(\log \max(u,v))\) من القسَمات، تمامًا كما في خوارزمية إقليدس المعتادة على الأعداد الصحيحة. وكل قسمة تنفّذ عددًا ثابتًا من عمليات الأسّ المعياري وحسابًا واحدًا لمتسلسلة هندسية بطريقة التقسيم والسيطرة، وكلاهما لوغاريتمي في أحجام الأسس المعنية. لذا يكون الزمن الكلي \(O(P^2\operatorname{polylog} U)\) حيث \(P=168\) و \(U=\max p^5\)، بينما يساوي استهلاك الذاكرة \(O(P)\) لقائمة الأعداد الأولية مع عمق استدعاء \(O(\log U)\).