Detailed mathematical approach

Problem Summary

For this problem, the drone-delivery expectation reduces to evaluating

$E(n)=n\sum_{1\le i\le j\le n}\frac{1}{ij}=\frac n2\left(H_n^2+H_n^{(2)}\right),$

at the target value \(n=10^8\), where \(H_n=\sum_{k=1}^n \frac1k\) and \(H_n^{(2)}=\sum_{k=1}^n \frac1{k^2}\). The implementations exploit this exact identity and then switch to Euler-Maclaurin expansions, because summing \(10^8\) terms directly is unnecessary.

Mathematical Approach

Once the underlying probability model is simplified, everything depends only on harmonic sums. That replaces a large discrete process with a short analytic formula.

Step 1: Write the exact closed form

Introduce the ordinary harmonic number

$H_n=\sum_{k=1}^n\frac1k$

and the second-order harmonic number

$H_n^{(2)}=\sum_{k=1}^n\frac1{k^2}.$

The reduced expectation is then

$E(n)=\frac n2\left(H_n^2+H_n^{(2)}\right).$

This identity is exact, not asymptotic. For moderate \(n\), one may therefore compute \(E(n)\) directly from the defining sums.

Step 2: See why the double sum becomes harmonic numbers

Expanding the square gives

$H_n^2=\left(\sum_{i=1}^n\frac1i\right)\left(\sum_{j=1}^n\frac1j\right)=\sum_{i=1}^n\frac1{i^2}+2\sum_{1\le i<j\le n}\frac1{ij}.$

Adding \(H_n^{(2)}\) duplicates the diagonal terms, so

$H_n^2+H_n^{(2)}=2\sum_{1\le i\le j\le n}\frac1{ij}.$

Hence

$E(n)=n\sum_{1\le i\le j\le n}\frac1{ij}.$

This triangular-sum viewpoint explains why the closed form contains both \(H_n^2\) and \(H_n^{(2)}\): the off-diagonal pairs come from the square, while the diagonal \(i=j\) contributes the second-order harmonic sum.

Step 3: Use Euler-Maclaurin for large \(n\)

For the target input, direct summation is wasteful. The implementations use the standard asymptotic expansions

$H_n=\log n+\gamma+\frac{1}{2n}-\frac{1}{12n^2}+\frac{1}{120n^4}-\frac{1}{252n^6}+\frac{1}{240n^8}+O\left(n^{-10}\right),$

$H_n^{(2)}=\frac{\pi^2}{6}-\frac1n+\frac{1}{2n^2}-\frac{1}{6n^3}+\frac{1}{30n^5}-\frac{1}{42n^7}+\frac{1}{30n^9}+O\left(n^{-11}\right).$

Here \(\gamma\) is the Euler-Mascheroni constant, and \(\pi^2/6\) is the limit of the second-order harmonic sum.

Step 4: Substitute the expansions into the exact formula

Plugging those series into the identity for \(E(n)\) yields an approximation that is effectively exact for \(n=10^8\):

$E(n)\approx \frac n2\left[\left(\log n+\gamma+\frac{1}{2n}-\frac{1}{12n^2}+\frac{1}{120n^4}-\frac{1}{252n^6}+\frac{1}{240n^8}\right)^2+\left(\frac{\pi^2}{6}-\frac1n+\frac{1}{2n^2}-\frac{1}{6n^3}+\frac{1}{30n^5}-\frac{1}{42n^7}+\frac{1}{30n^9}\right)\right].$

The dominant growth is

$E(n)=\frac n2\left((\log n+\gamma)^2+\frac{\pi^2}{6}\right)+O(\log n).$

So the quantity grows on the order of \(n\log^2 n\), and the omitted tail of the asymptotic series is far too small to affect the nearest integer when \(n=10^8\).

Step 5: Why rounding to the nearest integer is safe

After multiplication by \(n/2\), the first neglected terms are still tiny: the harmonic-number remainder contributes on the order of \(n^{-9}\), and the second-order remainder contributes on the order of \(n^{-10}\). At \(n=10^8\), those analytic errors are astronomically below \(1/2\).

That is why the implementations can evaluate the truncated series with high precision and then round the result to the nearest integer.

Worked Example: \(n=5\)

For \(n=5\), direct summation is easy:

$H_5=1+\frac12+\frac13+\frac14+\frac15=\frac{137}{60},$

$H_5^{(2)}=1+\frac14+\frac19+\frac1{16}+\frac1{25}=\frac{5269}{3600}.$

Therefore

$E(5)=\frac52\left[\left(\frac{137}{60}\right)^2+\frac{5269}{3600}\right]=\frac{12019}{720}\approx 16.6930555556.$

This matches the exact checkpoint used by the implementations and confirms the closed form before any asymptotic approximation is applied.

How the Code Works

The C++, Python, and Java implementations all evaluate the same formula. They first prepare the constants \(\gamma\), \(\pi\), and the target input \(n=10^8\). Then they build the inverse powers \(n^{-1},n^{-2},\dots,n^{-9}\), because the Euler-Maclaurin expansions reuse them many times.

The C++ implementation also checks a few small values by direct summation, which verifies the exact harmonic identity numerically before the large-\(n\) path is used. For the target computation, the implementations form the truncated series for \(H_n\) and \(H_n^{(2)}\), substitute them into \(\frac n2(H_n^2+H_n^{(2)})\), and round to the nearest integer.

The Python and Java implementations use 50-digit decimal arithmetic. The Java version stores the natural logarithm of the target input as a high-precision constant, while the C++ version computes the logarithm directly in extended floating-point arithmetic.

Complexity Analysis

If one summed the defining series directly, the cost would be \(O(n)\) time and \(O(1)\) extra memory. The asymptotic method used for the real input evaluates only a fixed number of arithmetic operations on high-precision numbers, so with fixed precision it is \(O(1)\) time and \(O(1)\) memory. That is the point of replacing the huge sums with Euler-Maclaurin expansions.

Footnotes and References

1. Problem page: https://projecteuler.net/problem=724
2. Harmonic number: Wikipedia - Harmonic number
3. Generalized harmonic number: Wikipedia - Generalized harmonic number
4. Euler-Maclaurin formula: Wikipedia - Euler-Maclaurin formula
5. Euler's constant: Wikipedia - Euler's constant
6. Basel problem: Wikipedia - Basel problem