Detailed mathematical approach

Problem Summary

The ellipse geometry in Problem 471 reduces to the arithmetic sum

$G(n)=\sum_{a=3}^{n}\sum_{b=1}^{\lfloor (a-1)/2 \rfloor} r(a,b), \qquad r(a,b)=\frac{b(a-2b)}{a-b}.$

Each admissible pair \((a,b)\) contributes one radius term. The challenge is that the real input is enormous, so a direct double loop is far too slow. The implementations therefore reorganize the sum into two closed polynomial pieces and one harmonic interval.

The small checkpoints used by the implementations are \(r(3,1)=1/2\), \(r(6,2)=1\), \(r(12,3)=2\), \(G(10)=2.059722222\times 10^1\), and \(G(100)=1.922360980\times 10^4\).

Mathematical Approach

Start from

$G(n)=\sum_{a=3}^{n}\sum_{b=1}^{\lfloor (a-1)/2 \rfloor}\frac{b(a-2b)}{a-b}.$

The decisive simplification is to sum by the denominator \(a-b\) instead of summing by \(a\) first.

Step 1: Reparameterize the Summation Region

Set

$c=a-b,$

so \(a=b+c\). Then the summand becomes

$r(a,b)=\frac{b(c-b)}{c}=b-\frac{b^2}{c}.$

The original bounds are

$b\ge 1,\qquad b\le \left\lfloor\frac{a-1}{2}\right\rfloor,\qquad a\le n.$

After substituting \(a=b+c\), they become

$1\le b\le c-1,\qquad 1\le b\le n-c.$

Hence for fixed \(c\) the inner index runs over

$1\le b\le B(c),\qquad B(c)=\min(c-1,n-c),$

and the whole sum rewrites as

$G(n)=\sum_{c=2}^{n-1}\sum_{b=1}^{B(c)}\left(b-\frac{b^2}{c}\right).$

Step 2: Split at \(c=\lfloor n/2 \rfloor\)

Define

$q=\left\lfloor\frac{n}{2}\right\rfloor,\qquad m=\left\lfloor\frac{n-1}{2}\right\rfloor.$

If \(2\le c\le q\), then \(c-1\le n-c\), so \(B(c)=c-1\).

If \(q\lt c\le n-1\), then \(n-c\lt c\), so \(B(c)=n-c\).

Therefore

$G(n)=\sum_{c=2}^{q}\sum_{b=1}^{c-1}\left(b-\frac{b^2}{c}\right)+\sum_{c=q+1}^{n-1}\sum_{b=1}^{n-c}\left(b-\frac{b^2}{c}\right).$

The low-denominator block becomes a pure polynomial. The high-denominator block is also mostly polynomial, except for one harmonic contribution.

Step 3: Evaluate the Low-Denominator Block

For \(2\le c\le q\),

$\sum_{b=1}^{c-1} b=\frac{(c-1)c}{2},\qquad \sum_{b=1}^{c-1} b^2=\frac{(c-1)c(2c-1)}{6}.$

So the inner sum is

$\sum_{b=1}^{c-1}\left(b-\frac{b^2}{c}\right)=\frac{(c-1)c}{2}-\frac{(c-1)c(2c-1)}{6c}=\frac{c^2-1}{6}.$

Summing over \(c\) gives

$G_{\text{low}}(n)=\sum_{c=2}^{q}\frac{c^2-1}{6}=\frac{q(2q^2+3q-5)}{36}.$

This is the first closed polynomial used by the implementation.

Step 4: Evaluate the High-Denominator Block

For \(q\lt c\le n-1\), the upper limit is \(n-c\). Thus

$\sum_{b=1}^{n-c}\left(b-\frac{b^2}{c}\right)=\frac{(n-c)(n-c+1)}{2}-\frac{(n-c)(n-c+1)(2n-2c+1)}{6c}.$

The first term is polynomial. Writing \(x=n-c\), where \(x=1,2,\dots,m\), gives

$\sum_{c=n-m}^{n-1}\frac{(n-c)(n-c+1)}{2}=\sum_{x=1}^{m}\frac{x(x+1)}{2}=\frac{m(m+1)(m+2)}{6}.$

