We must count the prime values of
$$Q_n=2n^2+2n+1=n^2+(n+1)^2$$
under the bound \(Q_n \lt \text{LIMIT}=5\cdot 10^{15}\). Testing each value independently would be far too slow, so the solution uses a sieve specialized to this quadratic family.
The sequence is strictly increasing for \(n\ge 1\), so it is enough to search indices
$$1\le n\le m,\qquad m=\max\{n: 2n^2+2n+1\lt \text{LIMIT}\}.$$
Solving the quadratic inequality gives the approximation
$$m\approx \left\lfloor \frac{\sqrt{2\cdot \text{LIMIT}-1}-1}{2}\right\rfloor,$$
and the code then adjusts by a couple of integer checks to respect the strict inequality.
Assume \(p=Q_i\) is prime. Then, modulo \(p\),
$$Q_n-Q_i=2(n-i)(n+i+1).$$
Since \(Q_i\equiv 0 \pmod p\), we obtain
$$Q_n\equiv 0\pmod p \iff 2(n-i)(n+i+1)\equiv 0\pmod p.$$
Because \(p\) is odd, \(2\) is invertible modulo \(p\), so the only roots are
$$\boxed{n\equiv i\pmod p \quad \text{or}\quad n\equiv -i-1\pmod p.}$$
The second class is the same symmetry as
$$Q_{-n-1}=Q_n.$$
Take \(i=1\). Then \(Q_1=5\), so any term divisible by \(5\) must satisfy
$$n\equiv 1 \pmod 5 \quad \text{or}\quad n\equiv -2\equiv 3 \pmod 5.$$
Indeed,
$$Q_1=5,\qquad Q_3=25,\qquad Q_6=85,\qquad Q_8=145,$$
and all of them are divisible by \(5\). This is exactly what the sieve marks: not arbitrary indices, but two arithmetic progressions per prime factor.
The code initializes an array
$$V_n=Q_n.$$
When it reaches index \(i\), all prime factors already discovered earlier have already been divided out from every affected \(V_n\). So the current residual \(V_i\) is either \(1\) or a new prime factor that still has to be propagated.
This is why the program can safely use
$$p:=V_i$$
as the next sieving factor. It then walks through the two progressions
$$n=i+kp,\qquad n\equiv -i-1\pmod p,$$
marks those indices as composite, and repeatedly divides \(p\) out of their residual values. The repeated division is important: for example \(Q_3=25\) must lose both powers of \(5\), not just one.
If index \(i\) is still unmarked when visited, no earlier prime factor has hit it. Therefore the original number \(Q_i\) has no smaller prime divisor coming from the sequence sieve, so \(Q_i\) itself is prime and contributes to the answer. This is why the code increments the count exactly when composite[i] == 0.
From
$$4Q_n-1=(2n+1)^2$$
we get
$$-1\equiv (2n+1)^2 \pmod p$$
for every odd prime divisor \(p\mid Q_n\). So \(-1\) is a quadratic residue modulo \(p\), which implies \(p\equiv 1\pmod 4\). This explains why the prime divisors of the family are highly structured.
The implementation checks the sieve against brute force on smaller limits:
$$\text{solve}(1000)=10,\qquad \text{solve}(10^6)=175.$$
After those validations it computes the required value
$$\text{solve}(5\cdot 10^{15})=4037526.$$
The program first computes the maximum valid index \(m\), initializes values[n] = Q_n, and keeps a boolean composite array. Then it scans \(i=1,2,\dots,m\). If composite[i] is false, it counts \(Q_i\) as prime. Regardless of that flag, it looks at the current residual values[i]; if it is greater than \(1\), that residual is propagated through the two arithmetic progressions determined by \(i\) modulo that factor, and the factor is removed from every hit.
Memory usage is \(O(m)\), where \(m\approx \sqrt{\text{LIMIT}/2}\). The running time is sieve-like: each discovered factor updates two progressions of step \(p\), so the total work is vastly smaller than performing an independent primality test on every \(Q_n\). In practice this turns a prohibitively large quadratic-polynomial search into a manageable structured sieve.
Gesucht ist die Anzahl der Primwerte von
$$Q_n=2n^2+2n+1=n^2+(n+1)^2$$
unter der Schranke \(Q_n \lt \text{LIMIT}=5\cdot 10^{15}\). Ein separater Primtest für jeden Wert wäre viel zu teuer; stattdessen verwendet die Lösung ein speziell auf diese quadratische Familie zugeschnittenes Sieb.
Die Folge wächst für \(n\ge 1\) streng monoton, also genügt der Bereich
$$1\le n\le m,\qquad m=\max\{n: 2n^2+2n+1\lt \text{LIMIT}\}.$$
Aus der quadratischen Ungleichung erhält man näherungsweise
$$m\approx \left\lfloor \frac{\sqrt{2\cdot \text{LIMIT}-1}-1}{2}\right\rfloor,$$
wobei der Code anschließend mit ein paar Ganzzahltests auf die exakte strikte Schranke korrigiert.
Sei \(p=Q_i\) prim. Dann gilt modulo \(p\)
$$Q_n-Q_i=2(n-i)(n+i+1).$$
Da \(Q_i\equiv 0 \pmod p\), folgt
$$Q_n\equiv 0\pmod p \iff 2(n-i)(n+i+1)\equiv 0\pmod p.$$
Weil \(p\) ungerade ist, ist \(2\) invertierbar, also bleiben genau die beiden Wurzelklassen
$$\boxed{n\equiv i\pmod p \quad \text{oder}\quad n\equiv -i-1\pmod p.}$$
Die zweite Klasse ist dieselbe Symmetrie wie
$$Q_{-n-1}=Q_n.$$
Für \(i=1\) ist \(Q_1=5\). Also müssen alle durch \(5\) teilbaren Folgenglieder die Bedingung
$$n\equiv 1 \pmod 5 \quad \text{oder}\quad n\equiv -2\equiv 3 \pmod 5$$
erfüllen. Tatsächlich sind
$$Q_1=5,\qquad Q_3=25,\qquad Q_6=85,\qquad Q_8=145$$
alle durch \(5\) teilbar. Genau solche zwei arithmetischen Progressionen markiert das Sieb.
Der Code initialisiert
$$V_n=Q_n.$$
Wenn Index \(i\) erreicht wird, sind alle zuvor entdeckten Primfaktoren bereits aus allen betroffenen \(V_n\) herausdividiert worden. Der aktuelle Rest \(V_i\) ist also entweder \(1\) oder ein neuer Primfaktor, der noch propagiert werden muss.
Darum darf das Programm
$$p:=V_i$$
als nächsten Siebfaktor verwenden. Danach durchläuft es die beiden Progressionen
$$n=i+kp,\qquad n\equiv -i-1\pmod p,$$
markiert diese Indizes als zusammengesetzt und teilt \(p\) wiederholt aus den Restwerten heraus. Das wiederholte Teilen ist wichtig: \(Q_3=25\) muss beide Faktoren \(5\) verlieren, nicht nur einen.
Ist ein Index \(i\) beim Besuch noch unmarkiert, dann hat ihn kein früherer Primfaktor getroffen. Folglich besitzt die ursprüngliche Zahl \(Q_i\) keinen kleineren Primteiler aus der Sequenzstruktur und ist selbst prim. Genau deshalb erhöht der Code den Zähler genau dann, wenn composite[i] == 0 gilt.
Aus
$$4Q_n-1=(2n+1)^2$$
folgt für jeden ungeraden Primteiler \(p\mid Q_n\), dass
$$-1\equiv (2n+1)^2 \pmod p.$$
Also ist \(-1\) ein quadratischer Rest modulo \(p\), woraus \(p\equiv 1\pmod 4\) folgt. Das erklärt, warum die Primteiler dieser Familie stark strukturiert sind.
Die Implementierung vergleicht das Sieb zunächst mit Bruteforce bei kleineren Grenzen:
$$\text{solve}(1000)=10,\qquad \text{solve}(10^6)=175.$$
Danach wird der Zielwert berechnet:
$$\text{solve}(5\cdot 10^{15})=4037526.$$
Das Programm berechnet zunächst den maximal zulässigen Index \(m\), initialisiert values[n] = Q_n und verwaltet zusätzlich ein boolesches Feld composite. Danach läuft es \(i=1,2,\dots,m\) durch. Wenn composite[i] falsch ist, wird \(Q_i\) als Primzahl gezählt. Unabhängig davon betrachtet der Code den aktuellen Rest values[i]; ist er größer als \(1\), dann wird dieser Faktor entlang der beiden durch \(i\) modulo dieses Faktors bestimmten Progressionen propagiert und aus allen Treffern herausgeteilt.
Der Speicherbedarf ist \(O(m)\), mit \(m\approx \sqrt{\text{LIMIT}/2}\). Die Laufzeit verhält sich siebartig: Jeder entdeckte Faktor aktualisiert zwei Progressionen mit Schrittweite \(p\). Dadurch ist der Aufwand viel kleiner als ein unabhängiger Primtest für jedes einzelne \(Q_n\). Praktisch wird aus einer unbeherrschbaren Suche über ein quadratisches Polynom ein gut strukturiertes Siebproblem.
Şu dizide asal olan terimlerin sayısını bulmamız gerekiyor:
$$Q_n=2n^2+2n+1=n^2+(n+1)^2$$
ve koşul \(Q_n \lt \text{LIMIT}=5\cdot 10^{15}\). Her \(Q_n\) için ayrı asal testi yapmak çok pahalı olduğu için çözüm, bu kuadratik aileye özel bir eleme kullanır.
Dizi \(n\ge 1\) için sıkı artandır; dolayısıyla yalnızca
$$1\le n\le m,\qquad m=\max\{n: 2n^2+2n+1\lt \text{LIMIT}\}$$
aralığını taramak yeterlidir. Kuadratik eşitsizlikten
$$m\approx \left\lfloor \frac{\sqrt{2\cdot \text{LIMIT}-1}-1}{2}\right\rfloor$$
yaklaşımı gelir; kod da ardından birkaç tamsayı kontrolüyle tam sınırı ayarlar.
\(p=Q_i\) asal olsun. O zaman mod \(p\) altında
$$Q_n-Q_i=2(n-i)(n+i+1)$$
kimliği vardır. \(Q_i\equiv 0 \pmod p\) olduğundan
$$Q_n\equiv 0\pmod p \iff 2(n-i)(n+i+1)\equiv 0\pmod p.$$
\(p\) tek olduğu için \(2\) terslenebilir; bu yüzden tek kök sınıfları
$$\boxed{n\equiv i\pmod p \quad \text{veya}\quad n\equiv -i-1\pmod p}$$
olur. İkinci sınıf, aslında
$$Q_{-n-1}=Q_n$$
simetrisinin modüler yansımasıdır.
\(i=1\) için \(Q_1=5\) olur. O halde \(5\)'e bölünen tüm terimler
$$n\equiv 1 \pmod 5 \quad \text{veya}\quad n\equiv -2\equiv 3 \pmod 5$$
koşulunu sağlamalıdır. Gerçekten de
$$Q_1=5,\qquad Q_3=25,\qquad Q_6=85,\qquad Q_8=145$$
hep \(5\)'e bölünür. Eleğin işaretlediği şey tam olarak budur: rastgele indeksler değil, her asal için iki aritmetik dizi.
Kod başlangıçta
$$V_n=Q_n$$
dizisini kurar. İndeks \(i\)'ye gelindiğinde, daha önce keşfedilmiş tüm asal çarpanlar ilgili \(V_n\) değerlerinden çoktan bölünerek çıkarılmış olur. Dolayısıyla mevcut artık değer \(V_i\), ya \(1\)'dir ya da henüz yayılması gereken yeni bir asal çarpandır.
Bu yüzden program güvenle
$$p:=V_i$$
değerini bir sonraki eleme çarpanı olarak kullanır. Ardından
$$n=i+kp,\qquad n\equiv -i-1\pmod p$$
biçimindeki iki ilerleyiş boyunca gider, bu indeksleri bileşik diye işaretler ve \(p\) faktörünü ilgili artık değerlerden tekrar tekrar böler. Tekrarlı bölme önemlidir; örneğin \(Q_3=25\) için iki tane \(5\) de atılmalıdır.
Eğer indeks \(i\) ziyaret edildiğinde hâlâ işaretlenmemişse, daha küçük hiçbir asal çarpan onu vurmamış demektir. O halde özgün değer \(Q_i\), bu sıra içinde daha küçük bir asal bölen taşımıyordur ve kendisi asaldır. Kodun tam olarak composite[i] == 0 iken sayacı artırmasının nedeni budur.
$$4Q_n-1=(2n+1)^2$$
olduğundan, \(Q_n\)'i bölen her tek asal \(p\) için
$$-1\equiv (2n+1)^2 \pmod p$$
elde edilir. Yani \(-1\), mod \(p\)'de kuadratik artık olur; bu da \(p\equiv 1\pmod 4\) sonucunu verir. Bu ailedeki asal bölenlerin neden özel bir yapıya sahip olduğunu bu açıklar.
Uygulama önce eleği küçük sınırlar için brute-force ile karşılaştırır:
$$\text{solve}(1000)=10,\qquad \text{solve}(10^6)=175.$$
Ardından istenen hedef değer hesaplanır:
$$\text{solve}(5\cdot 10^{15})=4037526.$$
Program önce en büyük geçerli indeksi \(m\) hesaplar, sonra values[n] = Q_n dizisini ve buna paralel bir composite işaret dizisini oluşturur. Daha sonra \(i=1,2,\dots,m\) boyunca ilerler. Eğer composite[i] false ise \(Q_i\) asal sayılır. Bundan bağımsız olarak kod, mevcut artık değer values[i]'yi okur; bu değer \(1\)'den büyükse, o çarpanı \(i\)'nin belirlediği iki aritmetik ilerleme boyunca taşır ve denk geldiği tüm terimlerden böler.
Bellek \(O(m)\) düzeyindedir; burada \(m\approx \sqrt{\text{LIMIT}/2}\). Çalışma süresi eleğe benzer davranır: her keşfedilen çarpan, adımı \(p\) olan iki ilerlemeyi günceller. Böylece maliyet, her \(Q_n\) için bağımsız asal testi yapmaktan çok daha küçük olur. Pratikte kuadratik bir polinom ailesinde asal arama, iyi yapılandırılmış bir sieve problemine dönüşür.
Debemos contar los valores primos de
$$Q_n=2n^2+2n+1=n^2+(n+1)^2$$
bajo la cota \(Q_n \lt \text{LIMIT}=5\cdot 10^{15}\). Hacer una prueba de primalidad independiente para cada término sería demasiado lento, así que la solución usa un tamiz especializado para esta familia cuadrática.
La sucesión es estrictamente creciente para \(n\ge 1\), así que basta recorrer
$$1\le n\le m,\qquad m=\max\{n: 2n^2+2n+1\lt \text{LIMIT}\}.$$
Resolver la desigualdad cuadrática da la aproximación
$$m\approx \left\lfloor \frac{\sqrt{2\cdot \text{LIMIT}-1}-1}{2}\right\rfloor,$$
y luego el código la ajusta con comprobaciones enteras para respetar la desigualdad estricta.
Supongamos que \(p=Q_i\) es primo. Entonces, módulo \(p\),
$$Q_n-Q_i=2(n-i)(n+i+1).$$
Como \(Q_i\equiv 0 \pmod p\), obtenemos
$$Q_n\equiv 0\pmod p \iff 2(n-i)(n+i+1)\equiv 0\pmod p.$$
Dado que \(p\) es impar, \(2\) es invertible, así que las únicas raíces son
$$\boxed{n\equiv i\pmod p \quad \text{o}\quad n\equiv -i-1\pmod p.}$$
La segunda clase es la misma simetría que
$$Q_{-n-1}=Q_n.$$
Tomemos \(i=1\). Entonces \(Q_1=5\), y cualquier término divisible por \(5\) debe cumplir
$$n\equiv 1 \pmod 5 \quad \text{o}\quad n\equiv -2\equiv 3 \pmod 5.$$
En efecto,
$$Q_1=5,\qquad Q_3=25,\qquad Q_6=85,\qquad Q_8=145$$
son todos divisibles por \(5\). Eso es exactamente lo que marca el tamiz: dos progresiones aritméticas por cada factor primo, no índices arbitrarios.
El código inicializa
$$V_n=Q_n.$$
Cuando alcanza el índice \(i\), todos los factores primos descubiertos anteriormente ya han sido divididos de todos los \(V_n\) afectados. Por tanto, el residuo actual \(V_i\) es o bien \(1\), o bien un nuevo factor primo que todavía debe propagarse.
Por eso el programa puede usar con seguridad
$$p:=V_i$$
como siguiente factor de tamizado. Después recorre las dos progresiones
$$n=i+kp,\qquad n\equiv -i-1\pmod p,$$
marca esos índices como compuestos y divide repetidamente \(p\) de sus valores residuales. La división repetida es esencial: por ejemplo, \(Q_3=25\) debe perder las dos potencias de \(5\), no solo una.
Si el índice \(i\) sigue sin marcar cuando se visita, entonces ningún factor primo anterior lo ha alcanzado. En consecuencia, el valor original \(Q_i\) no tiene un divisor primo menor dentro de la estructura del tamiz y por tanto es primo. Esa es la razón por la que el código incrementa el contador exactamente cuando composite[i] == 0.
De
$$4Q_n-1=(2n+1)^2$$
se deduce que para todo primo impar \(p\mid Q_n\),
$$-1\equiv (2n+1)^2 \pmod p.$$
Es decir, \(-1\) es residuo cuadrático módulo \(p\), lo que implica \(p\equiv 1\pmod 4\). Esto explica por qué los divisores primos de esta familia tienen una estructura tan rígida.
La implementación compara primero el tamiz con fuerza bruta en límites pequeños:
$$\text{solve}(1000)=10,\qquad \text{solve}(10^6)=175.$$
Después calcula el valor pedido:
$$\text{solve}(5\cdot 10^{15})=4037526.$$
El programa calcula primero el mayor índice válido \(m\), luego inicializa values[n] = Q_n y un arreglo booleano composite. A continuación recorre \(i=1,2,\dots,m\). Si composite[i] es falso, cuenta \(Q_i\) como primo. Independientemente de eso, inspecciona el residuo actual values[i]; si es mayor que \(1\), propaga ese factor a lo largo de las dos progresiones determinadas por \(i\) módulo dicho factor y lo divide de todos los términos alcanzados.
La memoria es \(O(m)\), con \(m\approx \sqrt{\text{LIMIT}/2}\). El tiempo tiene comportamiento de tamiz: cada factor descubierto actualiza dos progresiones de paso \(p\). Eso es muchísimo mejor que hacer una prueba de primalidad separada para cada \(Q_n\). En la práctica, una búsqueda inmanejable sobre un polinomio cuadrático se convierte en un problema de tamizado estructurado.
我们要统计
$$Q_n=2n^2+2n+1=n^2+(n+1)^2$$
在 \(Q_n \lt \text{LIMIT}=5\cdot 10^{15}\) 条件下取素数的次数。若对每个 \(Q_n\) 都单独做素性测试，代价会非常高，因此代码使用了专门针对这个二次序列的筛法。
当 \(n\ge 1\) 时，该序列严格递增，所以只需遍历
$$1\le n\le m,\qquad m=\max\{n: 2n^2+2n+1\lt \text{LIMIT}\}.$$
解这个二次不等式可得近似
$$m\approx \left\lfloor \frac{\sqrt{2\cdot \text{LIMIT}-1}-1}{2}\right\rfloor,$$
代码再用整数检查微调，确保满足严格小于号。
设 \(p=Q_i\) 是素数。则在模 \(p\) 下有
$$Q_n-Q_i=2(n-i)(n+i+1).$$
由于 \(Q_i\equiv 0 \pmod p\)，于是
$$Q_n\equiv 0\pmod p \iff 2(n-i)(n+i+1)\equiv 0\pmod p.$$
因为 \(p\) 是奇素数，\(2\) 可逆，所以唯一的根类就是
$$\boxed{n\equiv i\pmod p \quad \text{或}\quad n\equiv -i-1\pmod p.}$$
第二类本质上就是
$$Q_{-n-1}=Q_n$$
这个对称性的模意义版本。
取 \(i=1\)，则 \(Q_1=5\)。因此所有能被 \(5\) 整除的项必须满足
$$n\equiv 1 \pmod 5 \quad \text{或}\quad n\equiv -2\equiv 3 \pmod 5.$$
确实有
$$Q_1=5,\qquad Q_3=25,\qquad Q_6=85,\qquad Q_8=145$$
都被 \(5\) 整除。筛法标记的正是这样的两条等差索引序列。
代码一开始设
$$V_n=Q_n.$$
当扫描到索引 \(i\) 时，所有更早发现的素因子都已经从相应的 \(V_n\) 中除掉了。因此当前残值 \(V_i\) 不是 \(1\)，就是一个尚未传播的新素因子。
所以程序可以安全地把
$$p:=V_i$$
当作下一步筛因子。然后沿着两条进程
$$n=i+kp,\qquad n\equiv -i-1\pmod p$$
前进，把这些索引标成合数，并且把残值里的所有 \(p\) 因子全部除去。反复整除是必要的，例如 \(Q_3=25\) 必须把两个 \(5\) 都除掉。
如果访问索引 \(i\) 时它仍未被标记，说明没有任何更早的素因子命中过它。因此原始值 \(Q_i\) 没有更小的筛内素因子，于是它本身就是素数。这正是代码在 composite[i] == 0 时增加计数的原因。
由
$$4Q_n-1=(2n+1)^2$$
可知，对任何奇素因子 \(p\mid Q_n\)，都有
$$-1\equiv (2n+1)^2 \pmod p.$$
也就是说 \(-1\) 在模 \(p\) 下是二次剩余，因此 \(p\equiv 1\pmod 4\)。这解释了该序列素因子的特殊结构。
实现先用暴力核对小范围：
$$\text{solve}(1000)=10,\qquad \text{solve}(10^6)=175.$$
然后计算目标值
$$\text{solve}(5\cdot 10^{15})=4037526.$$
程序先求出最大合法索引 \(m\)，再初始化 values[n] = Q_n 和布尔数组 composite。随后按 \(i=1,2,\dots,m\) 递增扫描。如果 composite[i] 仍为假，就把 \(Q_i\) 计入素数。无论是否计数，程序都会查看当前残值 values[i]；若它大于 \(1\)，就把这个因子沿着由 \(i\) 模该因子决定的两条等差序列传播，并从命中的残值中不断除掉它。
空间复杂度为 \(O(m)\)，其中 \(m\approx \sqrt{\text{LIMIT}/2}\)。时间复杂度呈筛法行为：每个发现的因子只更新两条步长为 \(p\) 的序列。这比对每个 \(Q_n\) 单独做素性测试快得多。于是原本难以承受的二次多项式搜索，被转化成了一个结构良好的专用筛问题。
Нужно посчитать простые значения выражения
$$Q_n=2n^2+2n+1=n^2+(n+1)^2$$
при условии \(Q_n \lt \text{LIMIT}=5\cdot 10^{15}\). Независимая проверка простоты для каждого \(Q_n\) слишком дорога, поэтому решение использует решето, специально адаптированное к этой квадратичной семье.
Последовательность строго возрастает при \(n\ge 1\), поэтому достаточно рассматривать
$$1\le n\le m,\qquad m=\max\{n: 2n^2+2n+1\lt \text{LIMIT}\}.$$
Решение квадратного неравенства даёт приближение
$$m\approx \left\lfloor \frac{\sqrt{2\cdot \text{LIMIT}-1}-1}{2}\right\rfloor,$$
после чего код корректирует результат несколькими целочисленными проверками, чтобы сохранить строгое неравенство.
Пусть \(p=Q_i\) — простой. Тогда по модулю \(p\)
$$Q_n-Q_i=2(n-i)(n+i+1).$$
Так как \(Q_i\equiv 0 \pmod p\), получаем
$$Q_n\equiv 0\pmod p \iff 2(n-i)(n+i+1)\equiv 0\pmod p.$$
Поскольку \(p\) нечётен, число \(2\) обратимо, и остаются ровно два класса корней:
$$\boxed{n\equiv i\pmod p \quad \text{или}\quad n\equiv -i-1\pmod p.}$$
Второй класс — это та же симметрия, что и
$$Q_{-n-1}=Q_n.$$
Возьмём \(i=1\). Тогда \(Q_1=5\), и все члены, делящиеся на \(5\), должны удовлетворять
$$n\equiv 1 \pmod 5 \quad \text{или}\quad n\equiv -2\equiv 3 \pmod 5.$$
Действительно,
$$Q_1=5,\qquad Q_3=25,\qquad Q_6=85,\qquad Q_8=145$$
все делятся на \(5\). Именно такие две арифметические прогрессии и отмечает решето.
Код инициализирует массив
$$V_n=Q_n.$$
К моменту, когда индекс \(i\) достигается, все ранее найденные простые множители уже удалены из всех затронутых \(V_n\). Поэтому текущий остаток \(V_i\) либо равен \(1\), либо является новым простым множителем, который ещё нужно распространить.
Именно поэтому программа может безопасно взять
$$p:=V_i$$
как следующий фактор для решета. После этого она проходит по двум прогрессиям
$$n=i+kp,\qquad n\equiv -i-1\pmod p,$$
помечает эти индексы как составные и многократно делит \(p\) из остаточных значений. Многократное деление существенно: например, из \(Q_3=25\) надо убрать обе степени числа \(5\).
Если к моменту посещения индекс \(i\) всё ещё непомечен, значит ни один меньший простой множитель его не задел. Следовательно, исходное число \(Q_i\) не имеет меньшего простого делителя в структуре решета и само является простым. Именно поэтому код увеличивает счётчик тогда и только тогда, когда composite[i] == 0.
Из равенства
$$4Q_n-1=(2n+1)^2$$
следует, что для любого нечётного простого делителя \(p\mid Q_n\)
$$-1\equiv (2n+1)^2 \pmod p.$$
То есть \(-1\) является квадратичным вычетом по модулю \(p\), откуда следует \(p\equiv 1\pmod 4\). Это объясняет жёсткую структуру простых делителей этой семьи.
Сначала реализация сверяет решето с брутфорсом на меньших пределах:
$$\text{solve}(1000)=10,\qquad \text{solve}(10^6)=175.$$
После этого вычисляется искомое значение
$$\text{solve}(5\cdot 10^{15})=4037526.$$
Программа сначала находит максимальный допустимый индекс \(m\), затем инициализирует values[n] = Q_n и булев массив composite. Далее она проходит по \(i=1,2,\dots,m\). Если composite[i] ложно, то \(Q_i\) засчитывается как простое. Независимо от этого код смотрит на текущий остаток values[i]; если он больше \(1\), то этот множитель распространяется вдоль двух арифметических прогрессий, задаваемых индексом \(i\) по модулю этого множителя, и удаляется из всех попавших значений.
Память равна \(O(m)\), где \(m\approx \sqrt{\text{LIMIT}/2}\). Время работы ведёт себя как у решета: каждый найденный множитель обновляет две прогрессии с шагом \(p\). Это намного дешевле, чем независимо тестировать простоту каждого \(Q_n\). На практике громоздкий поиск по квадратичному полиному превращается в хорошо структурированное решето.
نريد عدّ القيم الأولية للتعبير
$$Q_n=2n^2+2n+1=n^2+(n+1)^2$$
تحت الشرط \(Q_n \lt \text{LIMIT}=5\cdot 10^{15}\). واختبار أولية كل قيمة \(Q_n\) على حدة سيكون مكلفًا جدًا، لذلك يستعمل الحل غربالًا مخصصًا لهذه العائلة التربيعية.
بما أن المتتالية تزداد بصرامة عندما \(n\ge 1\)، يكفي أن نبحث في المجال
$$1\le n\le m,\qquad m=\max\{n: 2n^2+2n+1\lt \text{LIMIT}\}.$$
ومن حل المتباينة التربيعية نحصل تقريبًا على
$$m\approx \left\lfloor \frac{\sqrt{2\cdot \text{LIMIT}-1}-1}{2}\right\rfloor,$$
ثم يقوم الكود بتعديل هذه القيمة بفحوص صحيحة صغيرة حتى يضمن شرط الأصغرية الصارم.
إذا كان \(p=Q_i\) عددًا أوليًا، فلدينا بترديد \(p\)
$$Q_n-Q_i=2(n-i)(n+i+1).$$
وبما أن \(Q_i\equiv 0 \pmod p\)، فإن
$$Q_n\equiv 0\pmod p \iff 2(n-i)(n+i+1)\equiv 0\pmod p.$$
ولأن \(p\) فردي، فإن \(2\) قابل للعكس، ومن ثم تكون فئتا الجذور الوحيدتان هما
$$\boxed{n\equiv i\pmod p \quad \text{أو}\quad n\equiv -i-1\pmod p.}$$
والفئة الثانية ليست إلا الصياغة المعيارية للتناظر
$$Q_{-n-1}=Q_n.$$
خذ \(i=1\)، فنحصل على \(Q_1=5\). عندئذٍ كل حد قابل للقسمة على \(5\) يجب أن يحقق
$$n\equiv 1 \pmod 5 \quad \text{أو}\quad n\equiv -2\equiv 3 \pmod 5.$$
وفعلًا نجد أن
$$Q_1=5,\qquad Q_3=25,\qquad Q_6=85,\qquad Q_8=145$$
كلها قابلة للقسمة على \(5\). هذا هو بالضبط ما يعلّمه الغربال: متتاليتان حسابيتان لكل عامل أولي.
يبدأ الكود بتعريف
$$V_n=Q_n.$$
وعندما يصل إلى الفهرس \(i\)، تكون جميع العوامل الأولية التي اكتُشفت سابقًا قد أزيلت بالفعل من كل القيم \(V_n\) المتأثرة. لذلك فإن الباقي الحالي \(V_i\) إما أن يكون \(1\)، أو عاملًا أوليًا جديدًا ما زال يحتاج إلى نشر.
ولهذا يستطيع البرنامج أن يستعمل بأمان
$$p:=V_i$$
بوصفه عامل الغربلة التالي. ثم يسير على التقدّمين
$$n=i+kp,\qquad n\equiv -i-1\pmod p,$$
ويعلّم هذه الفهارس كمركبة، ويقسم العامل \(p\) مرارًا من القيم المتبقية. والتقسيم المتكرر مهم جدًا؛ فمثلًا \(Q_3=25\) يجب أن يفقد قوتي \(5\)، لا واحدة فقط.
إذا بقي الفهرس \(i\) غير معلَّم عندما نزوره، فهذا يعني أن أي عامل أولي أصغر لم يضربه سابقًا. وبالتالي فإن القيمة الأصلية \(Q_i\) لا تملك قاسمًا أوليًا أصغر ضمن بنية الغربال، ومن ثم فهي أولية بنفسها. ولهذا يزيد الكود العداد تمامًا عندما تكون composite[i] == 0.
من العلاقة
$$4Q_n-1=(2n+1)^2$$
نستنتج أنه لكل قاسم أولي فردي \(p\mid Q_n\) لدينا
$$-1\equiv (2n+1)^2 \pmod p.$$
أي إن \(-1\) باقٍ تربيعي modulo \(p\)، ومنه يلزم \(p\equiv 1\pmod 4\). وهذا يفسر البنية المقيدة لعوامل هذه العائلة.
يقارن التنفيذ الغربال أولًا مع brute force عند حدود أصغر:
$$\text{solve}(1000)=10,\qquad \text{solve}(10^6)=175.$$
ثم يحسب القيمة المطلوبة:
$$\text{solve}(5\cdot 10^{15})=4037526.$$
يحسب البرنامج أولًا أكبر فهرس صالح \(m\)، ثم يهيّئ المصفوفة values[n] = Q_n مع مصفوفة منطقية اسمها composite. بعد ذلك يمسح \(i=1,2,\dots,m\). فإذا كانت composite[i] غير مفعلة، يُحسب \(Q_i\) على أنه أولي. وبغض النظر عن ذلك، ينظر الكود إلى الباقي الحالي values[i]؛ فإذا كان أكبر من \(1\)، ينشر هذا العامل على طول المتتاليتين الحسابيتين اللتين تحددهما بقايا \(i\) modulo ذلك العامل، ثم يزيله من كل القيم التي يصيبها.
التعقيد في الذاكرة هو \(O(m)\)، حيث \(m\approx \sqrt{\text{LIMIT}/2}\). أما الزمن فيتصرف كسلوك الغربال: كل عامل مكتشف يحدّث متتاليتين بخطوة \(p\). وهذا أسرع بكثير من إجراء اختبار أولية مستقل لكل \(Q_n\). عمليًا تتحول عملية بحث صعبة جدًا داخل كثير حدود تربيعي إلى غربلة منظمة يمكن تنفيذها.