Detailed mathematical approach

Problem Summary

We must count the number \(F(L)\) of positive integer triples \((x,y,n)\) satisfying

$x \lt y \le L,\qquad \frac1x+\frac1y=\frac1n.$

The symmetry between \(x\) and \(y\) means that the strict inequality \(x \lt y\) removes the mirror-image duplicate of every solution. The published checkpoints are \(F(15)=4\) and \(F(1000)=1069\), while the real target is \(F(10^{12})\).

Mathematical Approach

Step 1: Extract the gcd and obtain a coprime parametrization

Let

$g=\gcd(x,y),\qquad x=g\,u,\qquad y=g\,v,$

with \(u \lt v\) and \(\gcd(u,v)=1\). Substituting into the reciprocal equation gives

$n(u+v)=g\,u\,v.$

Because \(\gcd(u,u+v)=1\) and \(\gcd(v,u+v)=1\), we can read divisibility information directly from this identity. Since \(u\mid n(u+v)\) and \(u\) is coprime to \(u+v\), it follows that \(u\mid n\). The same argument gives \(v\mid n\). As \(\gcd(u,v)=1\), we conclude that \(uv\mid n\), so we may write

$n=tuv,\qquad t\in\mathbb Z_{>0}.$

Then the previous equation simplifies to \(g=t(u+v)\). Therefore every admissible solution has the form

$x=t\,u(u+v),\qquad y=t\,v(u+v),\qquad n=t\,u\,v,$

with \(u,v,t \in \mathbb Z_{>0}\), \(u \lt v\), and \(\gcd(u,v)=1\). Conversely, every triple of this form satisfies \(\frac1x+\frac1y=\frac1n\), so this is a bijection and not just a one-way construction.

Step 2: Rewrite the search space in the form used by the algorithm

Introduce the new variables

$a=v,\qquad s=u+v.$

Then \(u=s-a\). The condition \(1\le u \lt v\) becomes

$a+1\le s\le 2a-1.$

The coprimality condition is preserved because

$\gcd(a,s)=\gcd(v,u+v)=\gcd(u,v)=1.$

The upper bound on \(y\) is now

$y=t\,a\,s\le L.$

So for a fixed coprime pair \((a,s)\), the number of admissible scale factors \(t\) is exactly

$1\le t\le \left\lfloor\frac{L}{as}\right\rfloor,$

which contributes \(\left\lfloor L/(as)\right\rfloor\) solutions. Hence

$F(L)=\sum_{a=1}^{A_{\max}}\ \sum_{s=a+1}^{U(a)} \mathbf{1}_{\gcd(a,s)=1}\left\lfloor\frac{L}{as}\right\rfloor,$

where

$U(a)=\min\left(2a-1,\left\lfloor\frac{L}{a}\right\rfloor\right).$

The outer loop can stop once no valid \(s\) exists. Since the smallest allowed \(s\) is \(a+1\), we need \(a(a+1)\le L\). Therefore

$A_{\max}=\left\lfloor\frac{\sqrt{1+4L}-1}{2}\right\rfloor.$

Step 3: Count coprime values in an interval with Möbius inversion

For a fixed outer parameter \(a\), define

$C_a(\alpha,\beta)=\#\left\{s\in\mathbb Z:\alpha\le s\le \beta,\ \gcd(a,s)=1\right\}.$

The standard Möbius identity is

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

Summing this over an interval gives

$C_a(\alpha,\beta)=\sum_{d\mid a}\mu(d)\left(\left\lfloor\frac{\beta}{d}\right\rfloor-\left\lfloor\frac{\alpha-1}{d}\right\rfloor\right).$

Only squarefree divisors contribute, because \(\mu(d)=0\) whenever \(d\) contains a repeated prime factor. The implementation therefore precomputes, for each outer value, exactly those divisors with nonzero Möbius sign and reuses them for every interval query.

Step 4: Group equal quotients into harmonic blocks

Fix \(a\) and set

$M=\left\lfloor\frac{L}{a}\right\rfloor,\qquad q(s)=\left\lfloor\frac{M}{s}\right\rfloor.$

As \(s\) increases, \(q(s)\) remains constant across contiguous ranges. If the current block starts at \(\ell\) and

$q=\left\lfloor\frac{M}{\ell}\right\rfloor,$

then the same quotient remains valid up to

$r=\min\left(U(a),\left\lfloor\frac{M}{q}\right\rfloor\right).$

Every \(s\in[\ell,r]\) contributes the same multiplier \(q\), so the whole block contributes

$q\cdot C_a(\ell,r).$

This is the main speedup: instead of scanning each inner value one by one, the algorithm jumps directly between maximal ranges on which the floor quotient is constant.

Worked Example: \(L=15\)

Here

$A_{\max}=\left\lfloor\frac{\sqrt{61}-1}{2}\right\rfloor=3.$

For \(a=1\), the interval \(a+1\le s\le 2a-1\) is empty, so there is no contribution.

For \(a=2\), the only possible value is \(s=3\). It is coprime to \(2\), and it contributes

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

For \(a=3\), the admissible values are \(s=4\) and \(s=5\). Both are coprime to \(3\), giving

$\left\lfloor\frac{15}{3\cdot 4}\right\rfloor=1,\qquad \left\lfloor\frac{15}{3\cdot 5}\right\rfloor=1.$

Therefore

$F(15)=2+1+1=4,$

which matches the checkpoint exactly.

How the Code Works

The C++, Python, and Java implementations all follow this derivation. They first compute \(A_{\max}\) using an integer square root so the outer range is exact. Then they build the Möbius function up to that limit with a sieve, and from it they assemble, for every outer value, the signed list of squarefree divisors needed by the interval formula above.

During the main summation, the implementation determines the valid inner interval, advances through maximal quotient blocks, and for each block evaluates the coprime count using the precomputed signed divisors. The C++ and Java versions store this data in compact flattened form to reduce allocation overhead, while the Python version keeps direct per-value lists. The Java implementation also parallelizes independent outer ranges across available CPU cores.

Complexity Analysis

Let \(A=A_{\max}=O(\sqrt L)\). Computing the Möbius function is sieve-like and near-linear in \(A\). Building the stored squarefree-divisor data costs

$\sum_{\substack{d\le A\\ \mu(d)\ne 0}}\left\lfloor\frac{A}{d}\right\rfloor=O(A\log A),$

and the memory needed for that data has the same order.

The summation stage is much smaller than a naive search over all \(x\) and \(y\). It processes quotient blocks rather than individual inner values, and each block performs one Möbius interval count over the stored divisors of the current outer value. In practice this behaves close to \(O(\sqrt L\log L)\), which is why the method remains feasible for \(L=10^{12}\), whereas direct enumeration would be hopelessly large.

Footnotes and References

1. Problem page: https://projecteuler.net/problem=454
2. Möbius function: Wikipedia — Möbius function
3. Möbius inversion formula: Wikipedia — Möbius inversion formula
4. Unit-fraction background: Wikipedia — Egyptian fraction