Detailed mathematical approach

Problem Summary

We seek the smallest positive integer \(n\) for which the Diophantine equation \( \frac{1}{x} + \frac{1}{y} = \frac{1}{n} \) has more than \(4{,}000{,}000\) distinct positive integer solutions. The successful approach does not search over pairs \((x,y)\) directly. Instead, it rewrites the equation as a divisor-count condition on \(n^2\), then searches over the prime factorization of \(n\).

That reformulation turns a two-variable equation into a multiplicative optimization problem. Once the exponent pattern of \(n\) is known, the number of solutions is determined, and the remaining task is to find the smallest \(n\) whose exponent pattern forces the solution count past the required threshold.

Mathematical Approach

The key step is to convert the reciprocal equation into a factorization problem.

From reciprocals to a factorization of \(n^2\)

Starting from

$\frac{1}{x} + \frac{1}{y} = \frac{1}{n},$

multiply through by \(xyn\):

$n(x+y)=xy.$

Rearranging gives

$xy-nx-ny=0,$

and after adding \(n^2\) to both sides we obtain

$\boxed{(x-n)(y-n)=n^2.}$

So every positive solution corresponds to a factorization \(uv=n^2\) with

$u=x-n,\qquad v=y-n,$

and conversely every positive factor pair of \(n^2\) produces a solution

$x=n+u,\qquad y=n+v.$

Counting distinct solutions with divisors of a square

If ordered pairs were counted, then each divisor \(u\mid n^2\) would determine exactly one pair \((u,n^2/u)\), so the ordered count would be \(d(n^2)\), where \(d\) is the divisor function.

Project Euler counts distinct solutions, so \((x,y)\) and \((y,x)\) are the same unless they lie on the diagonal. Because \(n^2\) is a square, there is exactly one central factorization with \(u=v=n\). Therefore the number of distinct solutions is

$\boxed{N(n)=\frac{d(n^2)+1}{2}.}$

Hence the requirement

$N(n) \gt 4{,}000{,}000$

is equivalent to

$d(n^2)\ge 8{,}000{,}001.$

The odd threshold is not accidental: every exponent in \(n^2\) is even, so \(d(n^2)\) is always an odd number.

Prime exponents and the divisor product

Write the prime factorization of \(n\) as

$n=\prod_{i=1}^{k} p_i^{a_i}.$

Then

$n^2=\prod_{i=1}^{k} p_i^{2a_i},$

and the divisor formula becomes

$d(n^2)=\prod_{i=1}^{k}(2a_i+1).$

So the original problem is equivalent to finding the smallest number of the form

$n=\prod p_i^{a_i}$

such that

$\prod (2a_i+1)\ge 8{,}000{,}001.$

At this point the variables \(x\) and \(y\) have disappeared completely; everything is controlled by the exponent vector \((a_1,a_2,\dots)\).

Why only consecutive primes with non-increasing exponents matter

To minimize \(n\), the exponents must be arranged in a very specific way.

First, the used primes must form an initial segment \(2,3,5,\dots\). If some larger prime \(q\) appears while a smaller prime \(p\lt q\) does not, replacing \(q\) by \(p\) leaves the divisor product unchanged and makes \(n\) smaller, so such a pattern cannot be optimal.

Second, the exponents may be assumed to be non-increasing:

$a_1\ge a_2\ge a_3\ge \cdots \ge 1.$

Indeed, if \(p\lt q\) but \(a_p \lt a_q\), then swapping those two exponents changes

$p^{a_p}q^{a_q}$

into

$p^{a_q}q^{a_p},$

which is smaller because the larger exponent is now attached to the smaller prime. The divisor product \(\prod(2a_i+1)\) is unchanged by the swap, so an optimal candidate must respect this monotonic ordering. This is the main invariant enforced by the recursive search.

Worked example: the smaller case \(n=6\)

A concrete example makes the divisor-count formula visible. Let \(n=6=2^1\cdot 3^1\). Then

$n^2=36=2^2\cdot 3^2,$

so

$d(n^2)=(2\cdot 1+1)(2\cdot 1+1)=9.$

Therefore the number of distinct solutions is

$\frac{9+1}{2}=5.$

The corresponding factor pairs of \(36\) are

$ (1,36),\ (2,18),\ (3,12),\ (4,9),\ (6,6), $

which produce

$ (x,y)=(7,42),\ (8,24),\ (9,18),\ (10,15),\ (12,12). $

This is exactly why the implementations track the divisor product \(\prod(2a_i+1)\): it counts solutions without ever constructing them explicitly.

The optimization actually solved by the search

After all the algebra, the real task is now clear: search through exponent vectors assigned to the consecutive primes

$2,3,5,7,\dots,53,$

maintain the current value of \(n\) and the current value of \(\prod(2a_i+1)\), and keep the smallest \(n\) whose divisor product reaches \(8{,}000{,}001\). The entire problem is therefore a depth-first search over a structured exponent space, not a search over solutions \((x,y)\).

How the Code Works

The C++, Python, and Java implementations use the same recursive state: the index of the next prime, the largest exponent that the next prime is allowed to take, the current candidate value of \(n\), the current divisor product \(\prod(2a_i+1)\), and the best complete answer found so far.

The first prime receives a generous upper bound, and every later prime is restricted by the exponent chosen for the previous one. If the current exponent is \(e\), the recursive call passes \(e\) down as the next upper bound. That enforces the invariant \(a_1\ge a_2\ge \cdots\) automatically.

At each level the implementation multiplies the current \(n\) by the current prime one more time, so trying exponents \(1,2,3,\dots\) is done incrementally rather than recomputing powers from scratch. In parallel, it multiplies the running divisor product by \(2e+1\). As soon as that product reaches \(8{,}000{,}001\), the current \(n\) is a valid answer candidate.

The pruning rule is simple and effective: if the current \(n\) has already reached or exceeded the best answer found so far, the whole branch is discarded, because any deeper continuation would only make \(n\) larger. The languages differ only in numeric representation: C++ uses 128-bit integer arithmetic, Python uses arbitrary-precision integers, and Java uses big integers, but the search logic is the same.

Complexity Analysis

The algorithm never enumerates the actual solution pairs \((x,y)\), so the meaningful cost measure is the number of recursive exponent states it visits. If that number is \(T\), then the running time is \(O(T)\), with only small arithmetic work per visited state.

The recursion depth is bounded by the number of primes supplied to the search, which here is 16. Thus the auxiliary space is \(O(16)\), or more generally \(O(k)\) for \(k\) explored primes. In practice the search is fast because the non-increasing exponent constraint removes most combinations, and the current-best cutoff prunes the remaining branches aggressively.

Footnotes and References

1. Problem page: https://projecteuler.net/problem=110
2. Divisor function: Wikipedia - Divisor function
3. Fundamental theorem of arithmetic: Wikipedia - Fundamental theorem of arithmetic
4. Diophantine equation: Wikipedia - Diophantine equation
5. Branch and bound: Wikipedia - Branch and bound