For each integer \(n \ge 3\), define
$$I(n)=\max\{m\in\mathbb{Z}_{>0}: m<n-1,\ m^{-1}\equiv m \pmod n\}.$$
The condition \(m^{-1}\equiv m \pmod n\) means that \(m\) is its own multiplicative inverse modulo \(n\). Multiplying both sides by \(m\) gives
$$m^2\equiv 1 \pmod n.$$
So the task is to find, for every \(n\), the largest solution of \(x^2\equiv 1 \pmod n\) that is strictly smaller than \(n-1\), and then sum those values for \(3 \le n \le 2\cdot 10^7\). The checkpoints \(I(7)=1\) and \(I(100)=51\) confirm that we are looking for the largest nontrivial square root of \(1\) modulo \(n\), because \(n-1\equiv -1\pmod n\) is always a solution but is excluded.
If \(x^2\equiv 1 \pmod n\), then \(n\mid (x-1)(x+1)\). Every solution is automatically coprime to \(n\), because a common divisor of \(x\) and \(n\) would also divide \(1\). Therefore the original inverse condition and the quadratic congruence are completely equivalent.
This reduces the problem to understanding the solution set of
$$x^2\equiv 1 \pmod n.$$
A brute-force scan over all \(1\le x<n\) would be far too slow, so the implementation builds these roots from the prime-power factorization of \(n\).
Write the factorization of \(n\) as
$$n=\prod_{j=1}^{t} q_j,$$
where each \(q_j\) is a prime power and the factors are pairwise coprime.
For an odd prime power \(p^a\), the congruence \(x^2\equiv 1 \pmod{p^a}\) has exactly two solutions:
$$x\equiv \pm 1 \pmod{p^a}.$$
The reason is that \(p^a\mid (x-1)(x+1)\), while \(\gcd(x-1,x+1)\mid 2\). Since \(p\) is odd, the two factors cannot both absorb a power of \(p\), so one must contain all of \(p^a\).
The power of \(2\) is the only exceptional case:
$$x^2\equiv 1 \pmod{2^e}\quad\text{has}\quad \begin{cases} 1\text{ solution}, & e=1,\\ 2\text{ solutions}, & e=2,\\ 4\text{ solutions}, & e\ge 3. \end{cases}$$
For \(e\ge 3\), the four residues are
$$x\equiv \pm 1,\qquad x\equiv 2^{e-1}\pm 1 \pmod{2^e}.$$
Indeed,
$$\left(2^{e-1}\pm 1\right)^2=2^{2e-2}\pm 2^e+1\equiv 1 \pmod{2^e}.$$
This special \(2\)-power behaviour is exactly why the implementation has two extra branches only when the factor \(2^e\) satisfies \(e\ge 3\).
Because the prime-power factors \(q_1,\dots,q_t\) are pairwise coprime, the Chinese Remainder Theorem says that a solution modulo \(n\) is the same thing as choosing one local root modulo each \(q_j\). Therefore the full solution set is the Cartesian product of the local solution sets.
If \(r\) is the number of distinct odd prime divisors of \(n\), then the total number of roots is
$$R(n)=2^r\cdot \begin{cases} 1, & 2\nmid n\text{ or }2\parallel n,\\ 2, & 4\mid n\text{ but }8\nmid n,\\ 4, & 8\mid n. \end{cases}$$
So the implementation never searches through all residues modulo \(n\); it only enumerates these CRT combinations.
For each prime-power factor \(q_j\), define
$$M_j=\frac{n}{q_j}.$$
Since \(\gcd(M_j,q_j)=1\), there exists an inverse \(u_j\) such that
$$M_j u_j\equiv 1 \pmod{q_j}.$$
Now set
$$E_j\equiv M_j u_j \pmod n.$$
Then \(E_j\) acts like a CRT projector:
$$E_j\equiv 1 \pmod{q_j},\qquad E_j\equiv 0 \pmod{q_i}\quad(i\ne j).$$
The implementations begin with the global residue
$$x_0=n-1,$$
which corresponds to choosing the local root \(-1\) on every prime-power factor.
If a factor \(q_j\) only has the two roots \(\{-1,+1\}\), then switching the local choice from \(-1\) to \(+1\) changes that component by \(2\). The lifted global correction is therefore
$$x\longmapsto x+2E_j \pmod n.$$
Nothing changes on the other factors, because \(E_j\equiv 0\) there.
When \(q_j=2^e\) with \(e\ge 3\), the extra local roots are \(2^{e-1}-1\) and \(2^{e-1}+1\). Relative to the starting value \(-1\), their offsets are
$$2^{e-1}\qquad\text{and}\qquad 2^{e-1}+2.$$
So the two additional lifted corrections are
$$x\longmapsto x+2^{e-1}E_j \pmod n,\qquad x\longmapsto x+\left(2^{e-1}+2\right)E_j \pmod n.$$
Processing the prime powers one by one and applying these corrections to every current candidate generates every root of \(x^2\equiv 1 \pmod n\) exactly once.
This checkpoint illustrates the method cleanly. Since
$$100=4\cdot 25,$$
there are two prime-power factors.
For \(q_1=4\), we have \(M_1=25\), and \(25\equiv 1 \pmod 4\), so
$$E_1\equiv 25 \pmod{100}.$$
For \(q_2=25\), we have \(M_2=4\). The inverse of \(4\) modulo \(25\) is \(19\), hence
$$E_2\equiv 4\cdot 19=76 \pmod{100}.$$
Start from the trivial root
$$x_0=99.$$
Flipping the choice on the factor \(4\) gives
$$99+2E_1\equiv 99+50\equiv 49 \pmod{100}.$$
Flipping the choice on the factor \(25\) gives
$$99+2E_2\equiv 99+152\equiv 51 \pmod{100}.$$
Flipping both factors gives
$$99+2E_1+2E_2\equiv 1 \pmod{100}.$$
Thus the four roots are
$$1,\ 49,\ 51,\ 99.$$
The largest one below \(99\) is \(51\), so
$$I(100)=51.$$
The C++, Python, and Java implementations first build a smallest-prime-factor sieve up to the overall limit. This allows each \(n\) to be factored into prime powers quickly.
For a fixed \(n\), the implementation starts from the residue \(n-1\), factors \(n\), and for each prime-power factor constructs the CRT projector described above. The modular inverse needed for that projector is computed with the extended Euclidean algorithm.
Each odd prime power, and also the factor \(4\), contributes one extra branch obtained by replacing the local residue \(-1\) with \(+1\). A factor \(2^e\) with \(e\ge 3\) contributes two more branches coming from the additional roots \(2^{e-1}\pm 1\). The candidate list is expanded incrementally, and the largest candidate strictly below \(n-1\) is stored as \(I(n)\).
Because the algorithm works with the full solution set of \(x^2\equiv 1 \pmod n\), it is exact; there is no heuristic search and no dependence on trial residues.
Let \(N=2\cdot 10^7\). Building the smallest-prime-factor sieve costs \(O(N\log\log N)\) time and \(O(N)\) memory.
For one value of \(n\), extracting the prime powers is close to logarithmic in practice, and the enumeration cost is proportional to the number of roots of \(x^2\equiv 1 \pmod n\), not to \(n\) itself. Since that root count is determined only by the distinct prime-power factors of \(n\), it stays very small compared with \(n\). This is what makes the full range up to \(2\cdot 10^7\) feasible.
Für jedes \(n \ge 3\) ist definiert
$$I(n)=\max\{m\in\mathbb{Z}_{>0}: m<n-1,\ m^{-1}\equiv m \pmod n\}.$$
Die Bedingung \(m^{-1}\equiv m \pmod n\) bedeutet, dass \(m\) sein eigenes multiplikatives Inverses modulo \(n\) ist. Durch Multiplikation mit \(m\) erhält man sofort
$$m^2\equiv 1 \pmod n.$$
Gesucht ist also für jedes \(n\) die größte Lösung von \(x^2\equiv 1 \pmod n\), die strikt kleiner als \(n-1\) ist, und anschließend die Summe dieser Werte für \(3 \le n \le 2\cdot 10^7\). Die Kontrollwerte \(I(7)=1\) und \(I(100)=51\) zeigen, dass die triviale Lösung \(n-1\equiv -1\pmod n\) zwar immer existiert, aber ausdrücklich ausgeschlossen wird.
Aus \(x^2\equiv 1 \pmod n\) folgt \(n\mid (x-1)(x+1)\). Jede Lösung ist automatisch teilerfremd zu \(n\), denn ein gemeinsamer Teiler von \(x\) und \(n\) würde auch \(1\) teilen. Damit sind die ursprüngliche Inversenbedingung und die quadratische Kongruenz vollständig äquivalent.
Das Problem reduziert sich also auf die Struktur der Lösungsmenge von
$$x^2\equiv 1 \pmod n.$$
Ein naiver Test aller \(1\le x<n\) wäre viel zu langsam. Daher konstruiert die Implementierung die Lösungen aus der Primzahlpotenzzerlegung von \(n\).
Schreibe
$$n=\prod_{j=1}^{t} q_j,$$
wobei jedes \(q_j\) eine Primzahlpotenz ist und die Faktoren paarweise teilerfremd sind.
Für eine ungerade Primzahlpotenz \(p^a\) besitzt \(x^2\equiv 1 \pmod{p^a}\) genau die zwei Lösungen
$$x\equiv \pm 1 \pmod{p^a}.$$
Denn aus \(p^a\mid (x-1)(x+1)\) und \(\gcd(x-1,x+1)\mid 2\) folgt bei ungeradem \(p\), dass nicht beide Faktoren gleichzeitig eine \(p\)-Potenz tragen können. Also muss \(p^a\) vollständig in \(x-1\) oder vollständig in \(x+1\) liegen.
Die Zweierpotenz ist die einzige Ausnahme. Bezeichnet \(N_e\) die Anzahl der Lösungen von \(x^2\equiv 1 \pmod{2^e}\), so gilt
$$N_e= \begin{cases} 1, & e=1,\\ 2, & e=2,\\ 4, & e\ge 3. \end{cases}$$
Für \(e\ge 3\) lauten die vier Reste
$$x\equiv \pm 1,\qquad x\equiv 2^{e-1}\pm 1 \pmod{2^e}.$$
Tatsächlich gilt
$$\left(2^{e-1}\pm 1\right)^2=2^{2e-2}\pm 2^e+1\equiv 1 \pmod{2^e}.$$
Genau dieses Sonderverhalten erklärt, warum die Implementierung nur im Fall \(2^e\) mit \(e\ge 3\) zwei zusätzliche Verzweigungen erzeugt.
Da die Primzahlpotenzen \(q_1,\dots,q_t\) paarweise teilerfremd sind, besagt der Chinesische Restsatz: Eine Lösung modulo \(n\) entspricht genau einer Wahl lokaler Lösungen modulo jedes \(q_j\). Die globale Lösungsmenge ist also das kartesische Produkt der lokalen Lösungsmegen.
Ist \(r\) die Anzahl der verschiedenen ungeraden Primteiler von \(n\), so ist die Anzahl der Wurzeln
$$R(n)= \begin{cases} 2^r, & \nu_2(n)\le 1,\\ 2^{r+1}, & \nu_2(n)=2,\\ 2^{r+2}, & \nu_2(n)\ge 3. \end{cases}$$
Hier bezeichnet \(\nu_2(n)\) den Exponenten von \(2\) in \(n\).
Die Implementierung durchsucht daher niemals alle Reste modulo \(n\), sondern nur diese CRT-Kombinationen.
Für jeden Primzahlpotenzfaktor \(q_j\) setze
$$M_j=\frac{n}{q_j}.$$
Wegen \(\gcd(M_j,q_j)=1\) existiert ein Inverses \(u_j\) mit
$$M_j u_j\equiv 1 \pmod{q_j}.$$
Definiere dann
$$E_j\equiv M_j u_j \pmod n.$$
Dann gilt
$$E_j\equiv 1 \pmod{q_j},\qquad E_j\equiv 0 \pmod{q_i}\quad(i\ne j).$$
Die C++, Python- und Java-Implementierungen starten mit
$$x_0=n-1,$$
also mit der globalen Wahl \(-1\) auf jedem Faktor.
Hat \(q_j\) nur die beiden lokalen Lösungen \(\{-1,+1\}\), dann entspricht der Wechsel von \(-1\) nach \(+1\) einer Änderung um \(2\) auf diesem Faktor. Global hebt sich das zu
$$x\longmapsto x+2E_j \pmod n$$
an, ohne die übrigen Faktoren zu verändern.
Für \(q_j=2^e\) mit \(e\ge 3\) kommen die zusätzlichen lokalen Wurzeln \(2^{e-1}-1\) und \(2^{e-1}+1\) hinzu. Gegenüber dem Startwert \(-1\) haben sie die Offsets
$$2^{e-1},\qquad 2^{e-1}+2.$$
Daraus entstehen die zwei weiteren globalen Korrekturen
$$x\longmapsto x+2^{e-1}E_j \pmod n,\qquad x\longmapsto x+\left(2^{e-1}+2\right)E_j \pmod n.$$
Indem die Implementierung die Primzahlpotenzen nacheinander verarbeitet und diese Korrekturen auf alle aktuellen Kandidaten anwendet, erzeugt sie jede Lösung von \(x^2\equiv 1 \pmod n\) genau einmal.
Dieses Kontrollbeispiel zeigt die Methode sehr klar. Es gilt
$$100=4\cdot 25.$$
Für \(q_1=4\) ist \(M_1=25\), und \(25\equiv 1 \pmod 4\). Also
$$E_1\equiv 25 \pmod{100}.$$
Für \(q_2=25\) ist \(M_2=4\). Das Inverse von \(4\) modulo \(25\) ist \(19\), also
$$E_2\equiv 4\cdot 19=76 \pmod{100}.$$
Der Startwert ist
$$x_0=99.$$
Der Wechsel auf dem Faktor \(4\) liefert
$$99+2E_1\equiv 99+50\equiv 49 \pmod{100}.$$
Der Wechsel auf dem Faktor \(25\) liefert
$$99+2E_2\equiv 99+152\equiv 51 \pmod{100}.$$
Beide Wechsel zusammen ergeben
$$99+2E_1+2E_2\equiv 1 \pmod{100}.$$
Die vier Wurzeln sind also
$$1,\ 49,\ 51,\ 99.$$
Die größte unterhalb von \(99\) ist \(51\). Daher
$$I(100)=51.$$
Die C++, Python- und Java-Implementierungen bauen zunächst ein Sieb der kleinsten Primteiler bis zur Obergrenze auf. Damit lässt sich jedes \(n\) schnell in Primzahlpotenzen zerlegen.
Für ein festes \(n\) startet die Implementierung mit dem Rest \(n-1\), bestimmt die Primzahlpotenzen von \(n\) und konstruiert für jeden dieser Faktoren den beschriebenen CRT-Projektor. Das dafür benötigte modulare Inverse wird mit dem erweiterten Euklidischen Algorithmus berechnet.
Jede ungerade Primzahlpotenz und auch der Faktor \(4\) liefern genau einen zusätzlichen Ast durch den Wechsel von \(-1\) zu \(+1\). Ein Faktor \(2^e\) mit \(e\ge 3\) liefert wegen der zusätzlichen lokalen Wurzeln noch zwei weitere Äste. Die Kandidatenmenge wächst also kontrolliert mit der Anzahl lokaler Wahlen, und der größte Kandidat unter \(n-1\) ist genau \(I(n)\).
Weil das Verfahren die gesamte Lösungsmenge von \(x^2\equiv 1 \pmod n\) konstruiert, ist das Resultat exakt und nicht heuristisch.
Sei \(N=2\cdot 10^7\). Der Aufbau des Siebs der kleinsten Primteiler benötigt \(O(N\log\log N)\) Zeit und \(O(N)\) Speicher.
Für ein einzelnes \(n\) ist die Zerlegung in Primzahlpotenzen in der Praxis nahezu logarithmisch, und die Erzeugung der Kandidaten kostet nur so viel wie die Anzahl der Wurzeln von \(x^2\equiv 1 \pmod n\), nicht etwa \(O(n)\). Da diese Anzahl allein von den verschiedenen Primzahlpotenzfaktoren abhängt, bleibt sie gegenüber \(n\) sehr klein. Genau dadurch wird der gesamte Bereich bis \(2\cdot 10^7\) berechenbar.
Her \(n \ge 3\) için
$$I(n)=\max\{m\in\mathbb{Z}_{>0}: m<n-1,\ m^{-1}\equiv m \pmod n\}$$
tanımlanır. \(m^{-1}\equiv m \pmod n\) koşulu, \(m\)'nin mod \(n\)'de kendi çarpımsal tersi olduğu anlamına gelir. Her iki tarafı \(m\) ile çarpınca
$$m^2\equiv 1 \pmod n$$
elde edilir. Dolayısıyla her \(n\) için \(x^2\equiv 1 \pmod n\) denkleminin \(n-1\)'den küçük en büyük çözümünü bulup bunları \(3 \le n \le 2\cdot 10^7\) aralığında toplamamız gerekir. \(I(7)=1\) ve \(I(100)=51\) kontrol değerleri, \(n-1\equiv -1\pmod n\) çözümünün her zaman var olmasına rağmen özellikle dışarıda bırakıldığını gösterir.
\(x^2\equiv 1 \pmod n\) ise \(n\mid (x-1)(x+1)\) olur. Her çözüm otomatik olarak \(n\) ile aralarında asaldır; çünkü \(x\) ile \(n\)'nin ortak böleni olsaydı \(1\)'i de bölmesi gerekirdi. Böylece özgün ters koşulu ile kuadratik kongruans tam olarak aynı probleme indirgenir.
Artık anlamamız gereken şey,
$$x^2\equiv 1 \pmod n$$
denkliğinin tüm çözümleridir. \(1\le x<n\) aralığını baştan sona taramak çok pahalı olur, bu yüzden uygulama \(n\)'nin asal kuvvetlerine ayrışımını kullanır.
\(n\)'yi
$$n=\prod_{j=1}^{t} q_j$$
şeklinde yazalım; burada her \(q_j\) bir asal kuvvettir ve bu çarpanlar ikişerli aralarında asaldır.
Tek asal kuvvet \(p^a\) için \(x^2\equiv 1 \pmod{p^a}\) denkleminin yalnızca iki çözümü vardır:
$$x\equiv \pm 1 \pmod{p^a}.$$
Bunun nedeni, \(p^a\mid (x-1)(x+1)\) iken \(\gcd(x-1,x+1)\mid 2\) olmasıdır. \(p\) tek olduğundan iki çarpan aynı anda \(p\)'nin kuvvetini paylaşamaz; bütün \(p^a\) ya \(x-1\)'i ya da \(x+1\)'i bölmek zorundadır.
\(2\)'nin kuvvetleri tek istisnadır. \(N_e\), \(x^2\equiv 1 \pmod{2^e}\) denkleminin çözüm sayısı olsun. O zaman
$$N_e= \begin{cases} 1, & e=1,\\ 2, & e=2,\\ 4, & e\ge 3. \end{cases}$$
\(e\ge 3\) olduğunda dört yerel kök şunlardır:
$$x\equiv \pm 1,\qquad x\equiv 2^{e-1}\pm 1 \pmod{2^e}.$$
Gerçekten de
$$\left(2^{e-1}\pm 1\right)^2=2^{2e-2}\pm 2^e+1\equiv 1 \pmod{2^e}.$$
Uygulamanın yalnızca \(2^e\) ve \(e\ge 3\) durumunda iki ek dal üretmesinin sebebi tam olarak budur.
\(q_1,\dots,q_t\) çarpanları ikişerli aralarında asal olduğundan, Çin Kalan Teoremi şunu söyler: mod \(n\)'de bir çözüm seçmek, her \(q_j\) için bir yerel kök seçmekle aynı şeydir. Yani küresel çözüm kümesi, yerel çözüm kümelerinin Kartezyen çarpımıdır.
\(r\), \(n\)'nin farklı tek asal bölen sayısı olsun. O zaman toplam kök sayısı
$$R(n)= \begin{cases} 2^r, & \nu_2(n)\le 1,\\ 2^{r+1}, & \nu_2(n)=2,\\ 2^{r+2}, & \nu_2(n)\ge 3 \end{cases}$$
Burada \(\nu_2(n)\), \(n\) içindeki \(2\) üssünü gösterir.
olur. Dolayısıyla uygulama mod \(n\)'deki bütün kalıntıları denemez; sadece bu CRT kombinasyonlarını üretir.
Her asal kuvvet çarpanı \(q_j\) için
$$M_j=\frac{n}{q_j}$$
tanımını yapalım. \(\gcd(M_j,q_j)=1\) olduğundan
$$M_j u_j\equiv 1 \pmod{q_j}$$
sağlayan bir \(u_j\) vardır. Sonra
$$E_j\equiv M_j u_j \pmod n$$
alınır. Bu durumda
$$E_j\equiv 1 \pmod{q_j},\qquad E_j\equiv 0 \pmod{q_i}\quad(i\ne j)$$
olur; yani \(E_j\) yalnızca ilgili yerel bileşeni değiştirir.
C++, Python ve Java uygulamaları
$$x_0=n-1$$
başlangıç değeriyle, yani her yerde yerel kök \(-1\) seçimiyle başlar.
Eğer \(q_j\) çarpanında yalnızca \(\{-1,+1\}\) kökleri varsa, \(-1\)'den \(+1\)'e geçmek yerelde \(2\) fark yaratır. Küresel düzeltme bu yüzden
$$x\longmapsto x+2E_j \pmod n$$
şeklindedir.
\(q_j=2^e\) ve \(e\ge 3\) olduğunda ek kökler \(2^{e-1}-1\) ile \(2^{e-1}+1\) olur. Başlangıç değeri \(-1\)'e göre bunların farkları
$$2^{e-1},\qquad 2^{e-1}+2$$
olduğu için iki ek küresel düzeltme daha gelir:
$$x\longmapsto x+2^{e-1}E_j \pmod n,\qquad x\longmapsto x+\left(2^{e-1}+2\right)E_j \pmod n.$$
Uygulama asal kuvvetleri tek tek işler ve bu düzeltmeleri mevcut tüm adaylara uygular; böylece \(x^2\equiv 1 \pmod n\) denkleminin her çözümü tam bir kez üretilmiş olur.
Kontrol örneği yöntemi açık biçimde gösterir. Burada
$$100=4\cdot 25.$$
\(q_1=4\) için \(M_1=25\) ve \(25\equiv 1 \pmod 4\), dolayısıyla
$$E_1\equiv 25 \pmod{100}.$$
\(q_2=25\) için \(M_2=4\). \(4\)'ün mod \(25\)'te tersi \(19\) olduğundan
$$E_2\equiv 4\cdot 19=76 \pmod{100}.$$
Başlangıç kökü
$$x_0=99$$
olur.
\(4\) çarpanındaki seçimi çevirmek
$$99+2E_1\equiv 99+50\equiv 49 \pmod{100}$$
sonucunu verir.
\(25\) çarpanındaki seçimi çevirmek
$$99+2E_2\equiv 99+152\equiv 51 \pmod{100}$$
sonucunu verir.
Her ikisini birlikte çevirmek ise
$$99+2E_1+2E_2\equiv 1 \pmod{100}$$
olur. Böylece dört kök
$$1,\ 49,\ 51,\ 99$$
şeklindedir. \(99\)'dan küçük en büyük kök \(51\) olduğu için
$$I(100)=51$$
elde edilir.
C++, Python ve Java uygulamaları önce üst sınıra kadar en küçük asal bölen eleğini kurar. Böylece her \(n\) hızlıca asal kuvvet çarpanlarına ayrılabilir.
Sabit bir \(n\) için uygulama \(n-1\) kalıntısından başlar, \(n\)'yi çarpanlarına ayırır ve her asal kuvvet için yukarıdaki CRT projektörünü kurar. Gerekli modüler ters, genişletilmiş Öklid algoritması ile hesaplanır.
Her tek asal kuvvet ve ayrıca \(4\) çarpanı, yerel \(-1\) seçimini \(+1\)'e çeviren tek bir ek dal üretir. \(2^e\) ve \(e\ge 3\) durumunda ise \(2^{e-1}\pm 1\) köklerinden gelen iki ek dal daha oluşur. Aday listesi bu küçük yerel seçim kümesine göre büyür ve \(n-1\)'den küçük en büyük aday tam olarak \(I(n)\) olur.
Algoritma tüm çözüm kümesini doğrudan ürettiği için sonuç tamdır; deneme-yanılma veya rastgele arama yoktur.
\(N=2\cdot 10^7\) olsun. En küçük asal bölen eleğinin kurulması \(O(N\log\log N)\) zaman ve \(O(N)\) bellek gerektirir.
Tek bir \(n\) için asal kuvvetleri çıkarmak pratikte logaritmik ölçektedir; aday üretme maliyeti ise \(n\)'ye değil, yalnızca \(x^2\equiv 1 \pmod n\) denkleminin kök sayısına bağlıdır. Bu sayı da sadece farklı asal kuvvet çarpanlarının sayısıyla belirlenir ve \(n\)'ye kıyasla çok küçük kalır. Bu yüzden \(2\cdot 10^7\) sınırı uygulanabilir olur.
Para cada entero \(n \ge 3\), se define
$$I(n)=\max\{m\in\mathbb{Z}_{>0}: m<n-1,\ m^{-1}\equiv m \pmod n\}.$$
La condición \(m^{-1}\equiv m \pmod n\) significa que \(m\) es su propio inverso multiplicativo módulo \(n\). Al multiplicar por \(m\), obtenemos
$$m^2\equiv 1 \pmod n.$$
Así, para cada \(n\) debemos hallar la mayor solución de \(x^2\equiv 1 \pmod n\) que sea estrictamente menor que \(n-1\), y sumar esos valores para \(3 \le n \le 2\cdot 10^7\). Los puntos de control \(I(7)=1\) e \(I(100)=51\) muestran que \(n-1\equiv -1\pmod n\) siempre es una solución, pero no cuenta porque el problema pide la mayor solución no trivial.
Si \(x^2\equiv 1 \pmod n\), entonces \(n\mid (x-1)(x+1)\). Toda solución es automáticamente coprima con \(n\), porque cualquier divisor común de \(x\) y \(n\) también dividiría \(1\). Por lo tanto, la condición original sobre inversos y la congruencia cuadrática son exactamente el mismo problema.
El objetivo pasa a ser entender todas las soluciones de
$$x^2\equiv 1 \pmod n.$$
Probar todos los residuos \(1\le x<n\) sería demasiado lento, así que la implementación construye las soluciones a partir de la factorización de \(n\) en potencias primas.
Escribamos
$$n=\prod_{j=1}^{t} q_j,$$
donde cada \(q_j\) es una potencia prima y los factores son coprimos dos a dos.
Para una potencia prima impar \(p^a\), la congruencia \(x^2\equiv 1 \pmod{p^a}\) tiene exactamente dos soluciones:
$$x\equiv \pm 1 \pmod{p^a}.$$
La razón es que \(p^a\mid (x-1)(x+1)\), mientras que \(\gcd(x-1,x+1)\mid 2\). Como \(p\) es impar, los dos factores no pueden compartir una potencia de \(p\), así que toda la potencia \(p^a\) debe caer en \(x-1\) o en \(x+1\).
La única excepción es la potencia de \(2\). Si \(N_e\) denota el número de soluciones de \(x^2\equiv 1 \pmod{2^e}\), entonces
$$N_e= \begin{cases} 1, & e=1,\\ 2, & e=2,\\ 4, & e\ge 3. \end{cases}$$
Cuando \(e\ge 3\), las cuatro raíces locales son
$$x\equiv \pm 1,\qquad x\equiv 2^{e-1}\pm 1 \pmod{2^e}.$$
En efecto,
$$\left(2^{e-1}\pm 1\right)^2=2^{2e-2}\pm 2^e+1\equiv 1 \pmod{2^e}.$$
Este comportamiento especial de \(2^e\) explica por qué la implementación añade dos ramas extra solo cuando aparece un factor \(2^e\) con \(e\ge 3\).
Como \(q_1,\dots,q_t\) son coprimos dos a dos, el teorema chino del resto afirma que una solución módulo \(n\) equivale a elegir una raíz local para cada \(q_j\). En consecuencia, el conjunto global de soluciones es el producto cartesiano de los conjuntos locales.
Si \(r\) es el número de divisores primos impares distintos de \(n\), entonces el número total de raíces es
$$R(n)= \begin{cases} 2^r, & \nu_2(n)\le 1,\\ 2^{r+1}, & \nu_2(n)=2,\\ 2^{r+2}, & \nu_2(n)\ge 3. \end{cases}$$
Aquí \(\nu_2(n)\) denota el exponente de \(2\) en la factorización de \(n\).
Por eso la implementación no examina todos los residuos módulo \(n\); solo enumera estas combinaciones CRT.
Para cada factor de potencia prima \(q_j\), definimos
$$M_j=\frac{n}{q_j}.$$
Como \(\gcd(M_j,q_j)=1\), existe un inverso \(u_j\) tal que
$$M_j u_j\equiv 1 \pmod{q_j}.$$
Ahora ponemos
$$E_j\equiv M_j u_j \pmod n.$$
Entonces \(E_j\) actúa como un proyector CRT:
$$E_j\equiv 1 \pmod{q_j},\qquad E_j\equiv 0 \pmod{q_i}\quad(i\ne j).$$
Las implementaciones en C++, Python y Java comienzan con
$$x_0=n-1,$$
que corresponde a elegir la raíz local \(-1\) en todos los factores.
Si un factor \(q_j\) solo tiene las dos raíces \(\{-1,+1\}\), cambiar de \(-1\) a \(+1\) altera ese componente local en \(2\). La corrección global levantada es entonces
$$x\longmapsto x+2E_j \pmod n.$$
Como \(E_j\equiv 0\) en los demás factores, nada cambia allí.
Cuando \(q_j=2^e\) con \(e\ge 3\), aparecen además las raíces \(2^{e-1}-1\) y \(2^{e-1}+1\). Comparadas con el valor inicial \(-1\), producen los desplazamientos
$$2^{e-1},\qquad 2^{e-1}+2.$$
Por tanto, las dos correcciones globales adicionales son
$$x\longmapsto x+2^{e-1}E_j \pmod n,\qquad x\longmapsto x+\left(2^{e-1}+2\right)E_j \pmod n.$$
Procesando las potencias primas una por una y aplicando estas correcciones a cada candidato actual, la implementación genera exactamente una vez todas las raíces de \(x^2\equiv 1 \pmod n\).
Este punto de control ilustra el método de forma directa. Como
$$100=4\cdot 25,$$
tenemos dos factores de potencia prima.
Para \(q_1=4\), se tiene \(M_1=25\) y \(25\equiv 1 \pmod 4\), así que
$$E_1\equiv 25 \pmod{100}.$$
Para \(q_2=25\), se tiene \(M_2=4\). El inverso de \(4\) módulo \(25\) es \(19\), por lo que
$$E_2\equiv 4\cdot 19=76 \pmod{100}.$$
La raíz inicial es
$$x_0=99.$$
Cambiar la elección en el factor \(4\) da
$$99+2E_1\equiv 99+50\equiv 49 \pmod{100}.$$
Cambiar la elección en el factor \(25\) da
$$99+2E_2\equiv 99+152\equiv 51 \pmod{100}.$$
Cambiar ambos factores produce
$$99+2E_1+2E_2\equiv 1 \pmod{100}.$$
Por lo tanto, las cuatro raíces son
$$1,\ 49,\ 51,\ 99.$$
La mayor estrictamente menor que \(99\) es \(51\), así que
$$I(100)=51.$$
Las implementaciones en C++, Python y Java construyen primero una criba del menor factor primo hasta el límite total. Eso permite factorizar rápidamente cada \(n\) en potencias primas.
Para un \(n\) fijo, la implementación parte del residuo \(n-1\), descompone \(n\) y construye para cada potencia prima el proyector CRT descrito arriba. El inverso modular necesario para ese proyector se obtiene con el algoritmo extendido de Euclides.
Cada potencia prima impar, y también el factor \(4\), aporta una rama adicional al cambiar la raíz local \(-1\) por \(+1\). Un factor \(2^e\) con \(e\ge 3\) aporta además dos ramas extra procedentes de las raíces \(2^{e-1}\pm 1\). La lista de candidatos se amplía de manera incremental, y el mayor candidato estrictamente menor que \(n-1\) es precisamente \(I(n)\).
Como el algoritmo construye el conjunto completo de soluciones de \(x^2\equiv 1 \pmod n\), el resultado es exacto y no depende de una búsqueda heurística.
Sea \(N=2\cdot 10^7\). Construir la criba del menor factor primo cuesta \(O(N\log\log N)\) tiempo y \(O(N)\) memoria.
Para un valor concreto de \(n\), extraer las potencias primas es casi logarítmico en la práctica, y el coste de enumeración es proporcional al número de raíces de \(x^2\equiv 1 \pmod n\), no a \(n\). Como ese número depende solo de los distintos factores de potencia prima de \(n\), permanece muy pequeño en comparación con \(n\). Por eso el rango completo hasta \(2\cdot 10^7\) resulta viable.
对每个 \(n \ge 3\)，定义
$$I(n)=\max\{m\in\mathbb{Z}_{>0}: m<n-1,\ m^{-1}\equiv m \pmod n\}.$$
条件 \(m^{-1}\equiv m \pmod n\) 表示 \(m\) 在模 \(n\) 意义下是自己的乘法逆元。两边同乘 \(m\)，立刻得到
$$m^2\equiv 1 \pmod n.$$
因此，题目的本质是：对每个 \(n\)，找出方程 \(x^2\equiv 1 \pmod n\) 中严格小于 \(n-1\) 的最大解，然后把这些值在 \(3 \le n \le 2\cdot 10^7\) 上求和。已知校验点 \(I(7)=1\)、\(I(100)=51\)，说明 \(n-1\equiv -1\pmod n\) 这个恒成立的平凡解必须排除。
若 \(x^2\equiv 1 \pmod n\)，则有 \(n\mid (x-1)(x+1)\)。任何这样的解都自动与 \(n\) 互素，因为如果某个数同时整除 \(x\) 和 \(n\)，它也必须整除 \(1\)。所以原题中的“逆元等于自身”和二次同余完全等价。
问题因此化为研究
$$x^2\equiv 1 \pmod n$$
的全部解。如果对每个 \(n\) 都从 \(1\) 扫到 \(n-1\)，代价会过高，所以实现采用 \(n\) 的素数幂分解来直接构造解集。
把 \(n\) 写成
$$n=\prod_{j=1}^{t} q_j,$$
其中每个 \(q_j\) 都是一个素数幂，而且这些因子两两互素。
对于奇素数幂 \(p^a\)，同余 \(x^2\equiv 1 \pmod{p^a}\) 恰好只有两个解：
$$x\equiv \pm 1 \pmod{p^a}.$$
原因是 \(p^a\mid (x-1)(x+1)\)，而 \(\gcd(x-1,x+1)\mid 2\)。由于 \(p\) 是奇素数，\(x-1\) 与 \(x+1\) 不可能同时吸收 \(p\) 的幂，因此全部的 \(p^a\) 只能落在其中一个因子上。
\(2\) 的幂是唯一的例外。记 \(N_e\) 为 \(x^2\equiv 1 \pmod{2^e}\) 的解数，则
$$N_e= \begin{cases} 1, & e=1,\\ 2, & e=2,\\ 4, & e\ge 3. \end{cases}$$
当 \(e\ge 3\) 时，四个局部解为
$$x\equiv \pm 1,\qquad x\equiv 2^{e-1}\pm 1 \pmod{2^e}.$$
因为
$$\left(2^{e-1}\pm 1\right)^2=2^{2e-2}\pm 2^e+1\equiv 1 \pmod{2^e}.$$
这正是实现中只有在出现 \(2^e\) 且 \(e\ge 3\) 时才增加两条额外分支的原因。
由于 \(q_1,\dots,q_t\) 两两互素，中国剩余定理告诉我们：模 \(n\) 的一个解，等价于在每个 \(q_j\) 上各选一个局部解。因此，全局解集就是各局部解集的笛卡尔积。
若 \(r\) 表示 \(n\) 的不同奇素因子的个数，则根的总数为
$$R(n)= \begin{cases} 2^r, & \nu_2(n)\le 1,\\ 2^{r+1}, & \nu_2(n)=2,\\ 2^{r+2}, & \nu_2(n)\ge 3. \end{cases}$$
这里 \(\nu_2(n)\) 表示 \(n\) 中 \(2\) 的指数。
所以实现从不枚举模 \(n\) 的全部剩余类，而只枚举这些 CRT 组合。
对每个素数幂因子 \(q_j\)，定义
$$M_j=\frac{n}{q_j}.$$
因为 \(\gcd(M_j,q_j)=1\)，存在逆元 \(u_j\) 满足
$$M_j u_j\equiv 1 \pmod{q_j}.$$
再令
$$E_j\equiv M_j u_j \pmod n.$$
那么 \(E_j\) 就是一个 CRT 投影量：
$$E_j\equiv 1 \pmod{q_j},\qquad E_j\equiv 0 \pmod{q_i}\quad(i\ne j).$$
C++、Python 和 Java 实现都从
$$x_0=n-1$$
开始，这对应于在每个局部模数上都选择根 \(-1\)。
如果 \(q_j\) 只有 \(\{-1,+1\}\) 两个局部解，那么把该因子上的选择从 \(-1\) 改成 \(+1\)，局部变化量就是 \(2\)。提升到模 \(n\) 后，对应的全局修正为
$$x\longmapsto x+2E_j \pmod n.$$
由于 \(E_j\) 在其他因子上同余于 \(0\)，其余局部条件不会被改变。
当 \(q_j=2^e\) 且 \(e\ge 3\) 时，还有额外局部根 \(2^{e-1}-1\) 与 \(2^{e-1}+1\)。相对于起点 \(-1\)，它们的偏移量分别是
$$2^{e-1},\qquad 2^{e-1}+2.$$
因此还会多出两种全局修正：
$$x\longmapsto x+2^{e-1}E_j \pmod n,\qquad x\longmapsto x+\left(2^{e-1}+2\right)E_j \pmod n.$$
实现按素数幂逐个处理，并把这些修正施加到当前所有候选值上，于是 \(x^2\equiv 1 \pmod n\) 的每个解都会被恰好生成一次。
这个校验点很适合展示算法。因为
$$100=4\cdot 25,$$
所以有两个素数幂因子。
对 \(q_1=4\)，有 \(M_1=25\)，而 \(25\equiv 1 \pmod 4\)，因此
$$E_1\equiv 25 \pmod{100}.$$
对 \(q_2=25\)，有 \(M_2=4\)。\(4\) 在模 \(25\) 下的逆元是 \(19\)，于是
$$E_2\equiv 4\cdot 19=76 \pmod{100}.$$
起始根为
$$x_0=99.$$
若翻转模 \(4\) 上的选择，则得到
$$99+2E_1\equiv 99+50\equiv 49 \pmod{100}.$$
若翻转模 \(25\) 上的选择，则得到
$$99+2E_2\equiv 99+152\equiv 51 \pmod{100}.$$
若两个因子都翻转，则得到
$$99+2E_1+2E_2\equiv 1 \pmod{100}.$$
因此四个根是
$$1,\ 49,\ 51,\ 99.$$
严格小于 \(99\) 的最大值是 \(51\)，所以
$$I(100)=51.$$
C++、Python 和 Java 实现首先建立到总上界为止的最小素因子筛，这样每个 \(n\) 都能很快分解成素数幂。
对固定的 \(n\)，实现从剩余类 \(n-1\) 出发，分解 \(n\)，并为每个素数幂构造上面的 CRT 投影量。构造过程中所需的模逆由扩展欧几里得算法计算。
每个奇素数幂以及因子 \(4\) 都只会产生一条额外分支，对应把局部根 \(-1\) 改成 \(+1\)。若存在 \(2^e\) 且 \(e\ge 3\)，则还会因为额外局部根 \(2^{e-1}\pm 1\) 再产生两条分支。候选集合以这种受控方式扩张，而严格小于 \(n-1\) 的最大候选值就是 \(I(n)\)。
由于算法构造的是完整解集而不是试探性搜索，所以结果是严格正确的。
令 \(N=2\cdot 10^7\)。构造最小素因子筛需要 \(O(N\log\log N)\) 时间和 \(O(N)\) 空间。
对单个 \(n\) 而言，分解成素数幂在实践中接近对数级，而枚举成本与 \(x^2\equiv 1 \pmod n\) 的根数成正比，而不是与 \(n\) 成正比。根数只由不同素数幂因子的数量决定，因此相对于 \(n\) 本身非常小。这正是算法能够处理到 \(2\cdot 10^7\) 的关键。
Для каждого \(n \ge 3\) определяется
$$I(n)=\max\{m\in\mathbb{Z}_{>0}: m<n-1,\ m^{-1}\equiv m \pmod n\}.$$
Условие \(m^{-1}\equiv m \pmod n\) означает, что число \(m\) является собственным мультипликативным обратным по модулю \(n\). Умножая на \(m\), получаем
$$m^2\equiv 1 \pmod n.$$
Значит, для каждого \(n\) нужно найти наибольшее решение сравнения \(x^2\equiv 1 \pmod n\), строго меньшее \(n-1\), а затем просуммировать эти значения для \(3 \le n \le 2\cdot 10^7\). Контрольные значения \(I(7)=1\) и \(I(100)=51\) показывают, что тривиальный корень \(n-1\equiv -1\pmod n\) всегда существует, но по условию не учитывается.
Если \(x^2\equiv 1 \pmod n\), то \(n\mid (x-1)(x+1)\). Любое такое решение автоматически взаимно просто с \(n\), потому что общий делитель \(x\) и \(n\) делил бы и \(1\). Поэтому исходное условие про обратный элемент и квадратное сравнение полностью эквивалентны.
Тем самым задача сводится к описанию всех решений
$$x^2\equiv 1 \pmod n.$$
Перебирать все \(1\le x<n\) слишком дорого, поэтому реализация строит решения из разложения \(n\) на простые степени.
Запишем
$$n=\prod_{j=1}^{t} q_j,$$
где каждое \(q_j\) является степенью простого числа, а все множители попарно взаимно просты.
Для нечётной простой степени \(p^a\) сравнение \(x^2\equiv 1 \pmod{p^a}\) имеет ровно два решения:
$$x\equiv \pm 1 \pmod{p^a}.$$
Действительно, \(p^a\mid (x-1)(x+1)\), а \(\gcd(x-1,x+1)\mid 2\). Поскольку \(p\) нечётно, два сомножителя не могут одновременно содержать степень \(p\), значит вся степень \(p^a\) должна делить либо \(x-1\), либо \(x+1\).
Степени двойки составляют единственное исключение. Если \(N_e\) обозначает число решений сравнения \(x^2\equiv 1 \pmod{2^e}\), то
$$N_e= \begin{cases} 1, & e=1,\\ 2, & e=2,\\ 4, & e\ge 3. \end{cases}$$
При \(e\ge 3\) четыре локальных корня имеют вид
$$x\equiv \pm 1,\qquad x\equiv 2^{e-1}\pm 1 \pmod{2^e}.$$
Потому что
$$\left(2^{e-1}\pm 1\right)^2=2^{2e-2}\pm 2^e+1\equiv 1 \pmod{2^e}.$$
Именно это особое поведение объясняет, почему реализация добавляет две дополнительные ветви только для множителя \(2^e\) при \(e\ge 3\).
Так как \(q_1,\dots,q_t\) попарно взаимно просты, китайская теорема об остатках говорит, что решение по модулю \(n\) однозначно задаётся выбором локального корня по каждому модулю \(q_j\). Значит, глобальное множество решений есть декартово произведение локальных множеств.
Если \(r\) обозначает число различных нечётных простых делителей числа \(n\), то число корней равно
$$R(n)= \begin{cases} 2^r, & \nu_2(n)\le 1,\\ 2^{r+1}, & \nu_2(n)=2,\\ 2^{r+2}, & \nu_2(n)\ge 3. \end{cases}$$
Здесь \(\nu_2(n)\) означает показатель степени двойки в разложении числа \(n\).
Поэтому реализация не рассматривает все остатки по модулю \(n\), а перечисляет только эти CRT-комбинации.
Для каждого множителя \(q_j\) положим
$$M_j=\frac{n}{q_j}.$$
Поскольку \(\gcd(M_j,q_j)=1\), существует обратный элемент \(u_j\), такой что
$$M_j u_j\equiv 1 \pmod{q_j}.$$
Затем определим
$$E_j\equiv M_j u_j \pmod n.$$
Тогда \(E_j\) является CRT-проектором:
$$E_j\equiv 1 \pmod{q_j},\qquad E_j\equiv 0 \pmod{q_i}\quad(i\ne j).$$
Реализации на C++, Python и Java начинают с остатка
$$x_0=n-1,$$
что соответствует выбору локального корня \(-1\) на каждом простом степенном множителе.
Если для \(q_j\) локальные корни равны только \(\{-1,+1\}\), то переход от \(-1\) к \(+1\) меняет этот компонент на \(2\). Подъём на модуль \(n\) даёт поправку
$$x\longmapsto x+2E_j \pmod n.$$
На остальных множителях ничего не меняется, потому что там \(E_j\equiv 0\).
Когда \(q_j=2^e\) и \(e\ge 3\), возникают дополнительные локальные корни \(2^{e-1}-1\) и \(2^{e-1}+1\). Относительно стартового значения \(-1\) они дают смещения
$$2^{e-1},\qquad 2^{e-1}+2.$$
Соответственно появляются ещё две глобальные поправки:
$$x\longmapsto x+2^{e-1}E_j \pmod n,\qquad x\longmapsto x+\left(2^{e-1}+2\right)E_j \pmod n.$$
Последовательно обрабатывая простые степени и применяя такие поправки ко всем текущим кандидатам, реализация строит каждое решение сравнения \(x^2\equiv 1 \pmod n\) ровно один раз.
Этот контрольный пример хорошо демонстрирует метод. Имеем
$$100=4\cdot 25.$$
Для \(q_1=4\) получаем \(M_1=25\), причём \(25\equiv 1 \pmod 4\), поэтому
$$E_1\equiv 25 \pmod{100}.$$
Для \(q_2=25\) имеем \(M_2=4\). Обратный к \(4\) по модулю \(25\) равен \(19\), значит
$$E_2\equiv 4\cdot 19=76 \pmod{100}.$$
Стартовый корень равен
$$x_0=99.$$
Переключение на факторе \(4\) даёт
$$99+2E_1\equiv 99+50\equiv 49 \pmod{100}.$$
Переключение на факторе \(25\) даёт
$$99+2E_2\equiv 99+152\equiv 51 \pmod{100}.$$
Переключение обоих факторов приводит к
$$99+2E_1+2E_2\equiv 1 \pmod{100}.$$
Следовательно, четыре корня таковы:
$$1,\ 49,\ 51,\ 99.$$
Наибольший корень, меньший \(99\), равен \(51\). Значит,
$$I(100)=51.$$
Реализации на C++, Python и Java сначала строят решето минимального простого делителя до общей границы. Это позволяет быстро разлагать каждое \(n\) на простые степени.
Для фиксированного \(n\) реализация начинает с остатка \(n-1\), факторизует число и для каждой простой степени строит описанный выше CRT-проектор. Необходимый модульный обратный элемент вычисляется расширенным алгоритмом Евклида.
Каждая нечётная простая степень, а также множитель \(4\), добавляет одну новую ветвь, соответствующую переходу от локального корня \(-1\) к \(+1\). Множитель \(2^e\) при \(e\ge 3\) добавляет ещё две ветви из-за корней \(2^{e-1}\pm 1\). Множество кандидатов растёт только вместе с числом локальных вариантов, и наибольший кандидат, строго меньший \(n-1\), и есть \(I(n)\).
Поскольку алгоритм перечисляет полное множество решений \(x^2\equiv 1 \pmod n\), он даёт точный результат без эвристик.
Пусть \(N=2\cdot 10^7\). Построение решета минимального простого делителя требует \(O(N\log\log N)\) времени и \(O(N)\) памяти.
Для одного значения \(n\) выделение простых степеней на практике близко к логарифмическому, а стоимость перечисления пропорциональна числу корней сравнения \(x^2\equiv 1 \pmod n\), а не самому \(n\). Это число зависит только от различных простых степенных множителей числа \(n\) и остаётся очень малым по сравнению с \(n\). Именно поэтому диапазон до \(2\cdot 10^7\) оказывается вычислимым.
لكل عدد صحيح \(n \ge 3\) نعرّف
$$I(n)=\max\{m\in\mathbb{Z}_{>0}: m<n-1,\ m^{-1}\equiv m \pmod n\}.$$
الشرط \(m^{-1}\equiv m \pmod n\) يعني أن \(m\) هو معكوسه الضربي نفسه بترديد \(n\). وإذا ضربنا الطرفين في \(m\) نحصل على
$$m^2\equiv 1 \pmod n.$$
إذًا المطلوب هو إيجاد أكبر حل للمعادلة \(x^2\equiv 1 \pmod n\) يكون أصغر تمامًا من \(n-1\)، ثم جمع هذه القيم لكل \(3 \le n \le 2\cdot 10^7\). والقيمتان \(I(7)=1\) و\(I(100)=51\) توضحان أن الحل التافه \(n-1\equiv -1\pmod n\) موجود دائمًا لكنه مستبعد من التعريف.
إذا كان \(x^2\equiv 1 \pmod n\)، فهذا يعني أن \(n\mid (x-1)(x+1)\). وكل حل من هذه الحلول يكون تلقائيًا أوليًا نسبيًا مع \(n\)، لأن أي قاسم مشترك بين \(x\) و\(n\) سيلزم أن يقسم \(1\). لذلك فشرط المعكوس الضربي الأصلي يكافئ تمامًا المقارنة التربيعية.
وبالتالي تتحول المسألة إلى فهم جميع حلول
$$x^2\equiv 1 \pmod n.$$
الفحص المباشر لكل \(1\le x<n\) غير عملي، ولهذا تبني التطبيقات الحلول انطلاقًا من تحليل \(n\) إلى قوى أولية.
لنكتب
$$n=\prod_{j=1}^{t} q_j,$$
حيث إن كل \(q_j\) قوة أولية، وهذه العوامل أولية نسبيًا فيما بينها اثنين اثنين.
بالنسبة إلى قوة أولية فردية \(p^a\)، فإن المعادلة \(x^2\equiv 1 \pmod{p^a}\) لها حلان فقط:
$$x\equiv \pm 1 \pmod{p^a}.$$
والسبب هو أن \(p^a\mid (x-1)(x+1)\)، بينما \(\gcd(x-1,x+1)\mid 2\). وبما أن \(p\) فردي، فلا يمكن للعاملين أن يشتركا معًا في قوة من \(p\)، ولذلك يجب أن تقع كل القوة \(p^a\) بالكامل في \(x-1\) أو بالكامل في \(x+1\).
قوى العدد \(2\) هي الحالة الاستثنائية الوحيدة. إذا رمزنا بعدد حلول المعادلة \(x^2\equiv 1 \pmod{2^e}\) بالرمز \(N_e\)، فإن
$$N_e= \begin{cases} 1, & e=1,\\ 2, & e=2,\\ 4, & e\ge 3. \end{cases}$$
وعندما \(e\ge 3\)، تكون الجذور المحلية الأربعة هي
$$x\equiv \pm 1,\qquad x\equiv 2^{e-1}\pm 1 \pmod{2^e}.$$
وذلك لأن
$$\left(2^{e-1}\pm 1\right)^2=2^{2e-2}\pm 2^e+1\equiv 1 \pmod{2^e}.$$
وهذا السلوك الخاص هو بالضبط سبب إضافة التطبيق فرعين إضافيين فقط عندما يظهر العامل \(2^e\) مع \(e\ge 3\).
بما أن العوامل \(q_1,\dots,q_t\) أولية نسبيًا فيما بينها، فإن مبرهنة البواقي الصينية تقول إن كل حل بترديد \(n\) يقابل اختيار جذر محلي واحد بترديد كل \(q_j\). لذلك فإن مجموعة الحلول الكلية هي حاصل الضرب الديكارتي لمجموعات الحلول المحلية.
إذا كان \(r\) هو عدد القواسم الأولية الفردية المميزة للعدد \(n\)، فإن عدد الجذور يساوي
$$R(n)= \begin{cases} 2^r, & \nu_2(n)\le 1,\\ 2^{r+1}, & \nu_2(n)=2,\\ 2^{r+2}, & \nu_2(n)\ge 3. \end{cases}$$
وهنا تمثل \(\nu_2(n)\) أس \(2\) في تحليل \(n\).
ولهذا لا تفحص التطبيقات جميع البواقي بترديد \(n\)، بل تعدد فقط هذه التركيبات الناتجة من CRT.
لكل عامل قوة أولية \(q_j\) نعرّف
$$M_j=\frac{n}{q_j}.$$
وبما أن \(\gcd(M_j,q_j)=1\)، يوجد معكوس \(u_j\) يحقق
$$M_j u_j\equiv 1 \pmod{q_j}.$$
ثم نضع
$$E_j\equiv M_j u_j \pmod n.$$
وعندها يعمل \(E_j\) كمسقط CRT:
$$E_j\equiv 1 \pmod{q_j},\qquad E_j\equiv 0 \pmod{q_i}\quad(i\ne j).$$
تبدأ تطبيقات C++ وPython وJava من الباقي
$$x_0=n-1,$$
وهو ما يقابل اختيار الجذر المحلي \(-1\) في كل العوامل.
إذا كان \(q_j\) لا يملك إلا الجذرين المحليين \(\{-1,+1\}\)، فإن الانتقال من \(-1\) إلى \(+1\) يضيف فرقًا مقداره \(2\) على هذا العامل. والرفع إلى الترديد \(n\) يعطي التصحيح
$$x\longmapsto x+2E_j \pmod n.$$
أما العوامل الأخرى فلا تتغير لأن \(E_j\equiv 0\) عليها.
وعندما يكون \(q_j=2^e\) مع \(e\ge 3\)، تظهر الجذور المحلية الإضافية \(2^{e-1}-1\) و\(2^{e-1}+1\). وبالنسبة إلى القيمة الابتدائية \(-1\)، فإن الإزاحتين المقابلتين هما
$$2^{e-1},\qquad 2^{e-1}+2.$$
ومن ثم نحصل على تصحيحين عالميين إضافيين:
$$x\longmapsto x+2^{e-1}E_j \pmod n,\qquad x\longmapsto x+\left(2^{e-1}+2\right)E_j \pmod n.$$
ومن خلال معالجة القوى الأولية واحدًا بعد الآخر وتطبيق هذه التصحيحات على جميع المرشحين الحاليين، تنتج التطبيقات كل حل للمعادلة \(x^2\equiv 1 \pmod n\) مرة واحدة بالضبط.
يعرض هذا المثال القياسي الطريقة بوضوح. لدينا
$$100=4\cdot 25.$$
بالنسبة إلى \(q_1=4\)، نجد أن \(M_1=25\) و\(25\equiv 1 \pmod 4\)، وبالتالي
$$E_1\equiv 25 \pmod{100}.$$
وبالنسبة إلى \(q_2=25\)، لدينا \(M_2=4\). ومعكوس \(4\) بترديد \(25\) هو \(19\)، ومن ثم
$$E_2\equiv 4\cdot 19=76 \pmod{100}.$$
الجذر الابتدائي هو
$$x_0=99.$$
تبديل الاختيار على العامل \(4\) يعطي
$$99+2E_1\equiv 99+50\equiv 49 \pmod{100}.$$
وتبديل الاختيار على العامل \(25\) يعطي
$$99+2E_2\equiv 99+152\equiv 51 \pmod{100}.$$
أما تبديل العاملين معًا فيعطي
$$99+2E_1+2E_2\equiv 1 \pmod{100}.$$
إذن الجذور الأربعة هي
$$1,\ 49,\ 51,\ 99.$$
وأكبر جذر أصغر من \(99\) هو \(51\)، وبالتالي
$$I(100)=51.$$
تقوم تطبيقات C++ وPython وJava أولًا ببناء غربال لأصغر عامل أولي حتى الحد الأعلى العام، وبهذا يمكن تحليل كل \(n\) سريعًا إلى قوى أولية.
ولقيمة ثابتة من \(n\)، يبدأ التنفيذ من الباقي \(n-1\)، ثم يحلل \(n\) ويكوّن لكل قوة أولية المسقط CRT المذكور أعلاه. ويُحسب المعكوس المعياري اللازم لهذا المسقط بواسطة خوارزمية إقليدس الموسعة.
كل قوة أولية فردية، وكذلك العامل \(4\)، تضيف فرعًا واحدًا جديدًا ناتجًا من تبديل الجذر المحلي \(-1\) إلى \(+1\). وإذا ظهر عامل \(2^e\) مع \(e\ge 3\)، فإنه يضيف فرعين آخرين ناتجين من الجذرين \(2^{e-1}\pm 1\). وهكذا تنمو قائمة المرشحين وفق عدد الاحتمالات المحلية فقط، ويكون أكبر مرشح أصغر من \(n-1\) هو القيمة المطلوبة \(I(n)\).
وبما أن الخوارزمية تبني مجموعة الحلول الكاملة للمعادلة \(x^2\equiv 1 \pmod n\)، فهي خوارزمية دقيقة تمامًا وليست تقريبية.
إذا كان \(N=2\cdot 10^7\)، فإن بناء غربال أصغر عامل أولي يحتاج إلى زمن \(O(N\log\log N)\) وذاكرة \(O(N)\).
أما لكل قيمة منفردة من \(n\)، فإن استخراج القوى الأولية يكون شبه لوغاريتمي عمليًا، وتكلفة تعداد المرشحين تتناسب مع عدد جذور المعادلة \(x^2\equiv 1 \pmod n\) لا مع \(n\) نفسه. وهذا العدد يعتمد فقط على عدد العوامل الأولية المميزة في تحليل \(n\)، ولذلك يبقى صغيرًا جدًا مقارنة بـ \(n\). ومن هنا تأتي قابلية تنفيذ المسألة حتى الحد \(2\cdot 10^7\).