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 1000 distinct solutions in positive integers when symmetric pairs are counted once, equivalently with \(x \ge y\).

A brute-force search over \(x\) and \(y\) is not the right viewpoint, because the interesting values of \(n\) are already large. The decisive observation is that the number of admissible pairs depends only on the prime factorization of \(n\), so the problem becomes an optimization over exponent vectors rather than a scan over denominator pairs.

Mathematical Approach

The implementations follow a short but very powerful chain of ideas: convert the reciprocal equation into a factorization of \(n^2\), express the solution count through the divisor function, and then search for the smallest prime factorization whose divisor count is large enough.

From reciprocals to a product equation

Start with

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

Multiplying by \(xyn\) gives

$n(x+y)=xy.$

Rearranging and completing the rectangle yields

$xy-nx-ny=0 \quad\Longrightarrow\quad (x-n)(y-n)=n^2.$

Thus every solution corresponds to positive integers \(a=x-n\) and \(b=y-n\) with \(ab=n^2\). Conversely, every factorization \(ab=n^2\) with \(a,b>0\) produces a solution \(x=n+a\), \(y=n+b\). Because the original equation is symmetric in \(x\) and \(y\), counting only solutions with \(x \ge y\) is exactly the same as counting each unordered factor pair of \(n^2\) once.

Counting solutions with the divisor function

Let \(d(m)\) denote the number of positive divisors of \(m\). The factor pairs of \(n^2\) normally come in mirrored pairs \((a,b)\) and \((b,a)\). Since \(n^2\) is a perfect square, there is one central unpaired divisor \(a=b=n\). Therefore

$\#\{\text{solutions with } x \ge y\}=\frac{d(n^2)+1}{2}.$

The problem asks for more than 1000 solutions, so

$\frac{d(n^2)+1}{2}>1000 \quad\Longleftrightarrow\quad d(n^2)>1999.$

Because \(d(n^2)\) is odd for every square, this is equivalent to

$d(n^2)\ge 2001.$

This odd threshold is what the implementations test during the search.

Prime factorization turns the condition into a product

If

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

then

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

and the divisor function gives

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

So the original Diophantine question becomes: find the smallest \(n=\prod p_i^{a_i}\) such that

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

At this point the problem is no longer about solving for \(x\) and \(y\) directly. It is about selecting a good exponent vector \((a_1,a_2,\dots)\).

Why the exponents can be assumed non-increasing

For a fixed multiset of exponents, the smallest possible value of \(n\) is obtained by assigning the largest exponent to the smallest prime, the next largest exponent to the next smallest prime, and so on. Indeed, if \(p<q\) but the exponents are arranged with \(a<b\), then

$p^a q^b > p^b q^a,$

so swapping those two exponents decreases \(n\) without changing \(d(n^2)\). Therefore an optimal factorization may always be written as

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

using the ascending primes \(2,3,5,7,\dots\). This monotonicity is the key invariant enforced by the depth-first search: once an exponent \(e\) is chosen for one prime, every later exponent is restricted to be at most \(e\).

Worked example: \(n=4\)

Take \(n=4=2^2\). Then

$d(n^2)=d(16)=2\cdot 2+1=5,$

so the formula predicts

$\frac{d(16)+1}{2}=\frac{5+1}{2}=3$

distinct solutions. The unordered factor pairs of \(16\) are

$16\cdot 1,\qquad 8\cdot 2,\qquad 4\cdot 4.$

They yield

$ (x,y)=(20,5),\qquad (12,6),\qquad (8,8), $

which are exactly the three solutions with \(x \ge y\). This small example captures the full divisor-pair mechanism used in the general argument.

A simple upper bound for the search

Using exponent 1 on the first seven primes already gives

$n=2\cdot 3\cdot 5\cdot 7\cdot 11\cdot 13\cdot 17=510510,$

and then

$d(n^2)=3^7=2187\ge 2001.$

So a valid solution definitely exists below \(510510\). This does not identify the optimum yet, but it shows that only a short prefix of small primes matters and that once the search finds its first admissible candidate, pruning will become effective very quickly.

How the Code Works

Search state

The C++, Python, and Java implementations all search over exponent vectors attached to the ascending primes \(2,3,5,7,\dots\). At each recursive call, the state contains four mathematical quantities: the next prime position, the largest exponent allowed for that prime, the current partial value of \(n\), and the current value of \(\prod (2a_i+1)\).

Recursive expansion of exponent vectors

For the current prime, the implementation tries exponents \(1,2,\dots\) up to the permitted maximum. Multiplying by the prime updates the candidate \(n\), and choosing exponent \(e\) multiplies the divisor-count product by \(2e+1\). The recursion then continues to the next prime with the new exponent \(e\) as the next upper bound, which automatically preserves the non-increasing pattern \(a_1 \ge a_2 \ge \cdots\).

If the divisor-count product reaches \(2001\) or more, the current \(n\) is already a complete candidate and can update the best answer. There is no benefit in appending more primes, because any extra positive exponent would only make \(n\) larger.

Pruning and termination

The crucial pruning rule is simple: if the current partial \(n\) is already at least as large as the best solution found so far, the entire branch is abandoned. Since later recursion only multiplies by additional primes, no descendant of that state can improve the answer.

Because the recursion updates both \(n\) and \(\prod (2a_i+1)\) incrementally, each visited state requires only constant-time arithmetic. All three language versions implement the same mathematical search, just in different syntactic forms.

Complexity Analysis

The natural complexity parameter here is the number \(T\) of exponent states visited by the depth-first search. Each state performs only a small amount of arithmetic, so the running time is \(O(T)\). There is no hidden scan over all \((x,y)\) pairs and no brute-force sweep over every integer up to the answer.

The memory usage is \(O(r)\), where \(r\) is the maximum recursion depth, i.e. the number of primes that receive a positive exponent on the current branch. For the threshold 1000, both \(T\) and \(r\) stay small because the non-increasing exponent invariant and the best-so-far pruning eliminate almost all of the nominal search tree. In practice the program finishes essentially instantly.

Footnotes and References

1. Project Euler Problem 108: https://projecteuler.net/problem=108
2. Diophantine equation: Wikipedia - Diophantine equation
3. Divisor function: Wikipedia - Divisor function
4. Fundamental theorem of arithmetic: Wikipedia - Fundamental theorem of arithmetic
5. Backtracking: Wikipedia - Backtracking