We must evaluate the prime sum
$$D(a,b,k)=\sum_{\substack{a \le p \lt a+b\\ p\ \text{prime}}} d(p,p-1,k).$$
The key fact used by the implementation is that for every contributing prime \(p\), the required term collapses to a modular inverse of \(k-1\) modulo \(p\). After that simplification, the whole problem becomes: enumerate the primes in a long interval efficiently, compute one inverse for each of them, and add the results.
The number-theoretic core is simple once the prime-term identity is recognized. The difficult part is not symbolic algebra, but organizing the computation so that the interval \([a,a+b)\) can be processed without sieving every number up to \(a+b\).
The implementation uses the identity
$$d(p,p-1,k)\equiv (k-1)^{-1}\pmod{p}.$$
So for each prime \(p\) in the interval, the required contribution is the unique residue \(x\) satisfying
$$(k-1)x \equiv 1 \pmod{p}.$$
Because the modulus is prime, this inverse exists whenever \(p \nmid (k-1)\). In other words, the original arithmetic object \(d(p,p-1,k)\) does not need to be built directly; it is enough to solve one linear congruence per prime.
When \(\gcd(k-1,p)=1\), the extended Euclidean algorithm finds integers \(u\) and \(v\) such that
$$u(k-1)+vp=1.$$
Reducing that identity modulo \(p\) removes the second term and leaves
$$u(k-1)\equiv 1 \pmod{p}.$$
Therefore the required contribution is simply the least nonnegative residue of \(u\):
$$d(p,p-1,k)=u \bmod p.$$
This is why the three implementations use the extended Euclidean algorithm rather than repeated powering. It produces the inverse directly and works uniformly in all three languages.
The interval can start near \(10^9\), so ordinary sieving from \(2\) all the way to \(a+b\) would be wasteful. Instead, let
$$R=a+b-1,\qquad s=\left\lfloor \sqrt{R}\right\rfloor.$$
First compute all primes up to \(s\) with a standard sieve. These are the only base primes needed to mark composites in the target interval. For each base prime \(q\), the first multiple that must be removed from \([a,a+b)\) is
$$m_q=\max\left(q^2,\left\lceil\frac{a}{q}\right\rceil q\right).$$
Then every number
$$m_q,\ m_q+q,\ m_q+2q,\dots$$
inside the interval is composite and can be marked. After all base primes are processed, the unmarked entries are exactly the primes in \([a,a+b)\).
Once the interval primes are known, the sum is evaluated term by term:
$$D(a,b,k)=\sum_{\substack{a \le p \lt a+b\\ p\ \text{prime}}} \left((k-1)^{-1} \bmod p\right).$$
No further combinatorial structure is hidden here. The algorithm does precisely two things: identify the primes in the interval and attach to each prime the modular inverse of \(k-1\).
The interval \([101,102)\) contains only one prime, namely \(p=101\). Here
$$k-1=9.$$
So we solve
$$9x \equiv 1 \pmod{101}.$$
Since
$$9\cdot 45 = 405 = 4\cdot 101 + 1,$$
the inverse of \(9\) modulo \(101\) is \(45\). Therefore
$$D(101,1,10)=45,$$
which matches the checkpoint used by the implementation.
The C++, Python, and Java implementations follow the same pipeline. They first sieve all primes up to \(\lfloor\sqrt{a+b-1}\rfloor\). Next they allocate a Boolean segment of length \(b\) representing the numbers \(a,a+1,\dots,a+b-1\), and every composite hit by a base prime is marked. After that scan, each unmarked position corresponds to a prime \(p\) in the interval.
For every such prime, the implementation runs the extended Euclidean algorithm on \(k-1\) and \(p\), normalizes the Bézout coefficient into the range \(0\) to \(p-1\), and adds that value to the running total. The inverse routine is written defensively: if the reduced value of \(k-1\) is \(0\) modulo \(p\), or if the gcd is not \(1\), it returns \(0\). For the intended evaluations, the contributing primes are coprime to \(k-1\), so each genuine term is well defined.
Let \(R=a+b-1\). Generating all base primes up to \(\lfloor\sqrt{R}\rfloor\) costs \(O(\sqrt{R}\log\log R)\) time and \(O(\sqrt{R})\) memory with the simple sieve used in the implementations. Marking the interval of length \(b\) costs \(O(b\log\log R)\) time and \(O(b)\) memory. If \(\pi([a,a+b))\) denotes the number of primes in the interval, the inverse computations add \(O(\pi([a,a+b))\log R)\) time. So the total memory usage is linear in the interval length, and the method is efficient because it never stores all integers up to \(a+b\).
Gesucht ist die Primzahlsumme
$$D(a,b,k)=\sum_{\substack{a \le p \lt a+b\\ p\ \text{prim}}} d(p,p-1,k).$$
Die entscheidende Beobachtung der Implementierung lautet: Für jede beitragende Primzahl \(p\) reduziert sich der gesuchte Summand auf das modulare Inverse von \(k-1\) modulo \(p\). Danach besteht das Problem nur noch aus zwei Teilen: Primzahlen in einem langen Intervall effizient erzeugen und für jede dieser Primzahlen genau ein Inverses berechnen.
Sobald diese Vereinfachung erkannt ist, bleibt im Wesentlichen ein sauber organisiertes Zahlentheorie-Problem übrig. Nicht die Formel selbst ist teuer, sondern die effiziente Behandlung des Intervalls \([a,a+b)\).
Die Implementierung verwendet die Identität
$$d(p,p-1,k)\equiv (k-1)^{-1}\pmod{p}.$$
Für jede Primzahl \(p\) im Intervall ist also genau der Rest \(x\) gesucht, der
$$(k-1)x \equiv 1 \pmod{p}$$
erfüllt. Da der Modul prim ist, existiert dieses Inverse genau dann, wenn \(p\) kein Teiler von \(k-1\) ist. Das bedeutet: Man muss das arithmetische Objekt \(d(p,p-1,k)\) nicht direkt konstruieren, sondern nur eine lineare Kongruenz pro Primzahl lösen.
Gilt \(\gcd(k-1,p)=1\), so liefert der erweiterte euklidische Algorithmus ganze Zahlen \(u\) und \(v\) mit
$$u(k-1)+vp=1.$$
Reduziert man diese Gleichung modulo \(p\), verschwindet der zweite Term und es bleibt
$$u(k-1)\equiv 1 \pmod{p}.$$
Damit ist der gesuchte Summand einfach der kleinste nichtnegative Rest von \(u\):
$$d(p,p-1,k)=u \bmod p.$$
Deshalb benutzen die drei Implementierungen den erweiterten euklidischen Algorithmus statt modularer Potenzierung. Das Inverse wird direkt erzeugt und die Methode funktioniert in allen drei Sprachen identisch.
Das Intervall beginnt unter Umständen in der Nähe von \(10^9\). Ein vollständiges Sieb bis \(a+b\) wäre daher unnötig teuer. Setze
$$R=a+b-1,\qquad s=\left\lfloor \sqrt{R}\right\rfloor.$$
Zuerst werden alle Primzahlen bis \(s\) mit einem gewöhnlichen Sieb berechnet. Mehr Basisprimzahlen braucht man nicht, um im Zielintervall alle zusammengesetzten Zahlen zu markieren. Für jede Basisprimzahl \(q\) beginnt das Streichen bei
$$m_q=\max\left(q^2,\left\lceil\frac{a}{q}\right\rceil q\right).$$
Dann werden
$$m_q,\ m_q+q,\ m_q+2q,\dots$$
im Intervall markiert. Nach Verarbeitung aller Basisprimzahlen sind genau die unmarkierten Zahlen die Primzahlen in \([a,a+b)\).
Sobald die Intervallprimzahlen bekannt sind, wird die Summe direkt ausgewertet:
$$D(a,b,k)=\sum_{\substack{a \le p \lt a+b\\ p\ \text{prim}}} \left((k-1)^{-1} \bmod p\right).$$
Es steckt keine weitere kombinatorische Schicht dahinter. Der Algorithmus identifiziert nur die Primzahlen im Intervall und ordnet jeder davon das modulare Inverse von \(k-1\) zu.
Das Intervall \([101,102)\) enthält nur die Primzahl \(101\). Hier ist
$$k-1=9.$$
Also lösen wir
$$9x \equiv 1 \pmod{101}.$$
Wegen
$$9\cdot 45 = 405 = 4\cdot 101 + 1$$
ist das Inverse von \(9\) modulo \(101\) gleich \(45\). Daher
$$D(101,1,10)=45,$$
genau der Kontrollwert der Implementierung.
Die C++-, Python- und Java-Implementierungen folgen derselben Pipeline. Zuerst werden alle Primzahlen bis \(\lfloor\sqrt{a+b-1}\rfloor\) gesiebt. Danach wird ein Boolesches Segment der Länge \(b\) angelegt, das die Zahlen \(a,a+1,\dots,a+b-1\) repräsentiert; jede von einer Basisprimzahl getroffene zusammengesetzte Zahl wird markiert. Nach diesem Schritt entspricht jede unmarkierte Position einer Primzahl \(p\) im Intervall.
Für jede dieser Primzahlen wird der erweiterte euklidische Algorithmus auf \(k-1\) und \(p\) angewandt, der Bézout-Koeffizient in den Bereich \(0\) bis \(p-1\) normalisiert und zum laufenden Ergebnis addiert. Die Inversenroutine ist defensiv formuliert: Falls der reduzierte Wert von \(k-1\) gleich \(0\) modulo \(p\) ist oder der ggT nicht \(1\) ist, wird \(0\) zurückgegeben. Bei den vorgesehenen Auswertungen sind die beitragenden Primzahlen jedoch teilerfremd zu \(k-1\), daher ist jeder echte Summand wohldefiniert.
Mit \(R=a+b-1\) kostet das Erzeugen aller Basisprimzahlen bis \(\lfloor\sqrt{R}\rfloor\) im einfachen Sieb \(O(\sqrt{R}\log\log R)\) Zeit und \(O(\sqrt{R})\) Speicher. Das Markieren des Intervalls der Länge \(b\) benötigt \(O(b\log\log R)\) Zeit und \(O(b)\) Speicher. Bezeichnet \(\pi([a,a+b))\) die Anzahl der Primzahlen im Intervall, so tragen die Inversenberechnungen zusätzlich \(O(\pi([a,a+b))\log R)\) Zeit bei. Der Speicherverbrauch ist also linear in der Intervalllänge, und das Verfahren bleibt effizient, weil niemals alle Zahlen bis \(a+b\) gehalten werden müssen.
Hesaplanması gereken ifade
$$D(a,b,k)=\sum_{\substack{a \le p \lt a+b\\ p\ \text{asal}}} d(p,p-1,k).$$
Uygulamanın kullandığı temel gözlem şudur: katkı veren her asal \(p\) için ilgili terim, \(k-1\)'in \(p\)'ye göre modüler tersi olur. Bu sadeleşmeden sonra problem iki parçaya ayrılır: uzun bir aralıktaki asalları verimli biçimde bulmak ve her asal için bir tane modüler ters hesaplayıp toplamak.
Asıl sayı teorisi kısmı, bu indirgeme fark edildiğinde oldukça nettir. Zor olan şey cebir değil, \([a,a+b)\) aralığını gereksiz iş yapmadan işlemektir.
Uygulama şu özdeşliği kullanır:
$$d(p,p-1,k)\equiv (k-1)^{-1}\pmod{p}.$$
Dolayısıyla aralıktaki her asal \(p\) için bulunması gereken katkı,
$$(k-1)x \equiv 1 \pmod{p}$$
denklemini sağlayan tek kalandır. Modül asal olduğu için bu ters, ancak ve ancak \(p \nmid (k-1)\) ise vardır. Başka bir deyişle \(d(p,p-1,k)\) değerini doğrudan kurmaya gerek yoktur; her asal için bir doğrusal kongruens çözmek yeterlidir.
\(\gcd(k-1,p)=1\) olduğunda genişletilmiş Öklid algoritması \(u\) ve \(v\) tamsayılarını
$$u(k-1)+vp=1$$
eşitliğini sağlayacak şekilde verir. Bu bağıntı \(p\)'ye göre mod alınınca ikinci terim kaybolur ve
$$u(k-1)\equiv 1 \pmod{p}$$
kalır. O halde aranan katkı, \(u\)'nun en küçük negatif olmayan kalıntısıdır:
$$d(p,p-1,k)=u \bmod p.$$
Bu nedenle C++, Python ve Java uygulamaları tekrarlı üs alma yerine genişletilmiş Öklid kullanır. Modüler ters doğrudan elde edilir ve yöntem her dilde aynıdır.
Aralık \(10^9\) civarında başlayabildiği için, \(2\)'den \(a+b\)'ye kadar klasik eleme yapmak gereksizdir. Şunu tanımlayalım:
$$R=a+b-1,\qquad s=\left\lfloor \sqrt{R}\right\rfloor.$$
Önce \(s\)'ye kadar olan asallar klasik bir elekle bulunur. Hedef aralıktaki bileşikleri işaretlemek için bundan fazlasına gerek yoktur. Her taban asal \(q\) için, \([a,a+b)\) içinde elenmesi gereken ilk kat
$$m_q=\max\left(q^2,\left\lceil\frac{a}{q}\right\rceil q\right)$$
olur. Sonra
$$m_q,\ m_q+q,\ m_q+2q,\dots$$
değerleri aralık içinde işaretlenir. Tüm taban asallar işlendiğinde işaretsiz kalan sayılar tam olarak \([a,a+b)\) aralığındaki asallardır.
Aralıktaki asallar belirlendikten sonra toplam doğrudan hesaplanır:
$$D(a,b,k)=\sum_{\substack{a \le p \lt a+b\\ p\ \text{asal}}} \left((k-1)^{-1} \bmod p\right).$$
Burada gizli başka bir kombinatorik yapı yoktur. Algoritma yalnızca aralıktaki asalları tespit eder ve her asala \(k-1\)'in modüler tersini eşler.
\([101,102)\) aralığında yalnızca \(101\) asalı vardır. Bu durumda
$$k-1=9.$$
Dolayısıyla
$$9x \equiv 1 \pmod{101}$$
denklemini çözeriz. Şöyle ki
$$9\cdot 45 = 405 = 4\cdot 101 + 1,$$
yani \(9\)'un \(101\) modundaki tersi \(45\)'tir. O halde
$$D(101,1,10)=45,$$
ve bu değer uygulamadaki kontrol sonucu ile aynıdır.
C++, Python ve Java uygulamaları aynı işlem sırasını izler. Önce \(\lfloor\sqrt{a+b-1}\rfloor\)'ye kadar olan asallar elenir. Sonra \(a,a+1,\dots,a+b-1\) sayıları için uzunluğu \(b\) olan Boole dizisi ayrılır ve her taban asalın vurduğu bileşikler işaretlenir. Bu tarama tamamlanınca işaretlenmemiş her konum aralıktaki bir asal \(p\)'yi temsil eder.
Her böyle asal için uygulama, \(k-1\) ile \(p\) üzerinde genişletilmiş Öklid algoritmasını çalıştırır, elde edilen Bézout katsayısını \(0\) ile \(p-1\) arasına normalize eder ve toplama ekler. Ters alma yordamı savunmacı yazılmıştır: \(k-1\)'in indirgenmiş değeri \(p\) modunda \(0\) ise veya ortak bölen \(1\) değilse \(0\) döner. Hedef değerlendirmelerde katkı veren asallar \(k-1\) ile aralarında asaldır, dolayısıyla gerçek her terim tanımlıdır.
\(R=a+b-1\) olsun. \(\lfloor\sqrt{R}\rfloor\)'ye kadar taban asalları üretmek, kullanılan basit elekte \(O(\sqrt{R}\log\log R)\) zaman ve \(O(\sqrt{R})\) bellek gerektirir. Uzunluğu \(b\) olan hedef aralığı işaretlemek \(O(b\log\log R)\) zaman ve \(O(b)\) bellek kullanır. \([a,a+b)\) aralığındaki asal sayısını \(\pi([a,a+b))\) ile gösterirsek, modüler ters hesapları buna ek olarak \(O(\pi([a,a+b))\log R)\) zaman getirir. Dolayısıyla baskın bellek kullanımı aralık uzunluğunda lineerdir ve yöntem, \(a+b\)'ye kadar tüm sayıları tutmadığı için verimlidir.
Debemos evaluar la suma sobre primos
$$D(a,b,k)=\sum_{\substack{a \le p \lt a+b\\ p\ \text{primo}}} d(p,p-1,k).$$
La observación clave usada por la implementación es que, para cada primo \(p\) que contribuye, el término pedido se reduce al inverso modular de \(k-1\) módulo \(p\). Después de esa simplificación, el problema se convierte en dos tareas: enumerar eficientemente los primos de un intervalo largo y calcular un solo inverso para cada uno de ellos.
Una vez reconocida esta reducción, la parte aritmética es directa. La verdadera dificultad está en procesar el intervalo \([a,a+b)\) sin construir una criba completa hasta \(a+b\).
La implementación utiliza la identidad
$$d(p,p-1,k)\equiv (k-1)^{-1}\pmod{p}.$$
Así, para cada primo \(p\) del intervalo, la contribución buscada es el único residuo \(x\) que satisface
$$(k-1)x \equiv 1 \pmod{p}.$$
Como el módulo es primo, este inverso existe siempre que \(p \nmid (k-1)\). En otras palabras, no hace falta construir directamente el objeto aritmético \(d(p,p-1,k)\); basta con resolver una congruencia lineal para cada primo.
Si \(\gcd(k-1,p)=1\), el algoritmo extendido de Euclides produce enteros \(u\) y \(v\) tales que
$$u(k-1)+vp=1.$$
Al reducir esta identidad módulo \(p\), el segundo término desaparece y queda
$$u(k-1)\equiv 1 \pmod{p}.$$
Por lo tanto, la contribución pedida es simplemente el residuo no negativo de \(u\):
$$d(p,p-1,k)=u \bmod p.$$
Por esa razón, las implementaciones en C++, Python y Java usan Euclides extendido en lugar de exponenciación modular repetida. El inverso sale de forma directa y el procedimiento es uniforme.
El intervalo puede empezar cerca de \(10^9\), de modo que cribar desde \(2\) hasta \(a+b\) sería innecesario. Definamos
$$R=a+b-1,\qquad s=\left\lfloor \sqrt{R}\right\rfloor.$$
Primero se calculan todos los primos hasta \(s\) con una criba ordinaria. Esos son exactamente los primos base necesarios para marcar compuestos dentro del intervalo objetivo. Para cada primo base \(q\), el primer múltiplo que debe tacharse en \([a,a+b)\) es
$$m_q=\max\left(q^2,\left\lceil\frac{a}{q}\right\rceil q\right).$$
Luego se marcan
$$m_q,\ m_q+q,\ m_q+2q,\dots$$
mientras sigan dentro del intervalo. Al terminar con todos los primos base, las posiciones no marcadas son exactamente los primos de \([a,a+b)\).
Una vez conocidos los primos del intervalo, la suma final es
$$D(a,b,k)=\sum_{\substack{a \le p \lt a+b\\ p\ \text{primo}}} \left((k-1)^{-1} \bmod p\right).$$
No hay una capa combinatoria adicional. El algoritmo hace exactamente dos cosas: localizar los primos del intervalo y asociar a cada uno el inverso modular de \(k-1\).
El intervalo \([101,102)\) contiene un solo primo, \(p=101\). Aquí
$$k-1=9.$$
Debemos resolver
$$9x \equiv 1 \pmod{101}.$$
Como
$$9\cdot 45 = 405 = 4\cdot 101 + 1,$$
el inverso de \(9\) módulo \(101\) es \(45\). Por tanto,
$$D(101,1,10)=45,$$
que coincide con el valor de comprobación usado por la implementación.
Las implementaciones en C++, Python y Java siguen la misma secuencia. Primero generan todos los primos hasta \(\lfloor\sqrt{a+b-1}\rfloor\). Después reservan un segmento booleano de longitud \(b\) que representa los enteros \(a,a+1,\dots,a+b-1\), y marcan cada compuesto alcanzado por un primo base. Cuando termina ese barrido, cada posición no marcada corresponde a un primo \(p\) del intervalo.
Para cada uno de esos primos, la implementación ejecuta el algoritmo extendido de Euclides sobre \(k-1\) y \(p\), normaliza el coeficiente de Bézout al rango \(0\) a \(p-1\) y lo suma al acumulador. La rutina del inverso está escrita de forma defensiva: si el valor reducido de \(k-1\) es \(0\) módulo \(p\), o si el mcd no es \(1\), devuelve \(0\). En las evaluaciones previstas, los primos que contribuyen son coprimos con \(k-1\), así que cada término real está bien definido.
Sea \(R=a+b-1\). Generar todos los primos base hasta \(\lfloor\sqrt{R}\rfloor\) cuesta \(O(\sqrt{R}\log\log R)\) tiempo y \(O(\sqrt{R})\) memoria con la criba simple usada por las implementaciones. Marcar el intervalo de longitud \(b\) requiere \(O(b\log\log R)\) tiempo y \(O(b)\) memoria. Si \(\pi([a,a+b))\) es el número de primos del intervalo, los cálculos de inversos añaden \(O(\pi([a,a+b))\log R)\) tiempo. En consecuencia, la memoria dominante es lineal en la longitud del intervalo, y el método es eficiente porque no necesita almacenar todos los enteros hasta \(a+b\).
目标是计算下面这段素数求和：
$$D(a,b,k)=\sum_{\substack{a \le p \lt a+b\\ p\ \text{为素数}}} d(p,p-1,k).$$
实现所依赖的关键事实是：对每一个真正参与求和的素数 \(p\)，所需的项都可以化成 \(k-1\) 在模 \(p\) 下的乘法逆元。这样一来，原题就被拆成两个非常明确的任务：第一，在很大的区间 \([a,a+b)\) 中高效找出所有素数；第二，对每个区间素数计算一次模逆并累加。
只要接受这个核心恒等式，后面的推导就相当直接。真正需要仔细设计的是计算流程，而不是复杂的符号变换。因为区间可能接近 \(10^9\)，所以算法必须避免从 \(2\) 一直筛到 \(a+b\)。
实现使用的核心等式是
$$d(p,p-1,k)\equiv (k-1)^{-1}\pmod{p}.$$
因此，对区间中的每个素数 \(p\)，真正要找的是唯一的剩余类 \(x\)，满足
$$(k-1)x \equiv 1 \pmod{p}.$$
由于模数 \(p\) 是素数，只要 \(p\) 不整除 \(k-1\)，这个逆元就存在且唯一。也就是说，计算时根本不需要直接构造更复杂的 \(d(p,p-1,k)\) 形式，只需要对每个素数求解一个一次同余方程。
当 \(\gcd(k-1,p)=1\) 时，扩展欧几里得算法会给出整数 \(u\) 与 \(v\)，使得
$$u(k-1)+vp=1.$$
把这条等式对 \(p\) 取模，第二项会消失，于是得到
$$u(k-1)\equiv 1 \pmod{p}.$$
所以要求的贡献值就是 \(u\) 在模 \(p\) 下的最小非负代表元：
$$d(p,p-1,k)=u \bmod p.$$
这正是三个实现都采用扩展欧几里得算法，而不是重复平方快速幂的原因。快速幂当然也能在素数模下求逆，但扩展欧几里得会直接返回满足 Bézout 等式的系数，因此路线更直接，而且对 C++、Python、Java 三种写法都很统一。
设
$$R=a+b-1,\qquad s=\left\lfloor \sqrt{R}\right\rfloor.$$
先用普通筛法求出所有不超过 \(s\) 的素数。这些素数就是后续分段筛中需要的全部“基准素数”。原因很简单：任何合数如果落在 \([a,a+b)\) 中，它一定有一个不超过其平方根的素因子，而该平方根又不超过 \(\sqrt{R}\)。
接下来，对每个基准素数 \(q\)，都要在区间 \([a,a+b)\) 内找到第一个需要标记的倍数。这个起点是
$$m_q=\max\left(q^2,\left\lceil\frac{a}{q}\right\rceil q\right).$$
这里取最大值有两个目的。第一，\(\left\lceil\frac{a}{q}\right\rceil q\) 保证我们从区间内部的第一个 \(q\) 的倍数开始；第二，\(q^2\) 保证不会重复标记那些早就应该由更小素数处理的合数。于是区间中所有
$$m_q,\ m_q+q,\ m_q+2q,\dots$$
都可以标记为合数。把所有基准素数都处理完之后，仍然没有被标记的位置就恰好对应区间中的素数。
一旦区间内的素数已经枚举完毕，最终结果就直接写成
$$D(a,b,k)=\sum_{\substack{a \le p \lt a+b\\ p\ \text{为素数}}} \left((k-1)^{-1} \bmod p\right).$$
这里没有隐藏的二次转换，也没有额外的组合计数。整个算法做的事情非常纯粹：先筛出区间素数，再把 \(k-1\) 在每个素数模数下的逆元累加起来。
区间 \([101,102)\) 中只有一个素数，就是 \(p=101\)。此时
$$k-1=9.$$
于是我们要求解
$$9x \equiv 1 \pmod{101}.$$
因为
$$9\cdot 45 = 405 = 4\cdot 101 + 1,$$
所以 \(9\) 在模 \(101\) 下的逆元就是 \(45\)。区间里没有别的素数，因此
$$D(101,1,10)=45.$$
这个值与实现中的校验结果完全一致，也说明了整个方法的计算对象确实只是素数上的模逆。
C++、Python 和 Java 三个实现遵循完全相同的流程。第一步，先筛出所有不超过 \(\lfloor\sqrt{a+b-1}\rfloor\) 的基准素数。第二步，分配一个长度为 \(b\) 的布尔数组，对应整数 \(a,a+1,\dots,a+b-1\)。随后用每个基准素数去标记它在该区间里碰到的所有合数位置。第三步，扫描这个布尔数组，所有未被标记的位置就对应区间素数 \(p\)。
对每个这样的 \(p\)，实现都会运行一次扩展欧几里得算法，输入为 \(k-1\) 和 \(p\)。算法返回的 Bézout 系数经过一次模 \(p\) 归一化后，就变成区间求和中的单项贡献，并被加入总和。实现中的逆元过程写得比较稳健：如果 \(k-1\) 对 \(p\) 取模后为 \(0\)，或者最大公因数不是 \(1\)，就返回 \(0\)。而在目标计算里，真正参与求和的素数都与 \(k-1\) 互素，所以每个有效项都能得到唯一逆元。
记 \(R=a+b-1\)。用普通筛法生成所有不超过 \(\lfloor\sqrt{R}\rfloor\) 的基准素数，需要 \(O(\sqrt{R}\log\log R)\) 时间和 \(O(\sqrt{R})\) 级别的空间。对长度为 \(b\) 的目标区间做分段标记，需要 \(O(b\log\log R)\) 时间和 \(O(b)\) 空间。如果把区间中的素数个数记作 \(\pi([a,a+b))\)，那么所有模逆计算额外带来 \(O(\pi([a,a+b))\log R)\) 时间。因此，整个方法的主导空间复杂度与区间长度 \(b\) 成线性关系，而不是与 \(a+b\) 成线性关系；这正是分段筛在这里真正重要的地方。
Нужно вычислить сумму по простым числам
$$D(a,b,k)=\sum_{\substack{a \le p \lt a+b\\ p\ \text{простое}}} d(p,p-1,k).$$
Ключевой факт, на котором построена реализация, состоит в том, что для каждого простого \(p\), попадающего в суммирование, величина \(d(p,p-1,k)\) сводится к мультипликативному обратному для \(k-1\) по модулю \(p\). После этого задача распадается на две части: быстро перечислить все простые в длинном интервале и для каждого такого простого вычислить один модульный обратный.
После этой редукции арифметическая часть становится довольно прозрачной. Основная инженерная сложность связана не с формулами, а с тем, чтобы обрабатывать интервал \([a,a+b)\) без полного решета до \(a+b\).
Реализация использует тождество
$$d(p,p-1,k)\equiv (k-1)^{-1}\pmod{p}.$$
Значит, для каждого простого \(p\) из интервала нужно найти единственный класс вычетов \(x\), удовлетворяющий
$$(k-1)x \equiv 1 \pmod{p}.$$
Поскольку модуль прост, обратный элемент существует тогда и только тогда, когда \(p\) не делит \(k-1\). Иначе говоря, сам объект \(d(p,p-1,k)\) напрямую строить не нужно; достаточно решить одно линейное сравнение для каждого простого.
Если \(\gcd(k-1,p)=1\), расширенный алгоритм Евклида находит целые числа \(u\) и \(v\), для которых
$$u(k-1)+vp=1.$$
Если перейти по модулю \(p\), второй член исчезает, и получается
$$u(k-1)\equiv 1 \pmod{p}.$$
Следовательно, требуемый вклад равен наименьшему неотрицательному вычету числа \(u\):
$$d(p,p-1,k)=u \bmod p.$$
Именно поэтому реализации на C++, Python и Java используют расширенный алгоритм Евклида, а не модульное возведение в степень. Он напрямую возвращает коэффициент Бéзу, из которого сразу получается обратный элемент.
Интервал может начинаться около \(10^9\), поэтому просеивать все числа от \(2\) до \(a+b\) неэффективно. Обозначим
$$R=a+b-1,\qquad s=\left\lfloor \sqrt{R}\right\rfloor.$$
Сначала обычным решетом находятся все простые числа, не превосходящие \(s\). Это и есть все базовые простые, необходимые для пометки составных чисел в целевом интервале. Для каждого базового простого \(q\) первая кратная величина, которую нужно зачеркнуть внутри \([a,a+b)\), равна
$$m_q=\max\left(q^2,\left\lceil\frac{a}{q}\right\rceil q\right).$$
После этого помечаются числа
$$m_q,\ m_q+q,\ m_q+2q,\dots$$
пока они остаются внутри интервала. Когда все базовые простые обработаны, непомеченные позиции соответствуют ровно тем простым, которые лежат в \([a,a+b)\).
После перечисления всех простых итоговая формула принимает вид
$$D(a,b,k)=\sum_{\substack{a \le p \lt a+b\\ p\ \text{простое}}} \left((k-1)^{-1} \bmod p\right).$$
Никакой дополнительной комбинаторной структуры здесь нет. Алгоритм делает ровно две вещи: находит простые в интервале и ставит в соответствие каждому из них обратный элемент для \(k-1\) по модулю \(p\).
Интервал \([101,102)\) содержит только одно простое число, \(p=101\). Здесь
$$k-1=9.$$
Нужно решить сравнение
$$9x \equiv 1 \pmod{101}.$$
Так как
$$9\cdot 45 = 405 = 4\cdot 101 + 1,$$
обратный к \(9\) по модулю \(101\) равен \(45\). Следовательно,
$$D(101,1,10)=45,$$
что совпадает с контрольным значением реализации.
Реализации на C++, Python и Java используют один и тот же конвейер. Сначала они строят список всех простых до \(\lfloor\sqrt{a+b-1}\rfloor\). Затем выделяется булев сегмент длины \(b\), который представляет числа \(a,a+1,\dots,a+b-1\), и каждое составное число, достигаемое некоторым базовым простым, помечается. После этого каждая непомеченная позиция соответствует простому \(p\) из заданного интервала.
Для каждого такого простого реализация запускает расширенный алгоритм Евклида на паре \(k-1\) и \(p\), нормализует полученный коэффициент Бéзу в диапазон от \(0\) до \(p-1\) и прибавляет его к сумме. Процедура вычисления обратного написана с защитой от вырожденных случаев: если \(k-1\) после редукции равно \(0\) по модулю \(p\), либо если gcd не равен \(1\), возвращается \(0\). Для целевых вычислений все реально участвующие простые взаимно просты с \(k-1\), поэтому каждый корректный член суммы определён однозначно.
Пусть \(R=a+b-1\). Построение всех базовых простых до \(\lfloor\sqrt{R}\rfloor\) с помощью простого решета занимает \(O(\sqrt{R}\log\log R)\) времени и \(O(\sqrt{R})\) памяти. Пометка целевого интервала длины \(b\) требует \(O(b\log\log R)\) времени и \(O(b)\) памяти. Если \(\pi([a,a+b))\) обозначает число простых в интервале, то вычисления обратных элементов добавляют ещё \(O(\pi([a,a+b))\log R)\) времени. Итак, основной расход памяти линейный по длине интервала, а эффективность достигается за счёт того, что алгоритм не хранит все числа до \(a+b\).
المطلوب هو حساب مجموع على الأعداد الأولية:
$$D(a,b,k)=\sum_{\substack{a \le p \lt a+b\\ p\ \text{prime}}} d(p,p-1,k).$$
الفكرة الحاسمة التي تعتمد عليها التنفيذات هي أن كل حد فعلي في هذا المجموع يختزل إلى المعكوس الضربي للعدد \(k-1\) بترديد \(p\). بعد هذه الملاحظة تصبح المسألة عملية وواضحة: نولّد الأعداد الأولية داخل المجال \([a,a+b)\) بكفاءة، ثم نحسب معكوسًا واحدًا لكل عدد أولي ونجمع القيم الناتجة.
بمجرد قبول هذا الاختزال الأساسي، يصبح الجانب النظري مباشرًا نسبيًا. الصعوبة الحقيقية ليست في صياغة علاقة جديدة، بل في معالجة مجال قد يبدأ قرب \(10^9\) من دون غربلة كل الأعداد حتى \(a+b\).
تستخدم التنفيذات الهوية التالية:
$$d(p,p-1,k)\equiv (k-1)^{-1}\pmod{p}.$$
إذن المطلوب لكل عدد أولي \(p\) في المجال هو الباقي الوحيد \(x\) الذي يحقق
$$(k-1)x \equiv 1 \pmod{p}.$$
وبما أن \(p\) أولي، فإن هذا المعكوس يوجد متى كان \(p \nmid (k-1)\). وهذا يعني أن الكمية \(d(p,p-1,k)\) لا تحتاج إلى بناء مباشر؛ يكفي حل توافق خطي واحد لكل عدد أولي.
عندما يكون \(\gcd(k-1,p)=1\)، تعطي خوارزمية إقليدس الموسعة عددين صحيحين \(u\) و\(v\) بحيث
$$u(k-1)+vp=1.$$
وعند أخذ هذه المساواة بترديد \(p\) يختفي الحد الثاني فنحصل على
$$u(k-1)\equiv 1 \pmod{p}.$$
لذلك تكون مساهمة هذا العدد الأولي هي الباقي غير السالب الموافق لـ \(u\):
$$d(p,p-1,k)=u \bmod p.$$
ولهذا السبب تحديدًا تستخدم تنفيذات C++ وPython وJava خوارزمية إقليدس الموسعة بدل الاعتماد على رفع الأسّ بطريقة سريعة. المعكوس يخرج مباشرة من معاملات بéزو، والطريقة نفسها تعمل في اللغات الثلاث بلا فرق جوهري.
لأن المجال قد يبدأ قريبًا من \(10^9\)، فإن تنفيذ غربال كامل من \(2\) حتى \(a+b\) سيكون مكلفًا بلا داع. لنكتب
$$R=a+b-1,\qquad s=\left\lfloor \sqrt{R}\right\rfloor.$$
أولًا نولّد جميع الأعداد الأولية حتى \(s\) باستعمال غربال عادي. هذه هي كل الأوليات الأساسية المطلوبة لتمييز الأعداد المركبة داخل المجال الهدف. ولكل أولي أساسي \(q\)، فإن أول مضاعف يجب شطبه داخل \([a,a+b)\) هو
$$m_q=\max\left(q^2,\left\lceil\frac{a}{q}\right\rceil q\right).$$
بعد ذلك تُعلَّم القيم
$$m_q,\ m_q+q,\ m_q+2q,\dots$$
ما دامت تقع داخل المجال. وبعد معالجة جميع الأوليات الأساسية، تكون المواضع غير المعلّمة هي بالضبط الأعداد الأولية الواقعة في \([a,a+b)\).
بعد معرفة جميع أوليات المجال، تصبح الصيغة النهائية
$$D(a,b,k)=\sum_{\substack{a \le p \lt a+b\\ p\ \text{prime}}} \left((k-1)^{-1} \bmod p\right).$$
لا توجد هنا طبقة خفية إضافية من العدّ أو التركيبات. الخوارزمية تقوم بأمرين فقط: تعثر على الأعداد الأولية في المجال، ثم تربط بكل أولي المعكوس الضربي للعدد \(k-1\) بترديده.
المجال \([101,102)\) يحتوي على عدد أولي واحد فقط هو \(101\). وفي هذه الحالة
$$k-1=9.$$
إذًا نحل التوافق
$$9x \equiv 1 \pmod{101}.$$
ولأن
$$9\cdot 45 = 405 = 4\cdot 101 + 1,$$
فإن معكوس \(9\) بترديد \(101\) هو \(45\). ومن ثم
$$D(101,1,10)=45,$$
وهو نفس مقدار التحقق الذي تستعمله التنفيذات.
تنفيذات C++ وPython وJava تتبع خط سير واحدًا. في البداية تُولَّد جميع الأوليات حتى \(\lfloor\sqrt{a+b-1}\rfloor\). بعد ذلك يُنشأ مقطع منطقي طوله \(b\) يمثل الأعداد \(a,a+1,\dots,a+b-1\)، ثم تُعلَّم كل القيم المركبة التي تصيبها الأوليات الأساسية. بعد انتهاء هذه المرحلة يصبح كل موضع غير معلَّم ممثلًا لعدد أولي \(p\) داخل المجال المطلوب.
لكل عدد أولي من هذه الأعداد، تشغّل التنفيذات خوارزمية إقليدس الموسعة على الزوج \(k-1\) و\(p\)، ثم تُطبِّع معامل بéزو الناتج إلى المجال من \(0\) حتى \(p-1\) وتضيفه إلى المجموع. وروتين المعكوس مكتوب بشكل دفاعي: إذا كان \(k-1\) بعد الاختزال يساوي \(0\) بترديد \(p\)، أو إذا لم يكن القاسم المشترك الأكبر مساويًا لـ \(1\)، فإنه يعيد \(0\). أمّا في التقييمات المقصودة فالأوليات المساهمة كلها أوليات متباينة مع \(k-1\)، لذا يكون كل حد صحيحًا ومحددًا على نحو وحيد.
لنضع \(R=a+b-1\). توليد جميع الأوليات الأساسية حتى \(\lfloor\sqrt{R}\rfloor\) يحتاج إلى \(O(\sqrt{R}\log\log R)\) زمنًا و\(O(\sqrt{R})\) ذاكرة في الغربال البسيط المستخدم هنا. أمّا تعليم المجال ذي الطول \(b\) فيكلف \(O(b\log\log R)\) زمنًا و\(O(b)\) ذاكرة. وإذا رمزنا إلى عدد الأوليات داخل المجال بالرمز \(\pi([a,a+b))\)، فإن حساب المعكوسات يضيف \(O(\pi([a,a+b))\log R)\) زمنًا إضافيًا. وبالتالي فإن الذاكرة المسيطرة خطية في طول المجال، وهذه هي الفائدة الأساسية من الغربال المجزأ: لا حاجة إلى تخزين جميع الأعداد حتى \(a+b\).