For each positive integer \(k\), define
$$N_k=4k^2+1,$$
and let \(P(k)\) denote the largest prime factor of \(N_k\). The goal is to compute
$$S(K)=\sum_{k=1}^{K} P(k)\pmod{10^{18}},\qquad K=10^7.$$
Factoring each \(N_k\) independently would be far too slow, so the efficient approach is to sieve prime divisors across all \(k\) at once.
The key observation is that a prime divisor of \(4k^2+1\) forces \(k\) into very specific residue classes. That turns the problem into a structured sieve on arithmetic progressions rather than \(10^7\) separate factorizations.
If a prime \(p\) divides \(4k^2+1\), then
$$4k^2\equiv -1\pmod p,$$
so \(-1\) must be a quadratic residue modulo \(p\). Since \(4k^2+1\) is always odd, \(p=2\) never occurs. For an odd prime, the existence of a solution to
$$x^2\equiv -1\pmod p$$
implies
$$p\equiv 1\pmod 4.$$
Therefore only primes congruent to \(1 \bmod 4\) can appear as prime divisors of the numbers in this family.
Fix an odd prime \(p\equiv 1\pmod 4\). Once we find a square root \(x\) of \(-1\) modulo \(p\), we have
$$x^2\equiv -1\pmod p.$$
Because \(4k^2=(2k)^2\), the divisibility condition \(p\mid 4k^2+1\) is equivalent to
$$2k\equiv \pm x \pmod p.$$
Since \(2\) is invertible modulo every odd prime, this becomes
$$k\equiv \pm x\cdot 2^{-1}\pmod p.$$
So each eligible prime contributes at most two arithmetic progressions in \(k\). The implementations use Tonelli-Shanks to obtain the modular square root and then sweep exactly those two progressions.
When a progression reaches an index \(k\), the current value attached to that index is checked for divisibility by \(p\). If \(p\) divides it, the algorithm divides by \(p\) repeatedly until the factor is gone:
$$N_k=p^{v_p(N_k)}\cdot M_k,\qquad p\nmid M_k.$$
That repeated division matters because values such as \(325=5^2\cdot 13\) contain higher powers. Processing the complete \(p\)-adic valuation ensures the remaining cofactor is always correct after the sweep for \(p\).
The prime sweep runs in increasing order. For each \(k\), the implementation stores two pieces of information:
$$\text{current remainder of }N_k,\qquad \text{largest processed prime that has divided }N_k.$$
Because the primes are handled from small to large, whenever a new prime divides the current remainder it automatically becomes the largest processed prime seen so far for that \(k\). After all eligible small primes have been removed, the largest prime factor is simply the larger of those two stored quantities.
The implementations generate primes only up to \(2K\). This is sufficient because every prime factor \(p\le 2K\) of any \(N_k\) is removed during the sieve. After that, the remaining cofactor has no prime divisor \(\le 2K\).
Suppose the leftover cofactor were composite and greater than \(1\). Then it would have at least two prime factors, both strictly larger than \(2K\), so it would be at least
$$(2K+1)^2=4K^2+4K+1,$$
which is larger than every value in the family because
$$4k^2+1\le 4K^2+1.$$
This contradiction shows that the leftover is either \(1\) or one prime. Hence the final answer for each \(k\) is
$$P(k)=\max(\text{largest processed prime for }k,\ \text{leftover cofactor for }k).$$
For \(K=10\), the values are small enough to inspect directly:
$$\begin{aligned} 4(1)^2+1&=5 &&\Rightarrow P(1)=5,\\ 4(2)^2+1&=17 &&\Rightarrow P(2)=17,\\ 4(3)^2+1&=37 &&\Rightarrow P(3)=37,\\ 4(4)^2+1&=65=5\cdot 13 &&\Rightarrow P(4)=13,\\ 4(5)^2+1&=101 &&\Rightarrow P(5)=101,\\ 4(6)^2+1&=145=5\cdot 29 &&\Rightarrow P(6)=29,\\ 4(7)^2+1&=197 &&\Rightarrow P(7)=197,\\ 4(8)^2+1&=257 &&\Rightarrow P(8)=257,\\ 4(9)^2+1&=325=5^2\cdot 13 &&\Rightarrow P(9)=13,\\ 4(10)^2+1&=401 &&\Rightarrow P(10)=401. \end{aligned}$$
Therefore
$$S(10)=5+17+37+13+101+29+197+257+13+401=1070.$$
This is also a useful checkpoint for the implementation. The sieve interpretation is visible here: the prime \(5\) hits the residue classes \(k\equiv 1,4\pmod 5\), so among the first ten indices it strips factors from \(k=1,4,6,9\).
The C++, Python, and Java implementations all follow the same pipeline. First they build all primes up to \(2K\) with a linear sieve. Then they initialize an array containing the values \(4k^2+1\) and a second array that stores the largest prime factor already confirmed during the sweep.
For each prime \(p\equiv 1\pmod 4\), the implementation computes a square root of \(-1\) modulo \(p\), converts that root into the two valid residue classes of \(k\), and walks those arithmetic progressions with step \(p\). At each visited index, if the current remainder is divisible by \(p\), the code divides out all powers of \(p\) and updates the stored largest processed prime.
After every eligible prime up to \(2K\) has been handled, each index \(k\) has a remaining cofactor that is either \(1\) or prime. The implementation takes the maximum of that leftover value and the largest processed prime, adds it to the running total, and reduces the sum modulo \(10^{18}\).
Generating all primes up to \(2K\) with a linear sieve costs \(O(K)\) time and \(O(K)\) memory. The residue-class sweeps contribute about
$$\sum_{\substack{p\le 2K\\ p\equiv 1\!\!\!\pmod 4}} \frac{2K}{p},$$
which has the usual near-harmonic sieve growth, so the overall running time is \(O(K\log\log K)\) in practice. The repeated divisions by prime powers do not change that asymptotic bound, and the dominant storage remains the two arrays of length \(K+1\), so the memory usage is \(O(K)\).
Für jede positive ganze Zahl \(k\) definieren wir
$$N_k=4k^2+1,$$
und \(P(k)\) sei der größte Primfaktor von \(N_k\). Gesucht ist
$$S(K)=\sum_{k=1}^{K} P(k)\pmod{10^{18}},\qquad K=10^7.$$
Eine separate Faktorisierung jedes einzelnen \(N_k\) wäre viel zu langsam. Deshalb nutzt die Lösung ein Sieb über alle \(k\), das gemeinsame Kongruenzklassen der Primteiler ausnutzt.
Die zentrale Beobachtung lautet: Wenn eine Primzahl \(p\) einen Term \(4k^2+1\) teilt, dann ist \(k\) modulo \(p\) stark eingeschränkt. Dadurch wird aus vielen Einzel-Faktorisierungen ein strukturiertes Sieb auf arithmetischen Progressionen.
Gilt \(p\mid 4k^2+1\), dann folgt
$$4k^2\equiv -1\pmod p.$$
Also muss \(-1\) ein quadratischer Rest modulo \(p\) sein. Da \(4k^2+1\) stets ungerade ist, kommt \(p=2\) nie vor. Für eine ungerade Primzahl bedeutet die Lösbarkeit von
$$x^2\equiv -1\pmod p$$
genau, dass
$$p\equiv 1\pmod 4.$$
Somit können nur Primzahlen der Form \(4m+1\) als Primteiler dieser Zahlenfamilie auftreten.
Fixiere nun eine Primzahl \(p\equiv 1\pmod 4\). Sobald eine Quadratwurzel \(x\) von \(-1\) modulo \(p\) bekannt ist, gilt
$$x^2\equiv -1\pmod p.$$
Wegen \(4k^2=(2k)^2\) ist die Bedingung \(p\mid 4k^2+1\) äquivalent zu
$$2k\equiv \pm x \pmod p.$$
Da \(2\) modulo jeder ungeraden Primzahl invertierbar ist, folgt
$$k\equiv \pm x\cdot 2^{-1}\pmod p.$$
Für jedes zulässige \(p\) gibt es also höchstens zwei arithmetische Progressionen in \(k\), auf denen \(p\) als Teiler auftreten kann. Die Implementierungen bestimmen die Wurzel mit Tonelli-Shanks und durchlaufen anschließend genau diese beiden Progressionen.
Trifft das Sieb einen Index \(k\), wird geprüft, ob der aktuelle Rest dort durch \(p\) teilbar ist. Falls ja, wird wiederholt durch \(p\) dividiert, bis kein Faktor \(p\) mehr übrig ist:
$$N_k=p^{v_p(N_k)}\cdot M_k,\qquad p\nmid M_k.$$
Das ist wichtig, weil manche Werte höhere Potenzen enthalten, etwa \(325=5^2\cdot 13\). Nur durch vollständiges Entfernen der \(p\)-adischen Bewertung bleibt der gespeicherte Rest korrekt.
Die Primzahlen werden aufsteigend verarbeitet. Für jedes \(k\) speichert die Implementierung daher zwei Informationen:
$$\text{den aktuellen Rest von }N_k,\qquad \text{die größte bisher entfernte Primzahl.}$$
Wenn eine neue Primzahl \(p\) den aktuellen Rest teilt, ist sie wegen der aufsteigenden Reihenfolge automatisch die größte bisher verarbeitete Primzahl für dieses \(k\). Nach dem gesamten Sieb ist der größte Primfaktor also einfach das Maximum aus diesen beiden gespeicherten Größen.
Die Implementierungen erzeugen nur Primzahlen bis \(2K\). Das reicht aus, weil jede Primzahl \(p\le 2K\), die irgendein \(N_k\) teilt, während des Siebs vollständig entfernt wird. Danach besitzt der Rest also keinen Primteiler \(\le 2K\).
Angenommen, der verbleibende Rest wäre dennoch zusammengesetzt und größer als \(1\). Dann hätte er mindestens zwei Primfaktoren, beide strikt größer als \(2K\), also wäre er mindestens
$$(2K+1)^2=4K^2+4K+1,$$
doch für alle betrachteten Werte gilt
$$4k^2+1\le 4K^2+1.$$
Das ist unmöglich. Der verbleibende Rest ist daher entweder \(1\) oder selbst eine Primzahl. Somit gilt
$$P(k)=\max(\text{größte verarbeitete Primzahl für }k,\ \text{Restkoeffaktor für }k).$$
Für \(K=10\) kann man alles direkt ausrechnen:
$$\begin{aligned} 4(1)^2+1&=5 &&\Rightarrow P(1)=5,\\ 4(2)^2+1&=17 &&\Rightarrow P(2)=17,\\ 4(3)^2+1&=37 &&\Rightarrow P(3)=37,\\ 4(4)^2+1&=65=5\cdot 13 &&\Rightarrow P(4)=13,\\ 4(5)^2+1&=101 &&\Rightarrow P(5)=101,\\ 4(6)^2+1&=145=5\cdot 29 &&\Rightarrow P(6)=29,\\ 4(7)^2+1&=197 &&\Rightarrow P(7)=197,\\ 4(8)^2+1&=257 &&\Rightarrow P(8)=257,\\ 4(9)^2+1&=325=5^2\cdot 13 &&\Rightarrow P(9)=13,\\ 4(10)^2+1&=401 &&\Rightarrow P(10)=401. \end{aligned}$$
Daraus folgt
$$S(10)=5+17+37+13+101+29+197+257+13+401=1070.$$
Das ist zugleich ein guter Kontrollwert. Man sieht hier auch das Siebprinzip: Die Primzahl \(5\) trifft die Restklassen \(k\equiv 1,4\pmod 5\) und entfernt daher in den ersten zehn Indizes Faktoren bei \(k=1,4,6,9\).
Die C++-, Python- und Java-Implementierungen verwenden denselben Ablauf. Zuerst erzeugen sie mit einem linearen Sieb alle Primzahlen bis \(2K\). Anschließend wird ein Array mit den Werten \(4k^2+1\) aufgebaut sowie ein zweites Array, das den bisher sicher größten gefundenen Primfaktor pro Index speichert.
Für jede Primzahl \(p\equiv 1\pmod 4\) berechnet die Implementierung eine Quadratwurzel von \(-1\) modulo \(p\), leitet daraus die beiden gültigen Restklassen für \(k\) ab und durchläuft diese Progressionen mit Schrittweite \(p\). Ist der aktuelle Rest an einer solchen Stelle durch \(p\) teilbar, werden alle Potenzen von \(p\) entfernt und der gespeicherte größte verarbeitete Primteiler aktualisiert.
Nachdem alle zulässigen Primzahlen bis \(2K\) abgearbeitet sind, ist der verbleibende Rest an jedem Index entweder \(1\) oder prim. Die Implementierung nimmt das Maximum aus Rest und größtem bereits verarbeiteten Primteiler, addiert diesen Wert zur laufenden Summe und reduziert fortlaufend modulo \(10^{18}\).
Das Erzeugen aller Primzahlen bis \(2K\) mit einem linearen Sieb kostet \(O(K)\) Zeit und \(O(K)\) Speicher. Die Progressions-Sweeps tragen ungefähr
$$\sum_{\substack{p\le 2K\\ p\equiv 1\!\!\!\pmod 4}} \frac{2K}{p}$$
bei und zeigen damit das übliche nahezu harmonische Siebverhalten, also insgesamt \(O(K\log\log K)\) Laufzeit in der Praxis. Das wiederholte Herausdividieren von Primzahlpotenzen ändert diese Größenordnung nicht. Dominant bleiben die beiden Arrays der Länge \(K+1\), also \(O(K)\) Speicher.
Her pozitif \(k\) için
$$N_k=4k^2+1$$
tanımlanır ve \(P(k)\), bu sayının en büyük asal çarpanı olsun. İstenen büyüklük
$$S(K)=\sum_{k=1}^{K} P(k)\pmod{10^{18}},\qquad K=10^7.$$
Her \(N_k\)'yi bağımsız biçimde çarpanlara ayırmak çok yavaş olur. Verimli çözüm bunun yerine asal çarpanları tüm \(k\) değerleri üzerinde aynı anda eleyen bir yaklaşım kullanır.
Ana gözlem şudur: \(4k^2+1\)'i bölen bir asal \(p\), \(k\) için çok belirli artık sınıfları zorunlu kılar. Böylece problem, tek tek faktorizasyonlardan çıkıp aritmetik ilerlemeler üzerinde çalışan düzenli bir eleğe dönüşür.
Eğer bir asal \(p\), \(4k^2+1\)'i bölüyorsa
$$4k^2\equiv -1\pmod p$$
olur. Demek ki \(-1\), mod \(p\) altında kuadratik kalıntı olmalıdır. \(4k^2+1\) her zaman tek olduğu için \(p=2\) hiç görünmez. Tek bir asal için
$$x^2\equiv -1\pmod p$$
denkleminin çözülebilmesi
$$p\equiv 1\pmod 4$$
koşulunu verir. O halde bu sayı ailesinde asal çarpan olarak yalnızca \(4m+1\) biçimindeki asallar yer alabilir.
Şimdi \(p\equiv 1\pmod 4\) olan tek bir asal sabitleyelim. Eğer \(-1\)'in mod \(p\) altında bir karekökü \(x\) ise
$$x^2\equiv -1\pmod p$$
yazabiliriz. \(4k^2=(2k)^2\) olduğundan, \(p\mid 4k^2+1\) koşulu
$$2k\equiv \pm x \pmod p$$
ile aynıdır. \(2\), her tek asal modunda terslenebilir olduğu için
$$k\equiv \pm x\cdot 2^{-1}\pmod p$$
elde edilir. Yani her uygun asal için \(k\) en fazla iki aritmetik ilerleme üzerinde incelenmelidir. Uygulama, bu karekökü Tonelli-Shanks ile bulur ve sonra yalnızca bu iki ilerlemeyi tarar.
Eleğin ulaştığı bir \(k\) değeri için, o indiste tutulan güncel sayı \(p\)'ye bölünebiliyorsa algoritma tekrar tekrar bölme yapar:
$$N_k=p^{v_p(N_k)}\cdot M_k,\qquad p\nmid M_k.$$
Bu ayrıntı önemlidir; çünkü bazı terimler daha yüksek üsler içerir. Örneğin \(325=5^2\cdot 13\) için yalnızca tek bir \(5\) çıkarmak yeterli olmaz. Tüm \(p\)-adik kuvvetin temizlenmesi, kalan kofaktörün her aşamada doğru kalmasını sağlar.
Asallar küçükten büyüğe işlendiği için, her \(k\) için iki bilgiyi saklamak yeterlidir:
$$\text{\(N_k\)'nin o andaki kalanı},\qquad \text{şimdiye kadar kesinleşmiş en büyük işlenmiş asal çarpan.}$$
Eğer yeni bir asal \(p\) bu kalanı bölüyorsa, sıralama gereği bu \(p\) aynı zamanda o ana kadar görülen en büyük işlenmiş asal olur. Tüm elekten sonra en büyük asal çarpan, bu iki saklanan değerin maksimumudur.
Uygulama sadece \(2K\)'ye kadar olan asalları üretir. Bu yeterlidir; çünkü \(N_k\)'nin \(p\le 2K\) olan her asal böleni elekte bulunduğu anda tamamen çıkarılır. Dolayısıyla sonunda kalan kofaktörün \(\le 2K\) olan bir asal böleni kalmaz.
Şimdi kalan kofaktörün \(1\)'den büyük ve bileşik olduğunu varsayalım. Bu durumda en az iki asal çarpanı olmalı ve bunların her biri \(2K\)'den büyük olmalıdır. O zaman kalan değer en az
$$(2K+1)^2=4K^2+4K+1$$
olur. Oysa ailedeki tüm sayılar için
$$4k^2+1\le 4K^2+1$$
geçerlidir. Bu çelişkidir. Demek ki sonda kalan değer ya \(1\) ya da tek bir asaldır. Bu yüzden
$$P(k)=\max(\text{\(k\) için işlenmiş en büyük asal},\ \text{kalan kofaktör})$$
olur.
\(K=10\) için değerler doğrudan görülebilir:
$$\begin{aligned} 4(1)^2+1&=5 &&\Rightarrow P(1)=5,\\ 4(2)^2+1&=17 &&\Rightarrow P(2)=17,\\ 4(3)^2+1&=37 &&\Rightarrow P(3)=37,\\ 4(4)^2+1&=65=5\cdot 13 &&\Rightarrow P(4)=13,\\ 4(5)^2+1&=101 &&\Rightarrow P(5)=101,\\ 4(6)^2+1&=145=5\cdot 29 &&\Rightarrow P(6)=29,\\ 4(7)^2+1&=197 &&\Rightarrow P(7)=197,\\ 4(8)^2+1&=257 &&\Rightarrow P(8)=257,\\ 4(9)^2+1&=325=5^2\cdot 13 &&\Rightarrow P(9)=13,\\ 4(10)^2+1&=401 &&\Rightarrow P(10)=401. \end{aligned}$$
Böylece
$$S(10)=5+17+37+13+101+29+197+257+13+401=1070$$
elde edilir. Bu aynı zamanda iyi bir kontrol değeridir. Elek mantığı da burada açıkça görünür: \(5\) asalı, \(k\equiv 1,4\pmod 5\) sınıflarına vurur ve ilk on indiste \(k=1,4,6,9\) için çarpan temizler.
C++, Python ve Java uygulamalarının üçü de aynı matematiksel hattı izler. Önce lineer bir asal eleğiyle \(2K\)'ye kadar tüm asallar üretilir. Ardından bir dizide \(4k^2+1\) değerleri tutulur; ikinci bir dizide ise eleme sırasında o indeks için kesinleşmiş en büyük asal çarpan saklanır.
Her \(p\equiv 1\pmod 4\) asalı için uygulama, \(-1\)'in mod \(p\) karekökünü hesaplar, bundan \(k\) için iki geçerli artık sınıfını çıkarır ve adım uzunluğu \(p\) olan bu iki aritmetik ilerlemeyi dolaşır. Ulaşılan indisteki güncel kalan değer \(p\)'ye bölünüyorsa tüm kuvvetler çıkarılır ve saklanan en büyük işlenmiş asal güncellenir.
Tüm uygun asallar işlendiğinde her \(k\) için kalan kofaktör ya \(1\) ya da asaldır. Uygulama, bu kalan ile daha önce saklanan en büyük asalın maksimumunu toplama ekler ve toplamı sürekli \(10^{18}\) modunda tutar.
\(2K\)'ye kadar tüm asalları lineer elekle üretmek \(O(K)\) zaman ve \(O(K)\) bellek gerektirir. Artık sınıfları üzerinde yapılan taramalar yaklaşık olarak
$$\sum_{\substack{p\le 2K\\ p\equiv 1\!\!\!\pmod 4}} \frac{2K}{p}$$
büyüklüğünde iş üretir; bu da klasik elek davranışıyla pratikte \(O(K\log\log K)\) zamana karşılık gelir. Asal kuvvetleri tekrar tekrar bölerek temizlemek bu mertebeyi değiştirmez. Baskın bellek kullanımı uzunluğu \(K+1\) olan iki dizi olduğundan toplam bellek \(O(K)\)'dir.
Para cada entero positivo \(k\), definimos
$$N_k=4k^2+1,$$
y sea \(P(k)\) el mayor factor primo de \(N_k\). Hay que calcular
$$S(K)=\sum_{k=1}^{K} P(k)\pmod{10^{18}},\qquad K=10^7.$$
Factorizar cada \(N_k\) por separado sería demasiado lento. La idea eficiente es cribar divisores primos sobre todos los valores de \(k\) al mismo tiempo.
La observación clave es que, si un primo divide a \(4k^2+1\), entonces \(k\) debe pertenecer a clases de residuos muy concretas. Eso convierte el problema en una criba sobre progresiones aritméticas en lugar de millones de factorizaciones independientes.
Si un primo \(p\) divide \(4k^2+1\), entonces
$$4k^2\equiv -1\pmod p.$$
Por tanto, \(-1\) debe ser un residuo cuadrático módulo \(p\). Como \(4k^2+1\) siempre es impar, el primo \(2\) nunca aparece. Para un primo impar, que exista solución de
$$x^2\equiv -1\pmod p$$
implica
$$p\equiv 1\pmod 4.$$
Así, dentro de esta familia de números solo pueden aparecer como factores primos los congruentes con \(1 \bmod 4\).
Fijemos un primo \(p\equiv 1\pmod 4\). Si \(x\) es una raíz cuadrada de \(-1\) módulo \(p\), entonces
$$x^2\equiv -1\pmod p.$$
Como \(4k^2=(2k)^2\), la condición \(p\mid 4k^2+1\) equivale a
$$2k\equiv \pm x \pmod p.$$
Y como \(2\) es invertible módulo todo primo impar, obtenemos
$$k\equiv \pm x\cdot 2^{-1}\pmod p.$$
Es decir, cada primo admisible genera como mucho dos progresiones aritméticas en \(k\). Las implementaciones usan Tonelli-Shanks para hallar la raíz modular y luego recorren exactamente esas dos progresiones.
Cuando la criba alcanza un índice \(k\), comprueba si el valor actual asociado a ese índice es divisible por \(p\). Si lo es, divide repetidamente por \(p\) hasta eliminar toda la potencia:
$$N_k=p^{v_p(N_k)}\cdot M_k,\qquad p\nmid M_k.$$
Esto es importante porque algunos términos contienen potencias altas, por ejemplo \(325=5^2\cdot 13\). Si solo se eliminara un factor \(p\), el resto almacenado quedaría incorrecto.
Los primos se procesan en orden creciente. Por eso, para cada \(k\), la implementación solo necesita guardar dos datos:
$$\text{el resto actual de }N_k,\qquad \text{el mayor primo procesado que ya ha dividido }N_k.$$
Cuando un nuevo primo divide el resto actual, automáticamente pasa a ser el mayor primo pequeño confirmado para ese índice. Al final de la criba, el mayor factor primo se obtiene tomando el máximo entre esas dos cantidades.
Las implementaciones generan primos solo hasta \(2K\). Eso basta porque todo primo \(p\le 2K\) que divide algún \(N_k\) queda eliminado durante la criba. Por tanto, el cofactor final ya no tiene divisores primos \(\le 2K\).
Supongamos que ese cofactor restante fuera compuesto y mayor que \(1\). Entonces tendría al menos dos factores primos, ambos estrictamente mayores que \(2K\), así que valdría al menos
$$(2K+1)^2=4K^2+4K+1,$$
mientras que para toda la familia se cumple
$$4k^2+1\le 4K^2+1.$$
Es imposible. Luego el resto final es \(1\) o bien un primo. En consecuencia,
$$P(k)=\max(\text{mayor primo procesado para }k,\ \text{cofactor restante para }k).$$
Para \(K=10\), se puede ver todo directamente:
$$\begin{aligned} 4(1)^2+1&=5 &&\Rightarrow P(1)=5,\\ 4(2)^2+1&=17 &&\Rightarrow P(2)=17,\\ 4(3)^2+1&=37 &&\Rightarrow P(3)=37,\\ 4(4)^2+1&=65=5\cdot 13 &&\Rightarrow P(4)=13,\\ 4(5)^2+1&=101 &&\Rightarrow P(5)=101,\\ 4(6)^2+1&=145=5\cdot 29 &&\Rightarrow P(6)=29,\\ 4(7)^2+1&=197 &&\Rightarrow P(7)=197,\\ 4(8)^2+1&=257 &&\Rightarrow P(8)=257,\\ 4(9)^2+1&=325=5^2\cdot 13 &&\Rightarrow P(9)=13,\\ 4(10)^2+1&=401 &&\Rightarrow P(10)=401. \end{aligned}$$
Por tanto,
$$S(10)=5+17+37+13+101+29+197+257+13+401=1070.$$
Ese valor sirve también como comprobación. Además, aquí se ve bien la lógica de la criba: el primo \(5\) corresponde a las clases \(k\equiv 1,4\pmod 5\), de modo que entre los diez primeros índices afecta a \(k=1,4,6,9\).
Las implementaciones en C++, Python y Java siguen exactamente la misma idea matemática. Primero generan todos los primos hasta \(2K\) con una criba lineal. Después inicializan un arreglo con los valores \(4k^2+1\) y un segundo arreglo que guarda, para cada índice, el mayor factor primo ya confirmado durante el proceso.
Para cada primo \(p\equiv 1\pmod 4\), la implementación calcula una raíz cuadrada de \(-1\) módulo \(p\), la transforma en las dos clases de residuos válidas para \(k\) y recorre ambas progresiones aritméticas con paso \(p\). Cuando el resto actual en una posición es divisible por \(p\), divide todas las potencias de \(p\) y actualiza el mayor primo procesado para ese índice.
Al terminar con todos los primos admisibles hasta \(2K\), el cofactor que queda en cada posición es \(1\) o primo. Entonces la implementación toma el máximo entre ese resto y el mayor primo ya registrado, lo suma al acumulador y mantiene la suma reducida módulo \(10^{18}\).
Generar todos los primos hasta \(2K\) con una criba lineal cuesta \(O(K)\) tiempo y \(O(K)\) memoria. Los recorridos por clases de residuos aportan aproximadamente
$$\sum_{\substack{p\le 2K\\ p\equiv 1\!\!\!\pmod 4}} \frac{2K}{p},$$
que es el crecimiento casi armónico típico de una criba, así que el tiempo total es \(O(K\log\log K)\) en la práctica. Las divisiones repetidas por potencias primas no cambian ese orden asintótico. La memoria dominante sigue siendo la de los dos arreglos de longitud \(K+1\), por lo que el uso total es \(O(K)\).
对每个正整数 \(k\)，定义
$$N_k=4k^2+1,$$
并记 \(P(k)\) 为 \(N_k\) 的最大素因子。题目要求计算
$$S(K)=\sum_{k=1}^{K} P(k)\pmod{10^{18}},\qquad K=10^7.$$
如果把每个 \(N_k\) 都单独做一次完整分解，代价会非常高。真正可行的方法，是把“哪些素数会整除哪些 \(k\)”统一组织起来，像筛法一样批量处理。
核心思想是：一旦某个素数 \(p\) 整除 \(4k^2+1\)，那么 \(k\) 在模 \(p\) 意义下就只能落在极少数几个剩余类里。这样一来，问题就不再是 \(10^7\) 次互不相关的分解，而是若干条等差序列上的系统性筛除。
若素数 \(p\mid 4k^2+1\)，则有
$$4k^2\equiv -1\pmod p.$$
这说明 \(-1\) 必须是模 \(p\) 下的二次剩余。由于 \(4k^2+1\) 永远是奇数，所以 \(p=2\) 根本不可能出现。对奇素数来说，方程
$$x^2\equiv -1\pmod p$$
有解当且仅当
$$p\equiv 1\pmod 4.$$
因此，在整个数列 \(4k^2+1\) 中，可能出现的奇素因子只来自 \(4m+1\) 这一类素数。
固定一个满足 \(p\equiv 1\pmod 4\) 的素数。若 \(x\) 是 \(-1\) 在模 \(p\) 下的一个平方根，那么
$$x^2\equiv -1\pmod p.$$
因为 \(4k^2=(2k)^2\)，条件 \(p\mid 4k^2+1\) 就等价于
$$2k\equiv \pm x\pmod p.$$
又由于 \(2\) 在模奇素数下总是可逆，所以进一步得到
$$k\equiv \pm x\cdot 2^{-1}\pmod p.$$
这意味着每个可用素数 \(p\) 最多只对应 \(k\) 的两个剩余类。实现里使用 Tonelli-Shanks 算法求出 \(-1\) 的模平方根，然后只沿着这两条等差序列去检查和消去该素数。
筛到某个索引 \(k\) 时，如果当前位置保存的当前值可以被 \(p\) 整除，就不能只除一次，而要把所有 \(p\) 的幂都除掉：
$$N_k=p^{v_p(N_k)}\cdot M_k,\qquad p\nmid M_k.$$
这样做非常重要，因为有些项会含有更高的素数幂。例如
$$325=5^2\cdot 13.$$
如果只删掉一个 \(5\)，那么后续保存的剩余部分就会出错；只有完整去掉 \(p\)-进指数，当前余因子才是正确的。
实现按素数从小到大处理。因此，对每个 \(k\)，只需要维护两类信息：
$$\text{当前位置尚未分解完的余因子},\qquad \text{到目前为止已确认的最大处理过的素数因子。}$$
由于筛的顺序是递增的，所以一旦某个新素数 \(p\) 成功整除了当前位置，它就自动成为该位置到目前为止最大的已处理素因子。等所有相关小素数处理结束后，最大素因子只可能在这两者之中。
实现只生成到 \(2K\) 为止的素数，这一点并不是巧合，而是由上界决定的。因为所有 \(\le 2K\) 的相关素因子都会在筛过程中被完整删掉，所以最后留下来的余因子不可能再含有任何 \(\le 2K\) 的素因子。
假设最后余下的值大于 \(1\) 且仍然是合数，那么它至少有两个素因子，并且这两个素因子都必须严格大于 \(2K\)。于是这个合数至少是
$$(2K+1)^2=4K^2+4K+1,$$
但整个数列中的项都满足
$$4k^2+1\le 4K^2+1.$$
这就矛盾了。因此最终剩下的余因子只能是 \(1\) 或者一个单独的素数。于是每个 \(k\) 的最大素因子就是
$$P(k)=\max(\text{该位置已确认的最大处理过的素因子},\ \text{最终余因子}).$$
取 \(K=10\)，可以直接把前十项写出来：
$$\begin{aligned} 4(1)^2+1&=5 &&\Rightarrow P(1)=5,\\ 4(2)^2+1&=17 &&\Rightarrow P(2)=17,\\ 4(3)^2+1&=37 &&\Rightarrow P(3)=37,\\ 4(4)^2+1&=65=5\cdot 13 &&\Rightarrow P(4)=13,\\ 4(5)^2+1&=101 &&\Rightarrow P(5)=101,\\ 4(6)^2+1&=145=5\cdot 29 &&\Rightarrow P(6)=29,\\ 4(7)^2+1&=197 &&\Rightarrow P(7)=197,\\ 4(8)^2+1&=257 &&\Rightarrow P(8)=257,\\ 4(9)^2+1&=325=5^2\cdot 13 &&\Rightarrow P(9)=13,\\ 4(10)^2+1&=401 &&\Rightarrow P(10)=401. \end{aligned}$$
因此
$$S(10)=5+17+37+13+101+29+197+257+13+401=1070.$$
这也是程序可以用来核对的小规模结果。筛法视角在这里非常直观：素数 \(5\) 对应的剩余类是 \(k\equiv 1,4\pmod 5\)，所以它会命中 \(k=1,4,6,9\) 这些位置；同样地，其他可行素数也只沿各自的两条等差序列传播。
C++、Python 和 Java 三个实现采用的是同一条数学路线。首先，它们用线性筛生成所有不超过 \(2K\) 的素数。然后建立一个数组保存每个 \(k\) 对应的值 \(4k^2+1\)，再建立另一个数组保存“到目前为止已经确认的最大素因子”。
接下来，对每个满足 \(p\equiv 1\pmod 4\) 的素数，实现会先求出 \(-1\) 在模 \(p\) 下的一个平方根，再把它转换成 \(k\) 的两个合法剩余类。随后程序按步长 \(p\) 扫描这两条等差序列；凡是当前位置仍被 \(p\) 整除，就把 \(p\) 的所有幂次全部除尽，并把该位置记录的最大已处理素因子更新为 \(p\)。
当所有相关素数都处理完之后，每个位置上剩余的值一定是 \(1\) 或一个素数。最后，程序对每个 \(k\) 取“最终余因子”和“已确认最大素因子”两者中的较大者，累加到总和里，并始终对 \(10^{18}\) 取模。
用线性筛生成不超过 \(2K\) 的全部素数需要 \(O(K)\) 时间和 \(O(K)\) 空间。之后，对所有满足条件的素数做剩余类扫描，其总工作量大致是
$$\sum_{\substack{p\le 2K\\ p\equiv 1\!\!\!\pmod 4}} \frac{2K}{p},$$
也就是典型筛法中的近似调和增长，所以整体时间复杂度在实践中为 \(O(K\log\log K)\)。对某个位置反复除去同一素数的高次幂，并不会改变这个渐近量级。主要内存占用仍然来自两个长度为 \(K+1\) 的数组，因此空间复杂度为 \(O(K)\)。
Для каждого положительного целого \(k\) определим
$$N_k=4k^2+1,$$
а через \(P(k)\) обозначим наибольший простой делитель числа \(N_k\). Требуется вычислить
$$S(K)=\sum_{k=1}^{K} P(k)\pmod{10^{18}},\qquad K=10^7.$$
Факторизовать каждое \(N_k\) отдельно слишком дорого. Эффективное решение организует общий проход по простым делителям сразу для всех значений \(k\).
Главное наблюдение состоит в том, что простой делитель числа \(4k^2+1\) заставляет \(k\) лежать в очень узких классах вычетов. Благодаря этому задача превращается не в миллионы независимых факторизаций, а в решето по арифметическим прогрессиям.
Если простой \(p\) делит \(4k^2+1\), то
$$4k^2\equiv -1\pmod p.$$
Значит, \(-1\) должно быть квадратичным вычетом по модулю \(p\). Так как \(4k^2+1\) всегда нечетно, число \(2\) никогда не встречается среди делителей. Для нечетного простого существование решения уравнения
$$x^2\equiv -1\pmod p$$
эквивалентно условию
$$p\equiv 1\pmod 4.$$
Следовательно, в этой семье чисел могут появляться только простые делители вида \(4m+1\).
Зафиксируем простой \(p\equiv 1\pmod 4\). Если \(x\) является квадратным корнем из \(-1\) по модулю \(p\), то
$$x^2\equiv -1\pmod p.$$
Поскольку \(4k^2=(2k)^2\), условие \(p\mid 4k^2+1\) равносильно
$$2k\equiv \pm x\pmod p.$$
Так как \(2\) обратимо по модулю любого нечетного простого, получаем
$$k\equiv \pm x\cdot 2^{-1}\pmod p.$$
Итак, каждому допустимому простому соответствуют не более двух арифметических прогрессий по \(k\). Реализация находит корень с помощью алгоритма Тонелли-Шенкса и затем проходит ровно по этим двум прогрессиям.
Когда решето доходит до некоторого индекса \(k\), текущий остаток в этой позиции проверяется на делимость на \(p\). Если делимость есть, программа делит на \(p\) до тех пор, пока фактор не исчезнет полностью:
$$N_k=p^{v_p(N_k)}\cdot M_k,\qquad p\nmid M_k.$$
Это важно, потому что некоторые значения содержат высокие степени простых. Например,
$$325=5^2\cdot 13.$$
Если удалить только один множитель \(5\), оставшийся кофактор будет неверным. Полное удаление \(p\)-адической степени поддерживает корректное состояние на каждом шаге.
Простые перебираются по возрастанию. Поэтому для каждого \(k\) достаточно хранить две величины:
$$\text{текущий неразложенный остаток},\qquad \text{наибольший уже обработанный простой делитель.}$$
Когда новый простой \(p\) действительно делит текущий остаток, он автоматически становится крупнейшим среди всех уже обработанных делителей для данного индекса. После завершения просеивания ответ для этого \(k\) равен максимуму из этих двух величин.
Реализации строят список простых только до \(2K\). Этого достаточно, потому что любой простой делитель \(p\le 2K\) будет найден и полностью удален во время соответствующего прохода. Значит, итоговый остаток уже не содержит простых множителей \(\le 2K\).
Предположим теперь, что этот остаток больше \(1\) и при этом составной. Тогда у него есть по меньшей мере два простых множителя, оба строго больше \(2K\). Следовательно, сам остаток не меньше
$$(2K+1)^2=4K^2+4K+1,$$
тогда как для всех членов нашей семьи верно
$$4k^2+1\le 4K^2+1.$$
Получаем противоречие. Значит, оставшийся кофактор либо равен \(1\), либо сам является простым числом. Поэтому
$$P(k)=\max(\text{наибольший обработанный простой для }k,\ \text{оставшийся кофактор}).$$
При \(K=10\) можно выписать все значения явно:
$$\begin{aligned} 4(1)^2+1&=5 &&\Rightarrow P(1)=5,\\ 4(2)^2+1&=17 &&\Rightarrow P(2)=17,\\ 4(3)^2+1&=37 &&\Rightarrow P(3)=37,\\ 4(4)^2+1&=65=5\cdot 13 &&\Rightarrow P(4)=13,\\ 4(5)^2+1&=101 &&\Rightarrow P(5)=101,\\ 4(6)^2+1&=145=5\cdot 29 &&\Rightarrow P(6)=29,\\ 4(7)^2+1&=197 &&\Rightarrow P(7)=197,\\ 4(8)^2+1&=257 &&\Rightarrow P(8)=257,\\ 4(9)^2+1&=325=5^2\cdot 13 &&\Rightarrow P(9)=13,\\ 4(10)^2+1&=401 &&\Rightarrow P(10)=401. \end{aligned}$$
Следовательно,
$$S(10)=5+17+37+13+101+29+197+257+13+401=1070.$$
Это удобная контрольная точка. Здесь же видно, как работает решето: простой \(5\) соответствует классам вычетов \(k\equiv 1,4\pmod 5\), поэтому среди первых десяти индексов он срабатывает при \(k=1,4,6,9\).
Реализации на C++, Python и Java используют одну и ту же схему. Сначала линейным решетом строятся все простые числа до \(2K\). Затем создается массив значений \(4k^2+1\) и второй массив, в котором для каждого индекса хранится наибольший простой делитель, уже подтвержденный в ходе проходов.
Для каждого простого \(p\equiv 1\pmod 4\) программа находит квадратный корень из \(-1\) по модулю \(p\), переводит его в два допустимых класса вычетов для \(k\) и идет по этим арифметическим прогрессиям с шагом \(p\). Если текущий остаток в позиции делится на \(p\), удаляются все степени \(p\), а сохраненный максимальный уже обработанный простой обновляется.
После обработки всех допустимых простых до \(2K\) оставшийся кофактор в каждой позиции равен \(1\) или простому числу. Затем реализация берет максимум между этим остатком и ранее сохраненным простым, прибавляет его к сумме и постоянно редуцирует ответ по модулю \(10^{18}\).
Построение всех простых до \(2K\) линейным решетом занимает \(O(K)\) времени и \(O(K)\) памяти. Проходы по классам вычетов дают суммарно примерно
$$\sum_{\substack{p\le 2K\\ p\equiv 1\!\!\!\pmod 4}} \frac{2K}{p},$$
что соответствует обычному почти гармоническому росту ситовых методов, то есть общей практической сложности \(O(K\log\log K)\). Повторное деление на степени простых не меняет этот порядок. Основная память по-прежнему приходится на два массива длины \(K+1\), поэтому итоговая пространственная сложность равна \(O(K)\).
لكل عدد صحيح موجب \(k\) نعرّف
$$N_k=4k^2+1,$$
ولنرمز بـ \(P(k)\) إلى أكبر عامل أولي للعدد \(N_k\). المطلوب هو حساب
$$S(K)=\sum_{k=1}^{K} P(k)\pmod{10^{18}},\qquad K=10^7.$$
تحليل كل قيمة \(N_k\) على حدة سيكون بطيئًا جدًا. لذلك تعتمد الفكرة الفعالة على غربلة العوامل الأولية عبر جميع قيم \(k\) معًا.
الفكرة الجوهرية هي أن وجود عدد أولي \(p\) يقسم \(4k^2+1\) يفرض على \(k\) أن يقع في أصناف توافقية محددة جدًا. وهكذا تتحول المسألة من ملايين التحليلات المستقلة إلى غربال يعمل على متتاليات حسابية.
إذا كان \(p\) عددًا أوليًا ويقسم \(4k^2+1\)، فلدينا
$$4k^2\equiv -1\pmod p.$$
وهذا يعني أن \(-1\) يجب أن يكون باقيًا تربيعيًا modulo \(p\). وبما أن \(4k^2+1\) عدد فردي دائمًا، فإن \(p=2\) لا يظهر أبدًا. أما للعدد الأولي الفردي، فإن قابلية حل المعادلة
$$x^2\equiv -1\pmod p$$
تعطي الشرط
$$p\equiv 1\pmod 4.$$
إذن العوامل الأولية الممكنة في هذه العائلة لا تكون إلا من الشكل \(4m+1\).
لنثبّت عددًا أوليًا \(p\equiv 1\pmod 4\). إذا كان \(x\) جذرًا تربيعيًا لـ \(-1\) modulo \(p\)، فإن
$$x^2\equiv -1\pmod p.$$
وبما أن \(4k^2=(2k)^2\)، فإن الشرط \(p\mid 4k^2+1\) يكافئ
$$2k\equiv \pm x \pmod p.$$
ولأن \(2\) قابل للعكس modulo كل عدد أولي فردي، نحصل على
$$k\equiv \pm x\cdot 2^{-1}\pmod p.$$
أي أن كل عدد أولي مسموح به يحدد على الأكثر متتاليتين حسابيتين في قيم \(k\). تستخدم التطبيقات خوارزمية Tonelli-Shanks لاستخراج الجذر المعياري ثم تمسح هاتين المتتاليتين فقط.
عندما يصل الغربال إلى موضع \(k\)، فإنه يفحص القيمة الحالية المخزنة هناك. فإذا كانت قابلة للقسمة على \(p\)، فإنه يقسمها مرارًا حتى يزيل جميع قوى \(p\):
$$N_k=p^{v_p(N_k)}\cdot M_k,\qquad p\nmid M_k.$$
هذا التفصيل مهم لأن بعض الحدود تحتوي على قوى أعلى، مثل
$$325=5^2\cdot 13.$$
ولو أزيل عامل واحد فقط من \(5\) لبقيت القيمة المخزنة غير صحيحة. إزالة التقييم \(p\)-العددي كاملًا هي التي تحافظ على صحة العامل المتبقي.
تُعالَج الأعداد الأولية بترتيب تصاعدي، ولذلك يكفي أن يحفظ التنفيذ لكل \(k\) مقدارين:
نرمز إلى العامل المتبقي الحالي بـ \(R_k\)، ونرمز إلى أكبر عدد أولي مُعالَج ثبت أنه يقسم \(N_k\) بـ \(L_k\).
فإذا قسَمَ عدد أولي جديد \(p\) القيمة الحالية، فإنه يصبح تلقائيًا أكبر عدد أولي مُعالَج لهذا الموضع بسبب ترتيب المرور. وبعد انتهاء الغربلة، يكون أكبر عامل أولي هو الأكبر بين هاتين القيمتين.
التنفيذ لا يولد إلا الأعداد الأولية حتى \(2K\)، وهذا كافٍ تمامًا. فكل عامل أولي \(p\le 2K\) يقسم أحد الحدود سيُزال أثناء المسح الموافق له، ولذلك لا يمكن أن يحتوي العامل المتبقي النهائي على أي عامل أولي \(\le 2K\).
افترض أن العامل المتبقي النهائي أكبر من \(1\) ومع ذلك مركب. عندئذٍ لا بد أن يملك عاملين أوليين على الأقل، وكلاهما أكبر من \(2K\)، ومن ثم يكون هذا العدد على الأقل
$$(2K+1)^2=4K^2+4K+1,$$
لكن كل عنصر في العائلة يحقق
$$4k^2+1\le 4K^2+1.$$
وهذا تناقض. إذن العامل المتبقي النهائي إما \(1\) أو عدد أولي واحد. لذلك نحصل على
$$P(k)=\max(L_k,R_k).$$
عند \(K=10\) يمكن كتابة القيم مباشرة:
$$\begin{aligned} 4(1)^2+1&=5 &&\Rightarrow P(1)=5,\\ 4(2)^2+1&=17 &&\Rightarrow P(2)=17,\\ 4(3)^2+1&=37 &&\Rightarrow P(3)=37,\\ 4(4)^2+1&=65=5\cdot 13 &&\Rightarrow P(4)=13,\\ 4(5)^2+1&=101 &&\Rightarrow P(5)=101,\\ 4(6)^2+1&=145=5\cdot 29 &&\Rightarrow P(6)=29,\\ 4(7)^2+1&=197 &&\Rightarrow P(7)=197,\\ 4(8)^2+1&=257 &&\Rightarrow P(8)=257,\\ 4(9)^2+1&=325=5^2\cdot 13 &&\Rightarrow P(9)=13,\\ 4(10)^2+1&=401 &&\Rightarrow P(10)=401. \end{aligned}$$
ومن ثم
$$S(10)=5+17+37+13+101+29+197+257+13+401=1070.$$
وهذه قيمة فحص جيدة. كما يظهر منها منطق الغربال بوضوح: فالعدد الأولي \(5\) يضرب الأصناف \(k\equiv 1,4\pmod 5\)، ولذلك يزيل عوامل عند \(k=1,4,6,9\) ضمن أول عشرة مواضع.
تتبع تطبيقات C++ وPython وJava الخطة الرياضية نفسها. أولًا تُولَّد جميع الأعداد الأولية حتى \(2K\) باستخدام غربال خطي. ثم يُنشأ مصفوفة تحتوي القيم \(4k^2+1\)، ومصفوفة ثانية تحفظ لكل موضع أكبر عامل أولي تأكد وجوده أثناء المرور.
بعد ذلك، لكل عدد أولي \(p\equiv 1\pmod 4\)، يحسب التنفيذ جذرًا تربيعيًا لـ \(-1\) modulo \(p\)، ويحوّله إلى فئتي توافق صالحـتين لـ \(k\)، ثم يسير على هاتين المتتاليتين الحسابيتين بخطوة \(p\). وإذا كانت القيمة الحالية في موضع ما قابلة للقسمة على \(p\)، فإنه يزيل جميع قوى \(p\) ويحدّث أكبر عدد أولي مُعالَج لذلك الموضع.
وبعد الانتهاء من جميع الأعداد الأولية المناسبة حتى \(2K\)، يصبح العامل المتبقي في كل موضع مساويًا لـ \(1\) أو لعدد أولي. ثم يأخذ التنفيذ القيمة الكبرى بين هذا الباقي وبين أكبر عامل أولي سبق تسجيله، ويضيفها إلى المجموع مع الاختزال المستمر modulo \(10^{18}\).
توليد جميع الأعداد الأولية حتى \(2K\) بغربال خطي يحتاج إلى \(O(K)\) زمنًا و\(O(K)\) ذاكرة. أما المرور على أصناف التوافق فيضيف عملًا يقارب
$$\sum_{\substack{p\le 2K\\ p\equiv 1\!\!\!\pmod 4}} \frac{2K}{p},$$
وهو السلوك شبه التوافقي المعتاد في مسائل الغربلة، ولذلك تكون الكلفة الكلية عمليًا \(O(K\log\log K)\). كما أن إزالة قوى العوامل الأولية بالتكرار لا تغيّر هذا الترتيب. ويبقى الاستهلاك الأكبر للذاكرة هو مصفوفتين بطول \(K+1\)، أي \(O(K)\).