The second term is the only non-polynomial piece. Expand its numerator over the denominator \(c\):

$\frac{(n-c)(n-c+1)(2n-2c+1)}{c}=\frac{n(n+1)(2n+1)}{c}-(6n^2+6n+1)+(6n+3)c-2c^2.$

Now set

$L=n-m,\qquad R=n-1,\qquad H_{L,R}=\sum_{k=L}^{R}\frac{1}{k}.$

Also define the elementary sums

$S_1=\sum_{k=L}^{R}k=\frac{(L+R)m}{2},$

$S_2=\sum_{k=L}^{R}k^2=\frac{R(R+1)(2R+1)-(L-1)L(2L-1)}{6}.$

Then the rational tail becomes

$T(n)=n(n+1)(2n+1)H_{L,R}-(6n^2+6n+1)m+(6n+3)S_1-2S_2.$

Step 5: Assemble the Closed Form

Combining the low block, the polynomial part of the high block, and the rational tail gives

$\boxed{G(n)=\frac{q(2q^2+3q-5)}{36}+\frac{m(m+1)(m+2)}{6}-\frac{T(n)}{6}.}$

So the entire problem has been reduced to evaluating a single harmonic interval. For short intervals the implementations sum reciprocals directly. For large intervals they use the Euler-Maclaurin expansion

$H_t=\log t+\gamma+\frac{1}{2t}-\frac{1}{12t^2}+\frac{1}{120t^4}-\frac{1}{252t^6}+\frac{1}{240t^8}-\frac{5}{660t^{10}}+O(t^{-12}),$

and then compute \(H_{L,R}=H_R-H_{L-1}\).

Worked Example: \(n=10\)

Here

$q=\left\lfloor\frac{10}{2}\right\rfloor=5,\qquad m=\left\lfloor\frac{9}{2}\right\rfloor=4,\qquad L=6,\qquad R=9.$

The low block is

$G_{\text{low}}(10)=\frac{5(2\cdot 5^2+3\cdot 5-5)}{36}=\frac{25}{3}.$

The polynomial part of the high block is

$\frac{4\cdot 5\cdot 6}{6}=20.$

The harmonic interval is

$H_{6,9}=\frac16+\frac17+\frac18+\frac19=\frac{275}{504}.$

Also

$S_1=6+7+8+9=30,\qquad S_2=6^2+7^2+8^2+9^2=230.$

Therefore

$T(10)=10\cdot 11\cdot 21\cdot \frac{275}{504}-661\cdot 4+63\cdot 30-2\cdot 230=\frac{557}{12}.$

So

$G(10)=\frac{25}{3}+20-\frac{1}{6}\cdot\frac{557}{12}=\frac{1483}{72}=20.597222\ldots,$

which matches the checkpoint \(2.059722222\times 10^1\).

How the Code Works

The C++, Python, and Java implementations all evaluate the same closed form. They compute \(q\), \(m\), \(L\), and \(R\), evaluate the two polynomial pieces, and then evaluate the harmonic interval \(H_{L,R}\).

When the harmonic interval is small, the implementation sums reciprocals directly. When the interval is large, it evaluates the truncated Euler-Maclaurin series for \(H_R\) and \(H_{L-1}\) and subtracts them. The C++ implementation also keeps a direct double-sum routine for checkpoint validation on small inputs.

Finally, the implementation forms \(T(n)\), combines the three contributions, and prints the result in scientific notation with one digit before the decimal point and nine digits after it.

Complexity Analysis

The fast method uses \(O(1)\) memory. For very large inputs its running time is effectively \(O(1)\), because the only non-polynomial piece is replaced by a constant-length asymptotic expansion. For short harmonic intervals the code may still do a direct reciprocal sum up to a fixed cutoff, so the practical fast path is bounded independently of \(n\). The validation path that literally enumerates admissible pairs remains \(O(n^2)\) time and \(O(1)\) extra memory.

Footnotes and References

1. Problem page: https://projecteuler.net/problem=471
2. Harmonic number: Wikipedia — Harmonic number
3. Euler-Maclaurin formula: Wikipedia — Euler-Maclaurin formula
4. Faulhaber's formula: Wikipedia — Faulhaber's formula