Detailed mathematical approach

Problem Summary

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\).

Mathematical Approach

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.

Step 1: Define the Exceptional Set \(J_p\)

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.

Step 2: Decompose \(H_{pn+k}\) into a Core and a Tail

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.

Step 3: Represent the Block Term by an Even Polynomial

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\).

Step 4: Lift One Level and Lose One Power of \(p\)

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.

Step 5: Recover the Quantity \(M(p)\)

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).$

Worked Example: \(p=7\)

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.

How the Code Works

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.

Complexity Analysis

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).$

Footnotes and References

1. Problem page: https://projecteuler.net/problem=541
2. Harmonic numbers: Wikipedia - Harmonic number
3. p-adic numbers: Wikipedia - p-adic number
4. Hensel lifting background: Wikipedia - Hensel's lemma
5. Wolstenholme-type harmonic congruences: Wikipedia - Wolstenholme's theorem