Detailed mathematical approach

Problem Summary

For each integer \(M \ge 2\), define

$R(M)=\sum_{\substack{1\le p \lt q \le M\\ p+q \ge M\\ \gcd(p,q)=1}}\frac{1}{pq}, \qquad S(N)=\sum_{M=2}^{N}R(M).$

The checkpoints are \(S(2)=\tfrac12\), \(S(10)=6.9146825397\ldots\), and \(S(100)=58.2962380621\ldots\). The task is to evaluate \(S(10^7)\) to four decimal places. A direct evaluation over all triples \((M,p,q)\) is far too slow, so the implementation reorganizes the sum and replaces the coprimality test by Möbius inclusion-exclusion.

Mathematical Approach

Step 1: Exchange the Order of Summation

Fix a coprime pair \(1 \le p \lt q \le N\). Such a pair contributes to \(R(M)\) exactly when \(q \le M\) and \(M \le p+q\). Since \(M\) is also bounded by \(N\), the admissible values are

$q \le M \le \min(N,p+q).$

Therefore this pair appears with multiplicity

$w_N(p,q)=\min(N,p+q)-q+1.$

After switching the order of summation,

$S(N)=\sum_{\substack{1\le p \lt q \le N\\ \gcd(p,q)=1}}\frac{w_N(p,q)}{pq} =\sum_{q=2}^{N}\frac{1}{q}\sum_{\substack{1\le p \lt q\\ \gcd(p,q)=1}}\frac{w_N(p,q)}{p}.$

This is the key simplification: the problem becomes a sum of independent contributions indexed by the larger denominator \(q\).

Step 2: Encode Coprimality with Möbius Inversion

Define the coprime reciprocal sum and the coprime counting function

$A_q(m)=\sum_{\substack{1\le p \le m\\ \gcd(p,q)=1}}\frac{1}{p}, \qquad C_q(m)=\sum_{\substack{1\le p \le m\\ \gcd(p,q)=1}}1.$

Using the identity

$\mathbf{1}_{\gcd(p,q)=1}=\sum_{d \mid \gcd(p,q)}\mu(d),$

we can rewrite both quantities in terms of divisors of \(q\). If \(H_t=\sum_{j=1}^{t}\frac{1}{j}\) denotes the \(t\)-th harmonic number, then

$A_q(m)=\sum_{d \mid q}\frac{\mu(d)}{d}H_{\lfloor m/d \rfloor}, \qquad C_q(m)=\sum_{d \mid q}\mu(d)\left\lfloor \frac{m}{d}\right\rfloor.$

Only squarefree divisors matter, because \(\mu(d)=0\) whenever \(d\) contains a repeated prime factor. That is why the implementation factors \(q\) into distinct primes and enumerates all squarefree divisor subsets with their Möbius signs.

Step 3: The Easy Half \(q \le \lfloor N/2 \rfloor\)

If \(q \le \lfloor N/2 \rfloor\), then for every \(p \lt q\) we have

$p+q \le (q-1)+q = 2q-1 \le N-1,$

so the upper bound \(\min(N,p+q)\) is just \(p+q\). Hence

$w_N(p,q)=p+1.$

The contribution of this \(q\) becomes

$\frac{1}{q}\sum_{\substack{1\le p \lt q\\ \gcd(p,q)=1}}\frac{p+1}{p} =\frac{1}{q}\left(\sum_{\substack{1\le p \lt q\\ \gcd(p,q)=1}}1+\sum_{\substack{1\le p \lt q\\ \gcd(p,q)=1}}\frac{1}{p}\right).$

The first sum is \(\varphi(q)\), because among \(1,2,\dots,q-1\) there are exactly \(\varphi(q)\) integers coprime to \(q\). Therefore

$T_N(q)=\frac{\varphi(q)+A_q(q-1)}{q}\qquad \text{for } q \le \left\lfloor \frac{N}{2}\right\rfloor.$

Step 4: The Truncated Half \(q > \lfloor N/2 \rfloor\)

Now let \(q > \lfloor N/2 \rfloor\) and set \(a=N-q\). Then \(0 \le a \lt q\), and the weight changes at \(p=a\):

