Detailed mathematical approach

Problem Summary

The geometry of the incenter and circumcenter can be reduced to a primitive arithmetic parameterization. After that reduction, the quantity computed by the implementations is

$F(L)=\sum_{\substack{1\le r<s<2r\\ \gcd(r,s)=1\\ t\ge 1,\ rst\le L}} rst.$

Here \((r,s)\) describes a primitive configuration and \(t\) is the scaling factor. A direct scan over all admissible triples \((r,s,t)\) is too slow for \(L=10^9\), so the solution turns the sum into divisor-based range queries and batches equal floor values.

Mathematical Approach

The arithmetic form above is the key step. Once the geometry has been compressed into coprime integer parameters, the rest of the problem becomes a careful summation problem in elementary number theory.

Step 1: Describe the admissible domain

The reduced parameterization gives the conditions

$1\le r,\qquad r+1\le s\le 2r-1,\qquad \gcd(r,s)=1,\qquad rst\le L.$

For fixed \(r\) and \(s\), the scaling factor \(t\) can range from \(1\) up to \(\left\lfloor L/(rs)\right\rfloor\). Also, if \(r>\sqrt{L}\), then even the smallest allowed value \(s=r+1\) gives \(rs>r^2>L\), so no triple can contribute. Therefore the outer loop only needs

$1\le r\le R=\left\lfloor\sqrt{L}\right\rfloor.$

For each such \(r\), the valid interval for \(s\) is

$r+1\le s\le U_r=\min\!\left(2r-1,\left\lfloor\frac{L}{r}\right\rfloor\right).$

Step 2: Remove the scaling loop with triangular numbers

For one fixed coprime pair \((r,s)\), the total contribution of all admissible scales is

$\sum_{t=1}^{\lfloor L/(rs)\rfloor} rst=rs\sum_{t=1}^{\lfloor L/(rs)\rfloor} t.$

Introduce the triangular-number function

$T(n)=\frac{n(n+1)}{2}.$

Then the whole problem becomes

$F(L)=\sum_{r=1}^{R} r\sum_{\substack{r+1\le s\le U_r\\ \gcd(r,s)=1}} s\,T\!\left(\left\lfloor\frac{L}{rs}\right\rfloor\right).$

This is already a major simplification: the implementations never iterate over \(t\) explicitly.

Step 3: Encode coprimality by Möbius inversion

For a fixed \(r\), define the weighted coprime prefix sum

$C_r(x)=\sum_{\substack{1\le s\le x\\ \gcd(r,s)=1}} s.$

Using the standard identity

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

we obtain

$\begin{aligned} C_r(x) &=\sum_{1\le s\le x} s\sum_{d\mid \gcd(r,s)}\mu(d)\\ &=\sum_{d\mid r}\mu(d)\sum_{\substack{1\le s\le x\\ d\mid s}} s\\ &=\sum_{d\mid r}\mu(d)\,d\,T\!\left(\left\lfloor\frac{x}{d}\right\rfloor\right). \end{aligned}$

Only squarefree divisors matter, because \(\mu(d)=0\) for every non-squarefree \(d\). This is why the implementations precompute, for each \(r\), the squarefree divisors together with their Möbius signs.

Step 4: Turn interval sums into prefix differences

Once \(C_r(x)\) is available, any weighted coprime range sum follows immediately:

$\sum_{\substack{a\le s\le b\\ \gcd(r,s)=1}} s=C_r(b)-C_r(a-1).$

So the remaining challenge is not coprimality anymore. It is the changing floor value inside

$T\!\left(\left\lfloor\frac{L}{rs}\right\rfloor\right).$

Step 5: Batch all equal floor values

Fix \(r\) and suppose the current left endpoint in the \(s\)-interval is \(a\). Set

$m=\left\lfloor\frac{L}{ra}\right\rfloor.$

As \(s\) increases, the quantity \(\left\lfloor L/(rs)\right\rfloor\) stays equal to \(m\) on the whole block

$a\le s\le b=\min\!\left(U_r,\left\lfloor\frac{L}{rm}\right\rfloor\right).$

Therefore that entire block contributes

$r\,T(m)\left(C_r(b)-C_r(a-1)\right).$

Then the next block starts at \(b+1\). This is the same floor-grouping principle used in hyperbola-style divisor summation: one maximal interval replaces many identical floor evaluations.

Worked Example: \(L=15\)

Here

$R=\left\lfloor\sqrt{15}\right\rfloor=3.$

For \(r=1\), the interval \(r+1\le s\le U_r\) is empty, so there is no contribution.

For \(r=2\), the only admissible value is \(s=3\), and \(\gcd(2,3)=1\). Then

$\left\lfloor\frac{15}{2\cdot 3}\right\rfloor=2,$

so the contribution is

$2\cdot 3\cdot T(2)=6\cdot 3=18.$

For \(r=3\), the admissible values are \(s=4\) and \(s=5\). Both are coprime to \(3\), and both satisfy

$\left\lfloor\frac{15}{3s}\right\rfloor=1.$

Hence their contributions are

$3\cdot 4\cdot T(1)=12,\qquad 3\cdot 5\cdot T(1)=15.$

Adding everything gives

$18+12+15=45,$

which matches the checkpoint used by the implementations.

How the Code Works

The C++, Python, and Java implementations all follow the same structure. They first compute \(R=\lfloor\sqrt{L}\rfloor\), because no larger \(r\) can appear. Next they build a smallest-prime-factor sieve up to \(R\), which lets them factor every \(r\) quickly.

From those factorizations they generate, for each \(r\), the list of squarefree divisors and the corresponding Möbius sign. A coprime prefix query \(C_r(x)\) is then evaluated by summing \(d\,T(\lfloor x/d\rfloor)\) over that divisor list with the correct sign.

The main loop runs over \(r\). For each \(r\), it scans the valid \(s\)-interval in maximal blocks on which \(\left\lfloor L/(rs)\right\rfloor\) is constant. On each block, the implementation asks for one coprime interval sum, multiplies it by \(r\) and the corresponding triangular number, and adds the result to the total.

The numeric types are chosen to avoid overflow: the Python version uses arbitrary-precision integers automatically, the Java version uses big integers for the accumulated total, and the C++ version uses 128-bit arithmetic for the same reason.

Complexity Analysis

Let \(R=\lfloor\sqrt{L}\rfloor\). The smallest-prime-factor sieve costs \(O(R\log\log R)\) time and \(O(R)\) memory. Building the squarefree-divisor tables costs additional work proportional to the total number of stored squarefree divisors, namely

$O\!\left(\sum_{r\le R} 2^{\omega(r)}\right),$

where \(\omega(r)\) is the number of distinct prime factors of \(r\).

During the main summation, each floor-constant block for a given \(r\) needs two coprime prefix evaluations, and each such evaluation iterates over the same squarefree-divisor list. If \(B(r)\) denotes the number of blocks for that \(r\), then the main-loop cost is

$O\!\left(\sum_{r\le R} B(r)\,2^{\omega(r)}\right).$

This is far smaller in practice than iterating over every admissible triple \((r,s,t)\), because the \(t\)-sum has been collapsed into a triangular number and the \(s\)-loop is processed in batches rather than one value at a time.

Footnotes and References

1. Problem page: https://projecteuler.net/problem=496
2. Möbius function: Wikipedia — Möbius function
3. Möbius inversion formula: Wikipedia — Möbius inversion formula
4. Triangular number: Wikipedia — Triangular number
5. Dirichlet hyperbola method: Wikipedia — Dirichlet hyperbola method