Let
$$H_n=\sum_{m=1}^{n}\frac{1}{m}$$
be the nth harmonic number. For an odd prime \(p\), the important exceptional indices are those for which \(H_n\) is divisible by \(p\) in the \(p\)-adic sense. Equivalently, after reducing \(H_n\) to lowest terms, its denominator is coprime to \(p\) and its numerator is divisible by \(p\). Denote this set by \(J_p\). Once \(\max J_p\) is known, the target quantity is
$$M(p)=p\,\max J_p+(p-1).$$
The implementations first verify the checkpoint \(M(7)=719102\), then use the same method for \(p=137\).
The search is performed inside the rings \(\mathbb{Z}/p^r\mathbb{Z}\) rather than with huge rational numbers. That makes the divisibility by \(p\) explicit and lets the algorithm lift surviving indices level by level.
Write the reduced harmonic number as
$$H_n=\frac{A_n}{B_n}, \qquad \gcd(A_n,B_n)=1.$$
If \(p\nmid B_n\), then \(H_n\) has a well-defined residue modulo \(p\). The relevant indices are exactly those for which
$$H_n\in p\mathbb{Z}_p,$$
or equivalently
$$p\mid A_n \qquad \text{and} \qquad p\nmid B_n.$$
These are the members of \(J_p\). The denominator problem is therefore converted into a \(p\)-adic divisibility problem for harmonic numbers.
For \(n\ge 1\) and \(0\le k<p\), split the sum up to \(pn+k\) into multiples of \(p\), non-multiples up to \(pn\), and the final tail:
$$H_{pn+k}=\sum_{a=1}^{n}\frac{1}{ap}+\sum_{\substack{1\le m\le pn \\ p\nmid m}}\frac{1}{m}+\sum_{t=1}^{k}\frac{1}{pn+t}.$$
The first sum is \(H_n/p\), so the exact recurrence is
$$H_{pn+k}=\frac{H_n}{p}+B_p(n)+T_{p,k}(n),$$
where
$$B_p(n)=\sum_{\substack{1\le m\le pn \\ p\nmid m}}\frac{1}{m}, \qquad T_{p,k}(n)=\sum_{t=1}^{k}\frac{1}{pn+t}.$$
This is the backbone of the lifting tree. First compute the value at \(pn\), then generate \(pn+1,\dots,pn+p-1\) by adding one inverse at a time.
The only expensive part is \(B_p(n)\). Rewriting the non-multiples of \(p\) in complete blocks gives
$$B_p(n)=\sum_{a=0}^{n-1}\sum_{t=1}^{p-1}\frac{1}{ap+t}.$$
For fixed precision \(p^S\), each reciprocal admits a \(p\)-adic expansion
$$\frac{1}{ap+t}=\frac{1}{t}\cdot\frac{1}{1+ap/t}=\frac{1}{t}\sum_{j\ge 0}(-1)^j\left(\frac{ap}{t}\right)^j \pmod{p^S}.$$
After summing over the full residue system \(t=1,2,\dots,p-1\), the odd contributions cancel, so modulo \(p^S\) the whole block term lies in the span of even powers of \(n\). Thus there exist coefficients \(c_1,\dots,c_{(S-1)/2}\) such that
$$B_p(n)\equiv \sum_{j=1}^{(S-1)/2} c_j n^{2j} \pmod{p^S}.$$
The implementations do not derive these coefficients symbolically. Instead, they evaluate \(B_p(n)\) for enough small sample values and solve a short linear system modulo \(p^S\).
Suppose a surviving node \(n\) is known together with \(H_n\bmod p^r\), and suppose \(H_n\equiv 0\pmod p\). Then \(H_n/p\) is meaningful modulo \(p^{r-1}\). Using the polynomial block correction, the first child satisfies
$$H_{pn}\equiv \frac{H_n}{p}+B_p(n) \pmod{p^{r-1}}.$$
The remaining children are then generated by cumulative updates
$$H_{pn+k}=H_{pn+k-1}+\frac{1}{pn+k} \pmod{p^{r-1}}, \qquad 1\le k<p.$$
Only children with
$$H_{pn+k}\equiv 0\pmod p$$
survive to the next layer. Each lift consumes one power of \(p\), so a starting precision \(S\) permits at most \(S-1\) meaningful layers.
If \(n\in J_p\), then \(H_n/p\) is still \(p\)-integral, and both \(B_p(n)\) and every short tail \(T_{p,k}(n)\) are sums of \(p\)-adic units. Hence the whole block
$$H_{pn},H_{pn+1},\dots,H_{pn+p-1}$$
remains \(p\)-integral. Therefore the largest index attached to a surviving parent \(n\) is \(pn+(p-1)\). Once \(\max J_p\) is known, the final answer is exactly
$$M(p)=p\,\max J_p+(p-1).$$
For \(p=7\), the seed search only needs \(1\le n<7\). We have
$$H_6=1+\frac12+\frac13+\frac14+\frac15+\frac16=\frac{49}{20},$$
so \(H_6\equiv 0\pmod 7\). Therefore
$$J_7^{(0)}=\{6\}.$$
After one lift the surviving children are
$$J_7^{(1)}=\{42,48\}.$$
The next layers found by the implementation are
$$J_7^{(2)}=\{295,299,337,341\},$$
$$J_7^{(3)}=\{2096,2390\},$$
$$J_7^{(4)}=\{14675,16731,16735\},$$
$$J_7^{(5)}=\{102728\}.$$
No further child survives, so
$$\max J_7=102728.$$
Hence
$$M(7)=7\cdot 102728+6=719102,$$
which matches the validation used by the implementations.
The C++, Python, and Java implementations follow the same sequence. First they fix the precision at \(S=8\) and precompute \(p,p^2,\dots,p^S\). Then they sample the block sum \(B_p(n)\) at a few small values of \(n\) and solve a modular linear system to recover the coefficients of the even polynomial.
Next they compute \(H_1,H_2,\dots,H_{p-1}\) modulo \(p^S\), keeping only the values that are congruent to \(0\pmod p\). This forms the initial frontier. For each surviving node, the implementation evaluates the block polynomial, computes the child \(pn\), scans the remaining \(p-1\) children by successive inverse additions, and keeps only those whose harmonic values remain divisible by \(p\). Throughout the search it updates the largest surviving index.
The C++ implementation additionally parallelizes frontier expansion when the frontier is large, but the arithmetic and survivor test are the same in all three languages.
Let \(L_r\) be the number of surviving nodes at lifting level \(r\), and let \(d=(S-1)/2\) be the number of polynomial coefficients. Recovering the block polynomial costs \(O(d^3)\) modular operations after the sample sums are computed. Since \(S=8\) is fixed, this precomputation is constant-size.
The main cost comes from expanding the survivor tree. Each surviving node requires one polynomial evaluation of cost \(O(d)\) and then a scan of \(p\) children, each using one modular inverse and one modular addition. Therefore the total search cost is
$$O\left(\sum_r L_r(d+p)\right),$$
which for fixed \(S\) is effectively
$$O\left(p\sum_r L_r\right).$$
The memory usage is linear in the largest frontier:
$$O\left(\max_r L_r\right).$$
Sei
$$H_n=\sum_{m=1}^{n}\frac{1}{m}$$
die \(n\)-te harmonische Zahl. Für eine ungerade Primzahl \(p\) sind die entscheidenden Ausnahmen genau jene Indizes, für die \(H_n\) im \(p\)-adischen Sinn durch \(p\) teilbar ist. Äquivalent dazu ist im vollständig gekürzten Bruch der Nenner zu \(p\) teilerfremd, während der Zähler durch \(p\) teilbar ist. Diese Menge nennen wir \(J_p\). Sobald \(\max J_p\) bekannt ist, lautet die gesuchte Größe
$$M(p)=p\,\max J_p+(p-1).$$
Die Implementierungen prüfen zuerst den Kontrollwert \(M(7)=719102\) und behandeln danach mit derselben Methode den Fall \(p=137\).
Die Suche arbeitet nicht mit riesigen rationalen Zahlen, sondern in den Restklassenringen \(\mathbb{Z}/p^r\mathbb{Z}\). Dadurch wird die Teilbarkeit durch \(p\) auf jeder Ebene der Hebung explizit verfolgt.
Schreibe
$$H_n=\frac{A_n}{B_n}, \qquad \gcd(A_n,B_n)=1.$$
Falls \(p\nmid B_n\), besitzt \(H_n\) einen wohldefinierten Rest modulo \(p\). Relevant sind genau die Indizes mit
$$H_n\in p\mathbb{Z}_p,$$
also äquivalent
$$p\mid A_n \qquad \text{und} \qquad p\nmid B_n.$$
Das sind die Elemente von \(J_p\). Damit wird die Frage nach der Nennerteilbarkeit in eine \(p\)-adische Teilbarkeitsfrage für harmonische Zahlen übersetzt.
Für \(n\ge 1\) und \(0\le k<p\) zerlegen wir die Summe bis \(pn+k\) in Vielfache von \(p\), Nichtvielfache bis \(pn\) und die letzten \(k\) Summanden:
$$H_{pn+k}=\sum_{a=1}^{n}\frac{1}{ap}+\sum_{\substack{1\le m\le pn \\ p\nmid m}}\frac{1}{m}+\sum_{t=1}^{k}\frac{1}{pn+t}.$$
Die erste Summe ist \(H_n/p\). Also folgt die exakte Rekursion
$$H_{pn+k}=\frac{H_n}{p}+B_p(n)+T_{p,k}(n),$$
wobei
$$B_p(n)=\sum_{\substack{1\le m\le pn \\ p\nmid m}}\frac{1}{m}, \qquad T_{p,k}(n)=\sum_{t=1}^{k}\frac{1}{pn+t}.$$
Das ist der Kern des Hebungsbaums: zuerst den Wert bei \(pn\) berechnen, danach \(pn+1,\dots,pn+p-1\) durch sukzessives Addieren einzelner Inversen erzeugen.
Der aufwendige Teil ist \(B_p(n)\). In Blockform gilt
$$B_p(n)=\sum_{a=0}^{n-1}\sum_{t=1}^{p-1}\frac{1}{ap+t}.$$
Für feste Präzision \(p^S\) besitzt jeder Kehrwert die \(p\)-adische Entwicklung
$$\frac{1}{ap+t}=\frac{1}{t}\cdot\frac{1}{1+ap/t}=\frac{1}{t}\sum_{j\ge 0}(-1)^j\left(\frac{ap}{t}\right)^j \pmod{p^S}.$$
Summiert man über das vollständige Restsystem \(t=1,2,\dots,p-1\), so heben sich die ungeraden Beiträge weg. Deshalb liegt der gesamte Blockterm modulo \(p^S\) im Spann der geraden Potenzen von \(n\). Also existieren Koeffizienten \(c_1,\dots,c_{(S-1)/2}\) mit
$$B_p(n)\equiv \sum_{j=1}^{(S-1)/2} c_j n^{2j} \pmod{p^S}.$$
Die Implementierungen leiten diese Koeffizienten nicht symbolisch her, sondern bestimmen sie durch Auswertung einiger kleiner Stichproben und durch Lösung eines kleinen linearen Gleichungssystems modulo \(p^S\).
Sei ein überlebender Knoten \(n\) zusammen mit \(H_n\bmod p^r\) bekannt, und sei \(H_n\equiv 0\pmod p\). Dann ist \(H_n/p\) modulo \(p^{r-1}\) sinnvoll. Mit der Polynomdarstellung für \(B_p(n)\) ergibt sich für das erste Kind
$$H_{pn}\equiv \frac{H_n}{p}+B_p(n) \pmod{p^{r-1}}.$$
Die übrigen Kinder entstehen dann durch kumulative Updates
$$H_{pn+k}=H_{pn+k-1}+\frac{1}{pn+k} \pmod{p^{r-1}}, \qquad 1\le k<p.$$
Nur Kinder mit
$$H_{pn+k}\equiv 0\pmod p$$
überleben. Jede Hebung verbraucht genau eine Potenz von \(p\), daher erlaubt eine Startpräzision \(S\) höchstens \(S-1\) aussagekräftige Ebenen.
Für \(n\in J_p\) bleibt \(H_n/p\) weiterhin \(p\)-integral, und sowohl \(B_p(n)\) als auch jeder kurze Nachlauf \(T_{p,k}(n)\) sind Summen \(p\)-adischer Einheiten. Also bleibt der ganze Block
$$H_{pn},H_{pn+1},\dots,H_{pn+p-1}$$
\(p\)-integral. Der größte zu \(n\) gehörige Index ist daher \(pn+(p-1)\). Sobald \(\max J_p\) bekannt ist, erhält man direkt
$$M(p)=p\,\max J_p+(p-1).$$
Für \(p=7\) muss die Startsuche nur \(1\le n<7\) betrachten. Es gilt
$$H_6=1+\frac12+\frac13+\frac14+\frac15+\frac16=\frac{49}{20},$$
also \(H_6\equiv 0\pmod 7\). Damit ist
$$J_7^{(0)}=\{6\}.$$
Nach einem Hebungsschritt bleiben
$$J_7^{(1)}=\{42,48\}.$$
Die folgenden Ebenen sind
$$J_7^{(2)}=\{295,299,337,341\},$$
$$J_7^{(3)}=\{2096,2390\},$$
$$J_7^{(4)}=\{14675,16731,16735\},$$
$$J_7^{(5)}=\{102728\}.$$
Danach überlebt kein weiteres Kind, also
$$\max J_7=102728.$$
Daraus folgt
$$M(7)=7\cdot 102728+6=719102,$$
genau der Kontrollwert aus den Implementierungen.
Die C++-, Python- und Java-Implementierungen folgen demselben Ablauf. Zuerst wird die Präzision \(S=8\) festgelegt und \(p,p^2,\dots,p^S\) werden vorab berechnet. Dann wird der Blockterm \(B_p(n)\) für einige kleine Stichproben ausgewertet, und aus einem modularen linearen Gleichungssystem werden die Koeffizienten des geraden Polynoms gewonnen.
Anschließend berechnet die Implementierung \(H_1,H_2,\dots,H_{p-1}\) modulo \(p^S\) und behält nur die Werte mit Rest \(0\pmod p\). Das ist die erste Front des Hebungsbaums. Für jeden überlebenden Knoten wird die Blockkorrektur ausgewertet, der Wert bei \(pn\) bestimmt, dann werden die restlichen \(p-1\) Kinder durch sukzessive Inversenadditionen erzeugt, und nur bei erneuter Teilbarkeit durch \(p\) bleibt das Kind erhalten. Der größte überlebende Index wird fortlaufend aktualisiert.
Die C++-Version parallelisiert große Fronten zusätzlich, aber Arithmetik und Überlebenstest bleiben in allen drei Sprachen gleich.
Sei \(L_r\) die Anzahl der überlebenden Knoten auf Ebene \(r\), und sei \(d=(S-1)/2\) die Zahl der Polynomkoeffizienten. Das Rekonstruieren des Blockpolynoms kostet nach den Stichproben \(O(d^3)\) modulare Operationen. Wegen \(S=8\) ist diese Vorarbeit konstant klein.
Der Hauptaufwand liegt im Expandieren des Überlebensbaums. Jeder Knoten braucht eine Polynomauswertung der Kosten \(O(d)\) und danach einen Durchlauf über \(p\) Kinder mit jeweils einer modularen Inversion und einer modularen Addition. Insgesamt ergibt sich
$$O\left(\sum_r L_r(d+p)\right),$$
also bei festem \(S\) effektiv
$$O\left(p\sum_r L_r\right).$$
Der Speicherverbrauch ist linear in der größten Front:
$$O\left(\max_r L_r\right).$$
Tanım olarak
$$H_n=\sum_{m=1}^{n}\frac{1}{m}$$
\(H_n\), \(n\). harmonik sayıdır. Tek bir asal \(p\) için önemli olan istisna indisleri, \(H_n\)'nin \(p\)-adik anlamda \(p\)'ye bölündüğü durumlardır. Eşdeğer olarak, sadeleştirilmiş kesirde payda \(p\) ile aralarında asal kalırken pay \(p\)'ye bölünür. Bu kümeyi \(J_p\) ile gösterelim. \(\max J_p\) belirlendiğinde hedef büyüklük
$$M(p)=p\,\max J_p+(p-1).$$
Uygulamalar önce \(M(7)=719102\) doğrulamasını yapar, sonra aynı yöntemi \(p=137\) için uygular.
Çözüm dev rasyonel sayıları doğrudan üretmek yerine \(\mathbb{Z}/p^r\mathbb{Z}\) halkalarında çalışır. Böylece \(p\)'ye bölünebilirlik, kaldırma ağacı boyunca katman katman izlenebilir.
Harmonik sayıyı
$$H_n=\frac{A_n}{B_n}, \qquad \gcd(A_n,B_n)=1$$
şeklinde sade halde yazalım. Eğer \(p\nmid B_n\) ise \(H_n\)'nin modulo \(p\) kalanı anlamlıdır. İlgilendiğimiz indisler tam olarak
$$H_n\in p\mathbb{Z}_p,$$
yani eşdeğer olarak
$$p\mid A_n \qquad \text{ve} \qquad p\nmid B_n$$
koşulunu sağlayanlardır. Bunlar \(J_p\)'yi oluşturur. Böylece payda bölünebilirliği sorunu, harmonik sayının \(p\)-adik bütün kalması ve ayrıca \(p\)'ye bölünmesi sorununa dönüşür.
Her \(n\ge 1\) ve \(0\le k<p\) için toplamı üç parçaya ayıralım:
$$H_{pn+k}=\sum_{a=1}^{n}\frac{1}{ap}+\sum_{\substack{1\le m\le pn \\ p\nmid m}}\frac{1}{m}+\sum_{t=1}^{k}\frac{1}{pn+t}.$$
İlk parça \(H_n/p\) olduğundan temel bağıntı
$$H_{pn+k}=\frac{H_n}{p}+B_p(n)+T_{p,k}(n),$$
olur; burada
$$B_p(n)=\sum_{\substack{1\le m\le pn \\ p\nmid m}}\frac{1}{m}, \qquad T_{p,k}(n)=\sum_{t=1}^{k}\frac{1}{pn+t}.$$
Bu eşitlik kaldırma ağacının merkezidir. Önce \(pn\) noktasındaki değer hesaplanır, sonra \(pn+1,\dots,pn+p-1\) çocukları birer ters eklenerek taranır.
Zor kısım \(B_p(n)\) blok toplamıdır. Bunu
$$B_p(n)=\sum_{a=0}^{n-1}\sum_{t=1}^{p-1}\frac{1}{ap+t}$$
şeklinde yazabiliriz. Sabit \(p^S\) hassasiyeti altında her ters için
$$\frac{1}{ap+t}=\frac{1}{t}\cdot\frac{1}{1+ap/t}=\frac{1}{t}\sum_{j\ge 0}(-1)^j\left(\frac{ap}{t}\right)^j \pmod{p^S}$$
açılımı geçerlidir. \(t=1,2,\dots,p-1\) tam kalan sistemi üzerinde toplandığında tek dereceli katkılar yok olur. Bu nedenle blok toplamı modulo \(p^S\), \(n\)'in yalnızca çift kuvvetlerinin doğrusal birleşimi olur. Yani bazı \(c_1,\dots,c_{(S-1)/2}\) katsayıları için
$$B_p(n)\equiv \sum_{j=1}^{(S-1)/2} c_j n^{2j} \pmod{p^S}$$
yazabiliriz. Uygulamalar bu katsayıları kapalı formdan değil, küçük örnek değerlerden kurulan modüler lineer sistemden çıkarır.
Bir düğüm \(n\) için \(H_n\bmod p^r\) biliniyor ve ayrıca \(H_n\equiv 0\pmod p\) olsun. O zaman \(H_n/p\), \(p^{r-1}\) modunda anlamlıdır. Polinom blok düzeltmesi ile ilk çocuk için
$$H_{pn}\equiv \frac{H_n}{p}+B_p(n) \pmod{p^{r-1}}$$
elde edilir. Kalan çocuklar ise
$$H_{pn+k}=H_{pn+k-1}+\frac{1}{pn+k} \pmod{p^{r-1}}, \qquad 1\le k<p$$
güncellemesiyle bulunur. Yalnızca
$$H_{pn+k}\equiv 0\pmod p$$
koşulunu koruyan çocuklar bir sonraki seviyeye taşınır. Her kaldırma adımı bir \(p\) kuvveti tüketir, dolayısıyla başlangıç hassasiyeti \(S\) ise en fazla \(S-1\) anlamlı seviye vardır.
Eğer \(n\in J_p\) ise \(H_n/p\) yine \(p\)-adik tam sayıdır; ayrıca \(B_p(n)\) ve her kısa kuyruk \(T_{p,k}(n)\) de \(p\)-adik birimlerin sonlu toplamıdır. Böylece
$$H_{pn},H_{pn+1},\dots,H_{pn+p-1}$$
değerlerinin tamamı \(p\)-adik bütün kalır. Bu yüzden \(n\)'den doğan bloktaki en büyük indis \(pn+(p-1)\) olur. \(\max J_p\) bulunduğu anda sonuç doğrudan
$$M(p)=p\,\max J_p+(p-1)$$
şeklindedir.
\(p=7\) için başlangıç araması yalnızca \(1\le n<7\) aralığındadır. Burada
$$H_6=1+\frac12+\frac13+\frac14+\frac15+\frac16=\frac{49}{20},$$
dolayısıyla \(H_6\equiv 0\pmod 7\) olur. İlk küme
$$J_7^{(0)}=\{6\}$$
şeklindedir. İlk kaldırma sonunda
$$J_7^{(1)}=\{42,48\}$$
kalır. Sonraki seviyeler
$$J_7^{(2)}=\{295,299,337,341\},$$
$$J_7^{(3)}=\{2096,2390\},$$
$$J_7^{(4)}=\{14675,16731,16735\},$$
$$J_7^{(5)}=\{102728\}$$
olur. Daha sonra yeni çocuk çıkmaz, yani
$$\max J_7=102728.$$
Buradan
$$M(7)=7\cdot 102728+6=719102$$
elde edilir; bu, uygulamaların kullandığı doğrulama değeridir.
C++, Python ve Java uygulamaları aynı hattı izler. Önce hassasiyet \(S=8\) seçilir ve \(p,p^2,\dots,p^S\) kuvvetleri hazırlanır. Ardından \(B_p(n)\) blok toplamı birkaç küçük örnek üzerinde hesaplanır ve çift kuvvetli polinomun katsayıları modüler bir lineer sistem çözülerek bulunur.
Sonra \(H_1,H_2,\dots,H_{p-1}\) değerleri modulo \(p^S\) hesaplanır; yalnızca \(0\pmod p\) verenler ilk sınır kümesine alınır. Her yaşayan düğüm için blok polinomu değerlendirilir, \(pn\) çocuğu elde edilir, kalan \(p-1\) çocuk ardışık ters eklemeleriyle taranır ve sadece harmonik toplamı tekrar \(p\)'ye bölünen çocuklar tutulur. Arama süresince en büyük yaşayan indis güncellenir.
C++ sürümü büyük sınır kümelerinde paralel genişletme de yapar, ancak aritmetik ve eleme ölçütü üç dilde aynıdır.
\(L_r\), \(r\). kaldırma seviyesindeki yaşayan düğüm sayısı ve \(d=(S-1)/2\), polinom katsayı sayısı olsun. Blok polinomunu kurmak için gereken maliyet örnekler hazırlandıktan sonra \(O(d^3)\) modüler işlemdir. \(S=8\) sabit olduğundan bu ön hesaplama sabit boyutludur.
Ana maliyet yaşayan ağaç düğümlerini genişletmektir. Her düğüm için \(O(d)\) maliyetli bir polinom değerlendirmesi ve ardından \(p\) çocuk üzerinde dolaşma gerekir; her çocukta bir modüler ters ve bir modüler toplama vardır. Toplam süre
$$O\left(\sum_r L_r(d+p)\right)$$
olur; sabit \(S\) altında bu pratikte
$$O\left(p\sum_r L_r\right)$$
şeklindedir. Bellek kullanımı ise en büyük sınır kümesinin boyutuyla orantılıdır:
$$O\left(\max_r L_r\right).$$
Sea
$$H_n=\sum_{m=1}^{n}\frac{1}{m}$$
el \(n\)-ésimo número armónico. Para un primo impar \(p\), los índices excepcionales son aquellos para los que \(H_n\) es divisible por \(p\) en sentido \(p\)-ádico. De forma equivalente, al reducir \(H_n\) a una fracción irreducible, el denominador queda coprimo con \(p\) y el numerador sí es múltiplo de \(p\). Denotemos ese conjunto por \(J_p\). Una vez conocido \(\max J_p\), la cantidad buscada es
$$M(p)=p\,\max J_p+(p-1).$$
Las implementaciones validan primero el valor \(M(7)=719102\) y después aplican el mismo procedimiento para \(p=137\).
En vez de manipular directamente fracciones enormes, el método trabaja dentro de \(\mathbb{Z}/p^r\mathbb{Z}\). Así la divisibilidad por \(p\) puede seguirse capa por capa dentro del árbol de lifting.
Escribimos
$$H_n=\frac{A_n}{B_n}, \qquad \gcd(A_n,B_n)=1.$$
Si \(p\nmid B_n\), entonces \(H_n\) tiene residuo bien definido módulo \(p\). Los índices relevantes son exactamente los que cumplen
$$H_n\in p\mathbb{Z}_p,$$
o, de forma equivalente,
$$p\mid A_n \qquad \text{y} \qquad p\nmid B_n.$$
Esos índices forman \(J_p\). Por tanto, la cuestión sobre el denominador se convierte en una cuestión de divisibilidad \(p\)-ádica del propio número armónico.
Para \(n\ge 1\) y \(0\le k<p\), descomponemos la suma hasta \(pn+k\) en múltiplos de \(p\), no múltiplos hasta \(pn\), y los últimos \(k\) términos:
$$H_{pn+k}=\sum_{a=1}^{n}\frac{1}{ap}+\sum_{\substack{1\le m\le pn \\ p\nmid m}}\frac{1}{m}+\sum_{t=1}^{k}\frac{1}{pn+t}.$$
La primera suma es \(H_n/p\), y por eso obtenemos la recurrencia exacta
$$H_{pn+k}=\frac{H_n}{p}+B_p(n)+T_{p,k}(n),$$
donde
$$B_p(n)=\sum_{\substack{1\le m\le pn \\ p\nmid m}}\frac{1}{m}, \qquad T_{p,k}(n)=\sum_{t=1}^{k}\frac{1}{pn+t}.$$
Esta identidad es el esqueleto del lifting: primero se calcula el valor en \(pn\) y luego se recorre \(pn+1,\dots,pn+p-1\) añadiendo un inverso cada vez.
La parte costosa es \(B_p(n)\). Reordenando por bloques completos,
$$B_p(n)=\sum_{a=0}^{n-1}\sum_{t=1}^{p-1}\frac{1}{ap+t}.$$
Para precisión fija \(p^S\), cada recíproco tiene la expansión \(p\)-ádica
$$\frac{1}{ap+t}=\frac{1}{t}\cdot\frac{1}{1+ap/t}=\frac{1}{t}\sum_{j\ge 0}(-1)^j\left(\frac{ap}{t}\right)^j \pmod{p^S}.$$
Al sumar sobre el sistema completo \(t=1,2,\dots,p-1\), las contribuciones impares se cancelan. Por ello, módulo \(p^S\), el término de bloque queda dentro del espacio generado por las potencias pares de \(n\). Existen entonces coeficientes \(c_1,\dots,c_{(S-1)/2}\) tales que
$$B_p(n)\equiv \sum_{j=1}^{(S-1)/2} c_j n^{2j} \pmod{p^S}.$$
Las implementaciones no usan una fórmula cerrada para esos coeficientes. Los reconstruyen evaluando \(B_p(n)\) en unas pocas muestras y resolviendo un pequeño sistema lineal modular.
Supongamos que un nodo \(n\) ya es conocido junto con \(H_n\bmod p^r\), y además \(H_n\equiv 0\pmod p\). Entonces \(H_n/p\) tiene sentido módulo \(p^{r-1}\). Usando la corrección polinómica, el primer hijo cumple
$$H_{pn}\equiv \frac{H_n}{p}+B_p(n) \pmod{p^{r-1}}.$$
Los demás hijos se generan con actualizaciones acumulativas
$$H_{pn+k}=H_{pn+k-1}+\frac{1}{pn+k} \pmod{p^{r-1}}, \qquad 1\le k<p.$$
Solo sobreviven los hijos que conservan
$$H_{pn+k}\equiv 0\pmod p.$$
Cada lifting consume una potencia de \(p\), de modo que una precisión inicial \(S\) permite como máximo \(S-1\) niveles informativos.
Si \(n\in J_p\), entonces \(H_n/p\) sigue siendo \(p\)-integral, y tanto \(B_p(n)\) como cada cola corta \(T_{p,k}(n)\) son sumas de unidades \(p\)-ádicas. Por eso, todo el bloque
$$H_{pn},H_{pn+1},\dots,H_{pn+p-1}$$
permanece \(p\)-integral. Así, el mayor índice asociado al padre \(n\) es \(pn+(p-1)\). Una vez hallado \(\max J_p\), la respuesta es
$$M(p)=p\,\max J_p+(p-1).$$
Para \(p=7\), la búsqueda inicial solo recorre \(1\le n<7\). Se obtiene
$$H_6=1+\frac12+\frac13+\frac14+\frac15+\frac16=\frac{49}{20},$$
luego \(H_6\equiv 0\pmod 7\). Por tanto,
$$J_7^{(0)}=\{6\}.$$
Tras un paso de lifting sobreviven
$$J_7^{(1)}=\{42,48\}.$$
Las capas siguientes son
$$J_7^{(2)}=\{295,299,337,341\},$$
$$J_7^{(3)}=\{2096,2390\},$$
$$J_7^{(4)}=\{14675,16731,16735\},$$
$$J_7^{(5)}=\{102728\}.$$
Después no aparece ningún hijo nuevo, así que
$$\max J_7=102728.$$
En consecuencia,
$$M(7)=7\cdot 102728+6=719102,$$
exactamente el valor de validación usado por las implementaciones.
Las implementaciones en C++, Python y Java siguen el mismo flujo. Primero fijan la precisión en \(S=8\) y precalculan \(p,p^2,\dots,p^S\). Después evalúan \(B_p(n)\) en unas pocas muestras pequeñas y resuelven un sistema lineal modular para reconstruir los coeficientes del polinomio par.
A continuación calculan \(H_1,H_2,\dots,H_{p-1}\) módulo \(p^S\) y conservan solo los valores congruentes con \(0\pmod p\). Esa es la primera frontera. Para cada nodo superviviente, la implementación evalúa el polinomio de bloque, construye el valor en \(pn\), recorre los \(p-1\) hijos restantes mediante sumas sucesivas de inversos y retiene solo los hijos cuyo valor armónico sigue siendo divisible por \(p\). Durante todo el proceso mantiene actualizado el mayor índice superviviente.
La versión en C++ además paraleliza la expansión de la frontera cuando esta es grande, pero la aritmética y el criterio de supervivencia son los mismos en los tres lenguajes.
Sea \(L_r\) el número de nodos supervivientes en el nivel \(r\), y sea \(d=(S-1)/2\) el número de coeficientes del polinomio. Reconstruir el polinomio de bloque cuesta \(O(d^3)\) operaciones modulares una vez calculadas las muestras. Como \(S=8\) es fijo, esta parte es de tamaño constante.
El coste principal viene de expandir el árbol de supervivientes. Cada nodo requiere una evaluación polinómica de coste \(O(d)\) y un barrido de \(p\) hijos, cada uno con una inversión modular y una suma modular. Por tanto, el coste total es
$$O\left(\sum_r L_r(d+p)\right),$$
que para \(S\) fijo se simplifica en la práctica a
$$O\left(p\sum_r L_r\right).$$
La memoria usada es lineal en el tamaño de la frontera más grande:
$$O\left(\max_r L_r\right).$$
设
$$H_n=\sum_{m=1}^{n}\frac{1}{m}$$
为第 \(n\) 个调和数。对于奇素数 \(p\)，关键的例外下标是那些使 \(H_n\) 在 \(p\)-进意义下被 \(p\) 整除的 \(n\)。等价地，把 \(H_n\) 约成最简分数后，它的分母与 \(p\) 互素，而分子能被 \(p\) 整除。把这些下标组成的集合记为 \(J_p\)。一旦求出 \(\max J_p\)，目标量就是
$$M(p)=p\,\max J_p+(p-1).$$
实现先用 \(M(7)=719102\) 做校验，再把同一套方法用于 \(p=137\)。
整个算法并不直接处理巨大的有理数，而是在 \(\mathbb{Z}/p^r\mathbb{Z}\) 中逐层推进。这样，“是否还被 \(p\) 整除”就能沿着提升树被稳定地追踪。
把调和数写成最简形式
$$H_n=\frac{A_n}{B_n}, \qquad \gcd(A_n,B_n)=1.$$
如果 \(p\nmid B_n\)，那么 \(H_n\) 在模 \(p\) 意义下有明确的值。我们关心的正是满足
$$H_n\in p\mathbb{Z}_p$$
的下标，也就是等价地满足
$$p\mid A_n \qquad \text{且} \qquad p\nmid B_n$$
的下标。这些 \(n\) 就构成 \(J_p\)。因此，题目中“分母是否含有素因子 \(p\)”的问题，被改写成了调和数的 \(p\)-进整性与 \(p\)-进可整除性问题。
对任意 \(n\ge 1\) 和 \(0\le k<p\)，把到 \(pn+k\) 为止的和分成三部分：所有 \(p\) 的倍数项、不被 \(p\) 整除且不超过 \(pn\) 的项、以及最后的 \(k\) 项：
$$H_{pn+k}=\sum_{a=1}^{n}\frac{1}{ap}+\sum_{\substack{1\le m\le pn \\ p\nmid m}}\frac{1}{m}+\sum_{t=1}^{k}\frac{1}{pn+t}.$$
第一部分恰好等于 \(H_n/p\)，于是得到精确递推
$$H_{pn+k}=\frac{H_n}{p}+B_p(n)+T_{p,k}(n),$$
其中
$$B_p(n)=\sum_{\substack{1\le m\le pn \\ p\nmid m}}\frac{1}{m}, \qquad T_{p,k}(n)=\sum_{t=1}^{k}\frac{1}{pn+t}.$$
这就是提升树的骨架。先算出 \(pn\) 处的值，再通过不断加上一个倒数，依次扫描 \(pn+1,\dots,pn+p-1\)。
真正昂贵的是块和 \(B_p(n)\)。把非 \(p\) 倍数按块重排后，有
$$B_p(n)=\sum_{a=0}^{n-1}\sum_{t=1}^{p-1}\frac{1}{ap+t}.$$
在固定精度 \(p^S\) 下，每个倒数都可写成 \(p\)-进展开
$$\frac{1}{ap+t}=\frac{1}{t}\cdot\frac{1}{1+ap/t}=\frac{1}{t}\sum_{j\ge 0}(-1)^j\left(\frac{ap}{t}\right)^j \pmod{p^S}.$$
当对完整剩余系 \(t=1,2,\dots,p-1\) 求和时，奇次项会互相抵消，所以最终只剩偶次贡献。于是模 \(p^S\) 下，整个块修正项可以写成 \(n\) 的偶次幂线性组合。也就是说，存在系数 \(c_1,\dots,c_{(S-1)/2}\)，使得
$$B_p(n)\equiv \sum_{j=1}^{(S-1)/2} c_j n^{2j} \pmod{p^S}.$$
实现并没有把这些系数的闭式写死，而是先对若干个小样本 \(n\) 直接计算 \(B_p(n)\)，再在模 \(p^S\) 下解一个小型线性方程组，把这个偶多项式恢复出来。
设某个幸存节点 \(n\) 已知，并且知道 \(H_n\bmod p^r\)，同时满足 \(H_n\equiv 0\pmod p\)。那么 \(H_n/p\) 在模 \(p^{r-1}\) 下有意义。利用上面的块多项式，可以先得到第一个子节点
$$H_{pn}\equiv \frac{H_n}{p}+B_p(n) \pmod{p^{r-1}}.$$
随后其余子节点通过累加生成：
$$H_{pn+k}=H_{pn+k-1}+\frac{1}{pn+k} \pmod{p^{r-1}}, \qquad 1\le k<p.$$
只有仍满足
$$H_{pn+k}\equiv 0\pmod p$$
的子节点才会进入下一层。每做一次提升，就会少一位 \(p\) 的精度，因此初始精度 \(S\) 最多支持 \(S-1\) 层真正有效的扩展。
如果 \(n\in J_p\)，那么 \(H_n/p\) 仍然是 \(p\)-进整数，而 \(B_p(n)\) 与每个短尾项 \(T_{p,k}(n)\) 也都是 \(p\)-进单位的有限和。因此整个区间
$$H_{pn},H_{pn+1},\dots,H_{pn+p-1}$$
都会保持 \(p\)-进整性。于是由父节点 \(n\) 生成的最大下标就是 \(pn+(p-1)\)。一旦找到 \(\max J_p\)，最终答案立刻写成
$$M(p)=p\,\max J_p+(p-1).$$
当 \(p=7\) 时，种子搜索只需要检查 \(1\le n<7\)。这时
$$H_6=1+\frac12+\frac13+\frac14+\frac15+\frac16=\frac{49}{20},$$
所以 \(H_6\equiv 0\pmod 7\)，从而
$$J_7^{(0)}=\{6\}.$$
做一次提升后，只剩下
$$J_7^{(1)}=\{42,48\}.$$
再往后，程序得到的各层为
$$J_7^{(2)}=\{295,299,337,341\},$$
$$J_7^{(3)}=\{2096,2390\},$$
$$J_7^{(4)}=\{14675,16731,16735\},$$
$$J_7^{(5)}=\{102728\}.$$
此后已经没有新的幸存子节点，因此
$$\max J_7=102728.$$
于是
$$M(7)=7\cdot 102728+6=719102,$$
这正是实现所使用的校验值。
C++、Python 和 Java 三个实现遵循同一条流水线。首先固定 \(S=8\) 的 \(p\)-进精度，并预先计算 \(p,p^2,\dots,p^S\)。然后程序对若干个很小的样本 \(n\) 直接求出 \(B_p(n)\)，再通过一个模 \(p^S\) 的小型线性系统恢复偶多项式的系数。
接下来，程序计算 \(H_1,H_2,\dots,H_{p-1}\) 在模 \(p^S\) 下的值，只保留那些模 \(p\) 余数为 \(0\) 的项，作为提升树的第一层前沿。对于每一个幸存节点，程序先用多项式求块修正，再得到 \(pn\) 处的值，随后通过连续加上 \((pn+1)^{-1},(pn+2)^{-1},\dots\) 扫描其余 \(p-1\) 个子节点，只有仍然满足模 \(p\) 可整除的节点才会进入下一层。搜索过程中会持续更新当前的最大幸存下标。
C++ 版本在前沿较大时还会并行展开节点，但三种语言的算术规则与筛选条件完全相同。
记第 \(r\) 层幸存节点数为 \(L_r\)，记多项式项数为 \(d=(S-1)/2\)。在样本值已经得到之后，恢复块多项式系数需要 \(O(d^3)\) 次模运算。由于这里固定 \(S=8\)，这部分预处理是常数级的。
主要开销来自幸存树的扩展。每个节点都需要一次 \(O(d)\) 的多项式求值，随后再扫描 \(p\) 个候选子节点；每个子节点对应一次模逆和一次模加。因此总时间复杂度为
$$O\left(\sum_r L_r(d+p)\right),$$
在 \(S\) 固定时，可视为
$$O\left(p\sum_r L_r\right).$$
空间复杂度由最大前沿大小决定：
$$O\left(\max_r L_r\right).$$
Пусть
$$H_n=\sum_{m=1}^{n}\frac{1}{m}$$
обозначает \(n\)-е гармоническое число. Для нечётного простого \(p\) важны те индексы, для которых \(H_n\) делится на \(p\) в \(p\)-адическом смысле. Эквивалентно, после сокращения дроби знаменатель оказывается взаимно простым с \(p\), а числитель делится на \(p\). Это множество обозначается \(J_p\). Когда найдено \(\max J_p\), искомая величина равна
$$M(p)=p\,\max J_p+(p-1).$$
Реализации сначала проверяют значение \(M(7)=719102\), а затем тем же методом обрабатывают случай \(p=137\).
Вместо прямой работы с огромными рациональными числами решение переходит в кольца \(\mathbb{Z}/p^r\mathbb{Z}\). Это позволяет отслеживать делимость на \(p\) по уровням \(p\)-адической точности.
Запишем
$$H_n=\frac{A_n}{B_n}, \qquad \gcd(A_n,B_n)=1.$$
Если \(p\nmid B_n\), то \(H_n\) корректно рассматривается по модулю \(p\). Нас интересуют именно те индексы, для которых
$$H_n\in p\mathbb{Z}_p,$$
то есть эквивалентно
$$p\mid A_n \qquad \text{и} \qquad p\nmid B_n.$$
Именно они составляют \(J_p\). Так вопрос о делимости знаменателя переводится в вопрос о \(p\)-адической делимости самого гармонического числа.
Для любых \(n\ge 1\) и \(0\le k<p\) разложим сумму до \(pn+k\) на кратные \(p\), некратные \(p\) до \(pn\) и последние \(k\) слагаемых:
$$H_{pn+k}=\sum_{a=1}^{n}\frac{1}{ap}+\sum_{\substack{1\le m\le pn \\ p\nmid m}}\frac{1}{m}+\sum_{t=1}^{k}\frac{1}{pn+t}.$$
Первая сумма равна \(H_n/p\), поэтому получаем точную рекурсию
$$H_{pn+k}=\frac{H_n}{p}+B_p(n)+T_{p,k}(n),$$
где
$$B_p(n)=\sum_{\substack{1\le m\le pn \\ p\nmid m}}\frac{1}{m}, \qquad T_{p,k}(n)=\sum_{t=1}^{k}\frac{1}{pn+t}.$$
Это основа дерева подъёма: сначала вычисляется значение в точке \(pn\), затем узлы \(pn+1,\dots,pn+p-1\) получаются последовательным добавлением одного обратного элемента.
Самая дорогая часть вычислений связана с \(B_p(n)\). После перегруппировки по блокам имеем
$$B_p(n)=\sum_{a=0}^{n-1}\sum_{t=1}^{p-1}\frac{1}{ap+t}.$$
При фиксированной точности \(p^S\) каждый обратный элемент имеет \(p\)-адическое разложение
$$\frac{1}{ap+t}=\frac{1}{t}\cdot\frac{1}{1+ap/t}=\frac{1}{t}\sum_{j\ge 0}(-1)^j\left(\frac{ap}{t}\right)^j \pmod{p^S}.$$
После суммирования по полному набору остатков \(t=1,2,\dots,p-1\) нечётные вклады исчезают. Поэтому весь блочный член по модулю \(p^S\) выражается через чётные степени \(n\). Значит, существуют коэффициенты \(c_1,\dots,c_{(S-1)/2}\), для которых
$$B_p(n)\equiv \sum_{j=1}^{(S-1)/2} c_j n^{2j} \pmod{p^S}.$$
Реализации не используют заранее выписанную формулу этих коэффициентов. Они вычисляют \(B_p(n)\) на нескольких малых образцах и затем восстанавливают коэффициенты из небольшого линейного модульного уравнения.
Пусть известен выживший узел \(n\) вместе с \(H_n\bmod p^r\), причём \(H_n\equiv 0\pmod p\). Тогда величина \(H_n/p\) корректно определена по модулю \(p^{r-1}\). Используя блочный многочлен, получаем для первого потомка
$$H_{pn}\equiv \frac{H_n}{p}+B_p(n) \pmod{p^{r-1}}.$$
Остальные потомки строятся накопительными обновлениями
$$H_{pn+k}=H_{pn+k-1}+\frac{1}{pn+k} \pmod{p^{r-1}}, \qquad 1\le k<p.$$
В следующий слой переходят только те потомки, которые сохраняют условие
$$H_{pn+k}\equiv 0\pmod p.$$
Каждый шаг подъёма уменьшает доступную точность на одну степень \(p\), поэтому стартовая точность \(S\) даёт не более \(S-1\) содержательных уровней.
Если \(n\in J_p\), то \(H_n/p\) остаётся \(p\)-адически целым, а величины \(B_p(n)\) и все короткие хвосты \(T_{p,k}(n)\) представляют собой суммы \(p\)-адических единиц. Следовательно, весь блок
$$H_{pn},H_{pn+1},\dots,H_{pn+p-1}$$
остаётся \(p\)-адически целым. Наибольший индекс в блоке, порождённом \(n\), равен \(pn+(p-1)\). Поэтому после нахождения \(\max J_p\) сразу получаем
$$M(p)=p\,\max J_p+(p-1).$$
Для \(p=7\) начальный поиск ограничен диапазоном \(1\le n<7\). Имеем
$$H_6=1+\frac12+\frac13+\frac14+\frac15+\frac16=\frac{49}{20},$$
откуда \(H_6\equiv 0\pmod 7\). Следовательно,
$$J_7^{(0)}=\{6\}.$$
После первого шага подъёма остаются
$$J_7^{(1)}=\{42,48\}.$$
Следующие уровни таковы:
$$J_7^{(2)}=\{295,299,337,341\},$$
$$J_7^{(3)}=\{2096,2390\},$$
$$J_7^{(4)}=\{14675,16731,16735\},$$
$$J_7^{(5)}=\{102728\}.$$
Дальше потомков уже нет, значит
$$\max J_7=102728.$$
Тем самым
$$M(7)=7\cdot 102728+6=719102,$$
что точно совпадает с проверочным значением в реализациях.
Реализации на C++, Python и Java используют один и тот же конвейер. Сначала фиксируется точность \(S=8\), после чего заранее строятся степени \(p,p^2,\dots,p^S\). Затем программа вычисляет \(B_p(n)\) на нескольких малых образцах и восстанавливает коэффициенты чётного многочлена, решая небольшую линейную систему по модулю \(p^S\).
После этого вычисляются стартовые значения \(H_1,H_2,\dots,H_{p-1}\) по модулю \(p^S\); сохраняются только те, что дают остаток \(0\pmod p\). Это первая граница дерева подъёма. Для каждого выжившего узла программа вычисляет блочную поправку, получает значение в \(pn\), затем проходит по оставшимся \(p-1\) потомкам последовательным добавлением обратных элементов и оставляет только те узлы, у которых гармоническое число снова делится на \(p\). Максимальный найденный индекс поддерживается на всём протяжении поиска.
Версия на C++ дополнительно распараллеливает расширение большой границы, однако арифметика и критерий отбора во всех трёх языках полностью совпадают.
Пусть \(L_r\) обозначает число выживших узлов на уровне \(r\), а \(d=(S-1)/2\) — число коэффициентов многочлена. Восстановление блочного многочлена после вычисления образцов требует \(O(d^3)\) модульных операций. Так как \(S=8\) фиксировано, эта подготовка имеет постоянный размер.
Основная стоимость сосредоточена в расширении дерева выживших узлов. Для каждого узла требуется одно вычисление многочлена стоимости \(O(d)\), а затем просмотр \(p\) потомков, каждый из которых использует одну модульную инверсию и одно модульное сложение. Поэтому суммарная стоимость поиска равна
$$O\left(\sum_r L_r(d+p)\right),$$
что при фиксированном \(S\) практически сводится к
$$O\left(p\sum_r L_r\right).$$
Память линейна по размеру наибольшей границы:
$$O\left(\max_r L_r\right).$$
ليكن
$$H_n=\sum_{m=1}^{n}\frac{1}{m}$$
هو العدد التوافقي رقم \(n\). بالنسبة إلى عدد أولي فردي \(p\)، فإن الفهارس الاستثنائية هي تلك التي يكون فيها \(H_n\) قابلًا للقسمة على \(p\) بالمعنى \(p\)-adic. وبصورة مكافئة، إذا كُتب \(H_n\) في صورة كسر مختزل فإن مقامه يكون أوليًا نسبيًا مع \(p\)، بينما يقبل البسط القسمة على \(p\). نرمز لهذه المجموعة بـ \(J_p\). وعند معرفة \(\max J_p\) تصبح الكمية المطلوبة
$$M(p)=p\,\max J_p+(p-1).$$
تتحقق التطبيقات أولًا من القيمة \(M(7)=719102\)، ثم تطبق المنهج نفسه على \(p=137\).
بدلًا من تكوين كسور ضخمة مباشرة، تعمل الطريقة داخل الحلقات \(\mathbb{Z}/p^r\mathbb{Z}\). بهذه الصياغة يمكن تتبع قابلية القسمة على \(p\) مستوى بعد مستوى داخل شجرة الرفع.
لنكتب
$$H_n=\frac{A_n}{B_n}, \qquad \gcd(A_n,B_n)=1.$$
إذا كان \(p\nmid B_n\)، فإن لـ \(H_n\) معنى واضحًا بترديد \(p\). الفهارس المهمة هي بالضبط التي تحقق
$$H_n\in p\mathbb{Z}_p,$$
أو بصورة مكافئة
$$p\mid A_n \qquad \text{و} \qquad p\nmid B_n.$$
هذه الفهارس هي عناصر \(J_p\). وهكذا تتحول مسألة قسمة المقام إلى مسألة قسمة \(p\)-adic للعدد التوافقي نفسه.
لكل \(n\ge 1\) و\(0\le k<p\)، نقسم المجموع حتى \(pn+k\) إلى حدود مضاعفات \(p\)، وحدود غير المضاعفات حتى \(pn\)، ثم آخر \(k\) حدود:
$$H_{pn+k}=\sum_{a=1}^{n}\frac{1}{ap}+\sum_{\substack{1\le m\le pn \\ p\nmid m}}\frac{1}{m}+\sum_{t=1}^{k}\frac{1}{pn+t}.$$
المجموع الأول يساوي تمامًا \(H_n/p\)، ومن ثم نحصل على العلاقة الدقيقة
$$H_{pn+k}=\frac{H_n}{p}+B_p(n)+T_{p,k}(n),$$
حيث
$$B_p(n)=\sum_{\substack{1\le m\le pn \\ p\nmid m}}\frac{1}{m}, \qquad T_{p,k}(n)=\sum_{t=1}^{k}\frac{1}{pn+t}.$$
هذه هي البنية الأساسية لشجرة الرفع: نحسب أولًا القيمة عند \(pn\)، ثم نمر على الأبناء \(pn+1,\dots,pn+p-1\) بإضافة مقلوب واحد في كل خطوة.
الجزء المكلف هو \(B_p(n)\). وبعد إعادة الترتيب على شكل كتل كاملة نحصل على
$$B_p(n)=\sum_{a=0}^{n-1}\sum_{t=1}^{p-1}\frac{1}{ap+t}.$$
وعند تثبيت الدقة \(p^S\)، يملك كل مقلوب التوسع \(p\)-adic
$$\frac{1}{ap+t}=\frac{1}{t}\cdot\frac{1}{1+ap/t}=\frac{1}{t}\sum_{j\ge 0}(-1)^j\left(\frac{ap}{t}\right)^j \pmod{p^S}.$$
وعند الجمع على البواقي الكاملة \(t=1,2,\dots,p-1\)، تختفي المساهمات الفردية، فلا يبقى إلا الأسس الزوجية. لذلك يمكن كتابة حد الكتلة modulo \(p^S\) على صورة تركيب خطي لقوى \(n\) الزوجية. أي توجد معاملات \(c_1,\dots,c_{(S-1)/2}\) بحيث
$$B_p(n)\equiv \sum_{j=1}^{(S-1)/2} c_j n^{2j} \pmod{p^S}.$$
لا تستخدم التطبيقات صيغة مغلقة لهذه المعاملات، بل تعيد بناءها من قيم عينية صغيرة لـ \(B_p(n)\) عبر حل نظام خطي صغير modulo \(p^S\).
افترض أن عقدة باقية \(n\) معروفة ومعها \(H_n\bmod p^r\)، وأنها تحقق أيضًا \(H_n\equiv 0\pmod p\). عندئذ تصبح الكمية \(H_n/p\) ذات معنى modulo \(p^{r-1}\). وباستخدام كثير الحدود الخاص بحد الكتلة نحصل على الابن الأول:
$$H_{pn}\equiv \frac{H_n}{p}+B_p(n) \pmod{p^{r-1}}.$$
ثم تتولد بقية الأبناء بالتحديثات التراكمية
$$H_{pn+k}=H_{pn+k-1}+\frac{1}{pn+k} \pmod{p^{r-1}}, \qquad 1\le k<p.$$
ولا يبقى إلى المستوى التالي إلا الأبناء الذين يحققون
$$H_{pn+k}\equiv 0\pmod p.$$
كل خطوة رفع تستهلك قوة واحدة من \(p\)، ولذلك فإن الدقة الابتدائية \(S\) تسمح بحد أقصى مقداره \(S-1\) من المستويات المفيدة.
إذا كان \(n\in J_p\)، فإن \(H_n/p\) يبقى عددًا صحيحًا \(p\)-adically، كما أن \(B_p(n)\) وكل ذيل قصير \(T_{p,k}(n)\) هما مجموع لوحدات \(p\)-adic. لذلك يبقى كامل المقطع
$$H_{pn},H_{pn+1},\dots,H_{pn+p-1}$$
صحيحًا \(p\)-adically. ومن ثم فإن أكبر فهرس في الكتلة المتولدة من \(n\) هو \(pn+(p-1)\). وعندما نعرف \(\max J_p\)، تصبح الإجابة النهائية
$$M(p)=p\,\max J_p+(p-1).$$
عند \(p=7\)، يكفي في البداية فحص \(1\le n<7\). نجد أن
$$H_6=1+\frac12+\frac13+\frac14+\frac15+\frac16=\frac{49}{20},$$
ومن ثم \(H_6\equiv 0\pmod 7\). وبالتالي
$$J_7^{(0)}=\{6\}.$$
بعد خطوة رفع واحدة لا يبقى إلا
$$J_7^{(1)}=\{42,48\}.$$
ثم تعطي الطبقات التالية
$$J_7^{(2)}=\{295,299,337,341\},$$
$$J_7^{(3)}=\{2096,2390\},$$
$$J_7^{(4)}=\{14675,16731,16735\},$$
$$J_7^{(5)}=\{102728\}.$$
بعد ذلك لا ينجو أي ابن آخر، ولهذا
$$\max J_7=102728.$$
وعليه
$$M(7)=7\cdot 102728+6=719102,$$
وهي بالضبط قيمة التحقق المستخدمة في التطبيقات.
تتبع تطبيقات C++ وPython وJava المسار الحسابي نفسه. أولًا تثبت الدقة عند \(S=8\) وتُحسب القوى \(p,p^2,\dots,p^S\) مسبقًا. بعد ذلك تُقيَّم \(B_p(n)\) لعدد قليل من العينات الصغيرة، ويُحل نظام خطي modulo \(p^S\) لاستخراج معاملات كثير الحدود الزوجي.
ثم تُحسب القيم الابتدائية \(H_1,H_2,\dots,H_{p-1}\) بترديد \(p^S\)، ولا يُحتفَظ إلا بالقيم الموافقة لـ \(0\pmod p\). هذه هي جبهة شجرة الرفع الأولى. لكل عقدة باقية، يقيّم التنفيذ التصحيح كثير الحدود، ويحصل على القيمة عند \(pn\)، ثم يمسح الأبناء الباقين بإضافة المقلوبات تباعًا، ولا يحتفظ إلا بالأبناء الذين يبقى مجموعهم التوافقي قابلًا للقسمة على \(p\). ويُحدَّث أكبر فهرس باقٍ طوال البحث.
نسخة C++ تضيف أيضًا توسيعًا متوازيًا عندما تصبح الجبهة كبيرة بما يكفي، لكن الحسابات واختبار البقاء متماثلان في اللغات الثلاث.
ليكن \(L_r\) عدد العقد الباقية في المستوى \(r\)، وليكن \(d=(S-1)/2\) عدد معاملات كثير الحدود. إعادة بناء كثير الحدود الخاص بحد الكتلة تكلف \(O(d^3)\) من العمليات المعيارية بعد حساب العينات. ومع تثبيت \(S=8\) تصبح هذه المرحلة ثابتة الحجم.
أما الكلفة الأساسية فتأتي من توسيع شجرة العقد الباقية. فكل عقدة تحتاج إلى تقييم كثير حدود بكلفة \(O(d)\)، ثم مسح \(p\) أبناء، وفي كل ابن توجد عملية معكوس معياري وعملية جمع معيارية. لذلك تكون الكلفة الكلية للبحث
$$O\left(\sum_r L_r(d+p)\right),$$
وهي، عند ثبات \(S\)، تكافئ عمليًا
$$O\left(p\sum_r L_r\right).$$
أما الذاكرة فخطية في حجم أكبر جبهة:
$$O\left(\max_r L_r\right).$$