$w_N(p,q)= \begin{cases} p+1, & 1 \le p \le a,\\ a+1, & a \lt p \lt q. \end{cases}$

So the contribution splits into two pieces:

$T_N(q)=\frac{1}{q}\left(\sum_{\substack{1\le p \le a\\ \gcd(p,q)=1}}\left(1+\frac{1}{p}\right) +(a+1)\sum_{\substack{a \lt p \lt q\\ \gcd(p,q)=1}}\frac{1}{p}\right).$

Writing the second reciprocal sum as \(A_q(q-1)-A_q(a)\), we obtain

$T_N(q)=\frac{C_q(a)+A_q(a)+(a+1)\bigl(A_q(q-1)-A_q(a)\bigr)}{q}.$

This is the exact branch used by the implementation once \(q\) passes the midpoint.

Step 5: Final Closed Form

Combining the two ranges gives the complete formula

$\boxed{ S(N)= \sum_{q=2}^{\lfloor N/2 \rfloor}\frac{\varphi(q)+A_q(q-1)}{q} + \sum_{q=\lfloor N/2 \rfloor+1}^{N} \frac{C_q(N-q)+A_q(N-q)+(N-q+1)\bigl(A_q(q-1)-A_q(N-q)\bigr)}{q} .}$

All expensive work has now been reduced to three precomputable ingredients: Euler's totient values, harmonic prefixes, and factorizations of the integers \(2,\dots,N\).

Worked Example: \(N=10\)

For \(q=5\), we are still in the first branch. Since \(\varphi(5)=4\) and

$A_5(4)=1+\frac12+\frac13+\frac14=\frac{25}{12},$

we get

$T_{10}(5)=\frac{4+\frac{25}{12}}{5}=\frac{73}{60}=1.216666\ldots$

For \(q=8\), we are in the second branch with \(a=10-8=2\). Here

$C_8(2)=1,\qquad A_8(2)=1,\qquad A_8(7)-A_8(2)=\frac13+\frac15+\frac17=\frac{71}{105},$

so

$T_{10}(8)=\frac{1+1+3\cdot \frac{71}{105}}{8}=\frac{141}{280}=0.503571\ldots$

Summing all contributions from \(q=2\) to \(10\) yields

$S(10)=\frac{3485}{504}=6.914682539682\ldots,$

which matches the published checkpoint.

How the Code Works

The C++, Python, and Java implementations all follow the same pipeline. First they build a smallest-prime-factor table up to \(N\), which makes factorization of each \(q\) fast and also supports a linear-time computation of \(\varphi(q)\). Next they precompute the harmonic prefix array \(H_0,H_1,\dots,H_N\). Then each value of \(q\) is processed independently: factor \(q\), enumerate the squarefree divisors induced by its distinct prime factors, evaluate the Möbius formulas for \(A_q\) and \(C_q\), apply the appropriate branch, and add the result to the global sum.

The C++ and Java implementations additionally use compensated summation to reduce floating-point drift, and the C++ version can split the outer loop into parallel chunks. None of those engineering choices changes the mathematics; they only improve numerical stability and runtime for the target size.

Complexity Analysis

Building the smallest-prime-factor table, the totient table, and the harmonic prefix array costs \(O(N)\) time and \(O(N)\) memory. For each \(q\), the factorization is recovered from the precomputed table, and the Möbius expansion enumerates exactly \(2^{\omega(q)}\) squarefree divisors, where \(\omega(q)\) is the number of distinct prime factors of \(q\). Therefore the total running time is

$O\!\left(N+\sum_{q=2}^{N}2^{\omega(q)}\right),$

with \(O(N)\) memory usage. For the specific target \(N=10^7\), every integer has at most eight distinct prime factors, so the subset enumeration remains very small in practice. Parallel execution reduces wall-clock time but does not change the asymptotic work.

Footnotes and References

1. Problem page: https://projecteuler.net/problem=441
2. Möbius function and inversion: Wikipedia — Möbius function
3. Harmonic numbers: Wikipedia — Harmonic number
4. Euler's totient function: Wikipedia — Euler's totient function
5. Compensated summation: Wikipedia — Kahan summation algorithm