Problem 108: Diophantine Reciprocals I

Project Euler Problem 108 Solution

EulerSolve provides an optimized solution for Project Euler Problem 108, Diophantine Reciprocals I, with C++, Python, Java, and a step-by-step mathematical explanation.

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$....

Detailed mathematical approach

Problem Summary

We look for positive integers $x$, $y$ and $n$ with

$$\frac1x+\frac1y=\frac1n.$$

For $n=4$ there are exactly three distinct solutions: $\frac15+\frac1{20}$, $\frac16+\frac1{12}$ and $\frac18+\frac18$. The example shows what "distinct" means: $\frac1{20}+\frac15$ is not listed separately, so solutions that differ only in the order of $x$ and $y$ count once, and we may assume $x\le y$. The task is to find the least $n$ for which the number of distinct solutions exceeds one thousand. ^[1] The statement calls the problem an easier version of Problem 110, which asks the same question for more than four million solutions. ^[2]

The whole solution in one sentence: subtracting $n$ from $x$ and from $y$ turns every solution into a way of splitting $n^2$ into two factors, so the number of solutions is $(d(n^2)+1)/2$, where $d(m)$ is the number of divisors of $m$; this number depends only on the exponents in the prime factorisation of $n$, and a short search over exponent patterns on the primes $2,3,5,7,\dots$ finds the least $n$ whose square has at least $2001$ divisors.

Mathematical Approach

We solve small cases by hand, find and prove a formula for the number of solutions, and then search for the least $n$ with enough of them. Three methods appear on the way, each much faster than the one before.

1. Small cases by hand

For $n=1$ the only solution is $\frac12+\frac12=1$. For $n=2$ there are two, $\frac13+\frac16$ and $\frac14+\frac14$, and for $n=3$ also two, $\frac14+\frac1{12}$ and $\frac16+\frac16$. For $n=6$ there are five; as pairs $(x,y)$ they are $(7,42)$, $(8,24)$, $(9,18)$, $(10,15)$ and $(12,12)$. Section 2 shows that $x$ only has to be tried from $n+1$ to $2n$, so each list is a short check. The numbers of solutions for $n=1,2,\dots,12$ are

$1, 2, 2, 3, 2, 5, 2, 4, 3, 5, 2, 8.$

The On-Line Encyclopedia of Integer Sequences lists these counts as the number of ways to write $1/n$ as a sum of two unit fractions. ^[3] Every prime $n$ has exactly two solutions, while numbers with many small divisors, such as $6$ and $12$, have many. The values of $n$ with more solutions than every smaller number are $1, 2, 4, 6, 12, 24,$ $30, 60, 120, 180, 210, 360, \dots$, all of them built from the smallest primes.

2. Where the solutions lie

Since $\frac1x\lt\frac1x+\frac1y=\frac1n$, we have $x\gt n$, and likewise $y\gt n$. If $x\le y$, then $\frac1n\le\frac2x$, so $x\le2n$. The smaller number of a solution therefore lies among $n+1,\dots,2n$, and it determines the other one: $\frac1y=\frac1n-\frac1x=\frac{x-n}{nx}$, so

$$y=\frac{nx}{x-n},$$

and $x$ belongs to a solution exactly when this is a whole number. For $n=4$: $x=5$ gives $y=20$, $x=6$ gives $y=12$, $x=7$ gives $28/3$, which is not whole, and $x=8$ gives $y=8$, the three solutions of the statement.

This is already an algorithm: for $n=1,2,3,\dots$ try all $n$ values of $x$, and stop at the first $n$ with more than $1000$ solutions. It is correct but slow. Section 7 shows that the answer is at most $510510$, and testing every $n$ up to there costs up to $1+2+\dots+510510$, about $1.3\cdot10^{11}$ divisions.

3. From a sum of fractions to a product

Multiplying $\frac1x+\frac1y=\frac1n$ by $nxy$ gives $ny+nx=xy$. Moving everything to one side and adding $n^2$ to both sides gives $xy-nx-ny+n^2=n^2$, and the left side factors:

$$(x-n)(y-n)=n^2.$$

Put $a=x-n$ and $b=y-n$; by Section 2 both are positive, and $ab=n^2$. Conversely, any positive integers with $ab=n^2$ give a solution $x=n+a$, $y=n+b$: since $(n+a)(n+b)=n^2+n(a+b)+ab$${}=n(2n+a+b)$,

$$\begin{aligned}\frac1{n+a}+\frac1{n+b}&=\frac{2n+a+b}{(n+a)(n+b)}\\&=\frac1n.\end{aligned}$$

So the solutions $(x,y)$, taken in both orders, correspond one to one to the ordered pairs $(a,b)$ of positive integers with $ab=n^2$. For $n=4$ we have $16=1\cdot16=2\cdot8=4\cdot4$, and adding $4$ to both factors gives $(5,20)$, $(6,12)$ and $(8,8)$. In general the factor $a=1$ gives the splitting $\frac1n=\frac1{n+1}+\frac1{n(n+1)}$, and the factor $a=n$ gives $x=y=2n$.

4. Counting the factor pairs

An ordered pair $(a,b)$ with $ab=n^2$ is fixed by $a$, which can be any divisor of $n^2$, so there are $d(n^2)$ ordered pairs, where $d(m)$ is the number of positive divisors of $m$. The pairs $(a,b)$ and $(b,a)$ give the same distinct solution, and they are different pairs unless $a=b$, that is $a^2=n^2$, or $a=n$. The other $d(n^2)-1$ divisors fall into pairs $\{a,\,n^2/a\}$, one member below $n$ and one above. Write $s(n)$ for the number of distinct solutions. Together with the solution $x=y=2n$, it is $(d(n^2)-1)/2+1$, that is

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

This also shows why a square always has an odd number of divisors: its square root is the only divisor without a partner. ^[4] For $n=6$, the nine divisors of $36$ pair up as $1\cdot36$, $2\cdot18$, $3\cdot12$ and $4\cdot9$, with $6\cdot6$ left over, which gives the five solutions of Section 1. For $n=12$, the $15$ divisors of $144$ give $s(12)=8$.

5. Counting divisors with the prime factorisation

By the fundamental theorem of arithmetic, every integer $n\ge2$ is a product of prime powers in exactly one way apart from the order of the factors, $n=p_1^{e_1}\cdots p_k^{e_k}$ with different primes $p_i$ and exponents $e_i\ge1$; $n=1$ is the empty product, with $k=0$. ^[5] By the same uniqueness, a divisor of $m=p_1^{E_1}\cdots p_k^{E_k}$ contains no other prime and no $p_i$ more often than $m$ does, so it has the form $p_1^{f_1}\cdots p_k^{f_k}$ with $0\le f_i\le E_i$; every choice of these exponents gives a divisor, and different choices give different divisors. So $d(m)=(E_1+1)\cdots(E_k+1)$. ^[4] The square $n^2$ has the exponents $2e_i$, and therefore

$$d(n^2)=(2e_1+1)\cdots(2e_k+1).$$

For example, $144=2^4\cdot3^2$ has $5\cdot3=15$ divisors $2^f3^g$, one for each $f\in\{0,\dots,4\}$ and $g\in\{0,1,2\}$. A prime $n=p$ has $d(p^2)=3$ and $s(p)=2$, the solutions $\bigl(p+1,\,p(p+1)\bigr)$ and $(2p,2p)$, which explains the many twos of Section 1. Together with Section 4 this gives a closed form for the number of solutions:

$$s(n)=\frac{(2e_1+1)\cdots(2e_k+1)+1}{2}.$$

For $n=1$ it gives $s(1)=1$, since an empty product is $1$.

6. What "more than one thousand" means

We need $s(n)\gt1000$, that is $d(n^2)+1\gt2000$, or $d(n^2)\gt1999$. Since $d(n^2)$ is odd, the smallest allowed value is $2001$:

$$s(n)\gt1000\iff d(n^2)\ge2001.$$

A common slip is to count the ordered pairs $(x,y)$, that is, to ask for $d(n^2)\gt1000$. Such a program stops too early: the least $n$ with $d(n^2)\gt1000$ is $55440=2^4\cdot3^2\cdot5\cdot7\cdot11$, with $9\cdot5\cdot3\cdot3\cdot3=1215$ ordered solutions but only $608$ distinct ones.

7. A first candidate and a faster scan

If all exponents are $1$, every prime contributes a factor $3$ to $d(n^2)$. Since $3^6=729\lt2001\le2187=3^7$, seven primes are then needed, and the smallest choice is

$$2\cdot3\cdot5\cdot7\cdot11\cdot13\cdot17=510510,$$

the seventh primorial, that is, the product of the first seven primes. ^[6] It has $d(n^2)=2187$ and $s(n)=1094$, so the least $n$ with more than $1000$ solutions is at most $510510$. The real question is whether a smaller one exists.

The closed form gives a second algorithm: scan $n=1,2,3,\dots$, factor each $n$, and stop at the first $n$ with $s(n)\gt1000$. Trial division tries the divisors $2,3,4,\dots$ up to $\sqrt n$, at most $713$ of them for $n\le510510$, because $714^2\le510510\lt715^2$, and it divides out at most $18$ prime factors; so the scan needs fewer than $4\cdot10^8$ divisions. Stopping at the square root is safe: if $m=ab$ with $1\lt a\le b$, then $a\le\sqrt m$, so a number without a divisor from $2$ up to its square root is $1$ or a prime, and so is the part of $n$ that remains at the end. The sieve of Eratosthenes, which finds the primes up to a limit by crossing out the multiples of each prime, ^[7] can also record for every composite number the first prime that crosses it out, its smallest prime factor, and then each factorisation takes only a few look-ups. Either scan is exact, but it examines hundreds of thousands of numbers, while only a few hundred can be the answer.

8. Only the exponents matter

The closed form does not care which primes carry the exponents: $24=2^3\cdot3$, $40=2^3\cdot5$, $54=2\cdot3^3$ and $250=2\cdot5^3$ all have the exponents $3$ and $1$, so all of them have $d(n^2)=7\cdot3=21$ and $11$ solutions. Among numbers with the same exponents, the smallest one puts the largest exponent on $2$, the next largest on $3$, and so on, without skipping a prime.

The proof is an exchange argument. Let $p\lt q$ be primes whose exponents in $n$ are $e\lt f$, where $e=0$ means that $p$ does not divide $n$. Swapping the two exponents does not change the list of exponents, so it does not change $s(n)$, but it makes $n$ smaller:

$$p^eq^f=p^fq^e\Bigl(\frac qp\Bigr)^{f-e}\gt p^fq^e.$$

If, along the primes $2,3,5,7,\dots$, some exponent were larger than an earlier one, counting $0$ for a missing prime, such a swap would give a smaller number with the same $s(n)$. So the least $n$ with $s(n)\gt1000$ has the form

$$n=2^{a_1}3^{a_2}5^{a_3}\cdots p_k^{a_k},$$

with the first $k$ primes and $a_1\ge a_2\ge\dots\ge a_k\ge1$. Ramanujan's highly composite numbers, which have more divisors than every smaller number, have this shape for the same reason. ^[8] The converse is false, and the two kinds of record differ: $30$ has more solutions than every smaller number (Section 1), but $d(30)=8=d(24)$, so $30$ is not highly composite, and the search has to consider every number of this shape. Numbers of this shape are products of primorials: $a_1-a_2$ factors $2$, $a_2-a_3$ factors $2\cdot3$, and so on; for example, $360=2^3\cdot3^2\cdot5=2\cdot6\cdot30$. The On-Line Encyclopedia lists them as the least numbers of each prime signature, that is, of each list of exponents: $1, 2, 4, 6, 8, 12, 16,$ $24, 30, 32, 36, 48, 60, \dots$ ^[9]

Up to $510510$ there are only $246$ of them, counting $1$, and only $10$ of these reach $d(n^2)\ge2001$; the answer is the smallest of those ten. Moreover, the product of the first eight primes is $9699690\gt510510$, so the answer uses at most the first seven primes, and $2^{19}=524288\gt510510$, so its largest exponent is at most $18$.

9. Searching the exponent patterns

The numbers of Section 8 form a tree. The root is $n=1$; going down one level chooses the exponent of the next prime, at least $1$ and, after the prime $2$, at most the exponent of the previous prime, and multiplies that prime power in. Every node knows its number $n$ and the product $D=(2a_1+1)\cdots(2a_j+1)$ for the exponents chosen so far, which is $d(n^2)$. The tree is explored by depth-first search: exponent $1$ for the current prime and everything below it, then exponent $2$, and so on. ^[10] Two rules keep the tree small.

Stop at enough divisors. If $D\ge2001$, the node is a candidate: it is recorded if it beats the best candidate found so far, and the search does not go deeper, since further primes only make $n$ larger.
Stop at a worse number. If multiplying in the current prime makes $n$ at least as large as the best candidate so far, nothing below can do better, since numbers only grow on the way down; higher exponents of the same prime are larger still, so the loop over the exponents ends too.

Discarding a part of a search tree that cannot beat the best solution found so far is the idea of branch and bound. ^[11] The two rules never lose the answer $n^*$. By Section 8 it is a node of the tree. No node before it on its path has $D\ge2001$, since that node would be a smaller number with enough divisors, so the first rule does not end the path early. At each node of the path, with next prime $p$, the loop forms the node times $p,p^2,\dots$ up to the power that gives the next node of the path; all these numbers are at most $n^*$. The best candidate stays above $n^*$ until $n^*$ itself is recorded, since every number appears only once in the tree. So the second rule never stops these loops before they reach the next node of the path.

10. Worked example: more than ten solutions

For more than $10$ solutions we need $D\ge21$. The search runs as follows:

$2$ with $D=3$, then $6$ with $D=9$, then $30$ with $D=27\ge21$: the first candidate, so the best is $30$.
Exponent $2$ for the prime $2$: $4$ with $D=5$, then $12$ with $D=15$. Multiplying by $5$ gives $60\ge30$: stop. Exponent $2$ for the prime $3$ gives $36\ge30$: stop.
Exponent $3$: $8$ with $D=7$, then $24$ with $D=21$, a candidate below $30$, so the best is now $24$. Exponent $2$ for the prime $3$ gives $72\ge24$: stop.
Exponent $4$: $16$ with $D=9$, then $48\ge24$: stop. Exponent $5$ gives $32\ge24$, which ends the search.

The search visited $9$ nodes, counting the root, found $2$ candidates and stopped $5$ times. The result $24=2^3\cdot3$ has $d(576)=21$ and the $11$ solutions $(25,600)$, $(26,312)$, $(27,216)$, $(28,168)$, $(30,120)$, $(32,96)$, $(33,88)$, $(36,72)$, $(40,60)$, $(42,56)$ and $(48,48)$, while every $n\lt24$ has at most $8$. The example also shows a pitfall: the first candidate, the primorial $30$, is not the least one, so the primorial of Section 7 need not be the answer either. In the same way the search gives $6$ for more than $4$ solutions, the checkpoint of the C++ program, and $1260=2^2\cdot3^2\cdot5\cdot7$, with $113$ solutions, for more than $100$, the value quoted in Problem 110. ^[2]

11. The search for one thousand

For $D\ge2001$ the first path of the search takes exponent $1$ for every prime: it passes $2,6,30,210,2310,30030$ with $D=3,9,27,81,243,729$ and reaches $510510$ with $D=2187$, the bound of Section 7. After that the second rule cuts almost everything: the whole search visits $187$ nodes, finds only $2$ candidates, stops $107$ times and never uses more than seven primes. Without the second rule, the tree that the programs search, in which the exponent of $2$ is at most $63$, would have $8197$ nodes and $6790$ candidates, with numbers of up to $50$ digits; without that cap the loop over the powers of $2$ would never end. The same C++ program, given the threshold of Problem 110, four million, still answers within a few milliseconds.

How the Code Works

The C++ program ^[12] starts in main. parse_arguments reads --target=K, which replaces the default threshold $1000$ (parse_int_after_prefix accepts only digits), and --skip-checkpoints; anything else, or a threshold below $1$, gives exit code $1$. run_checkpoints requires count_solutions(4) to be $3$, the example of the statement, and solve(4) to be $6$, the value of Section 10, and gives exit code $2$ otherwise. count_solutions, which only this checkpoint uses, evaluates the closed form of Section 5 by trial division (Section 7): it multiplies $2e+1$ for each prime power $p^e$ it divides out and $3$ for a prime left over at the end. Composite trial divisors never divide, because their prime factors are already gone.

solve(target) holds the first $15$ primes in a std::array and a SearchState with the best number found so far, best_n, and need_divisors, which is 2 * target + 1, that is $2001$, as in Section 6. It calls dfs_min_n(primes, 0, 63, 1, 1, state): first prime, exponents up to $63$, number $1$ and $D=1$. dfs_min_n is the search of Section 9. If divisor_count >= need_divisors, the first rule, it updates best_n when the number is smaller and returns; it also returns when all $15$ primes are used. Otherwise it multiplies the number by the current prime for the exponents $1,2,\dots$ up to max_exp, leaves the loop with break as soon as the number reaches best_n, the second rule, and else calls itself for the next prime, with the current exponent as the new max_exp and the divisor count multiplied by 2 * exp + 1. Passing the exponent down as the new limit is what keeps the exponents non-increasing.

The partial numbers are of type u128, unsigned 128-bit integers, the result is a u64, and best_n starts at $2^{64}-1$. They need 128 bits because each product is formed before it is compared with best_n: a number below $2^{64}$ times a prime up to $47$ stays below $2^{70}$, while a 64-bit product could wrap around to a small number below best_n, and from a threshold of about thirteen million on the search would return wrong numbers. The type unsigned __int128 is a GCC and Clang extension, not standard C++. Since $2\cdot3\cdots47\approx6.1\cdot10^{17}$ is below $2^{64}\approx1.8\cdot10^{19}$ and $2\cdot3\cdots53\approx3.3\cdot10^{19}$ is above it, a number of the shape of Section 8 below $2^{64}$ uses at most $15$ primes and an exponent of at most $63$: the two constants match the 64-bit result exactly.

The largest $d(n^2)$ for $n\lt2^{64}-1$ is $72335025$, so the program answers every threshold up to $36167512$, and the divisor counts it forms never exceed that value and fit into an int. For thresholds from $36167513$ up to $2^{30}-1$ it prints its starting value $2^{64}-1$ unchanged instead of reporting that the answer does not fit. With the usual 32-bit int, a threshold from $2^{30}$ to $2^{31}-1$ makes 2 * target + 1 overflow, and a digit string whose value exceeds $2^{31}-1$ already overflows in the digit loop of parse_int_after_prefix, which has no range check. Signed overflow is undefined behaviour in C++, so for these inputs neither the output nor the exit code is guaranteed.

The Python and Java versions run the same search with the divisor bound $2001$, that is, the threshold $1000$, built in, without options or checkpoints. In Python the last line calls solve, which sets best to float('inf'), runs the nested function dfs(0, 63, 1, 1), which updates best as a nonlocal variable, and prints the result; Python's integers have no size limit. The comment above the search, d(n²) > 2*1000-1 = 2001, contains a slip, since $2\cdot1000-1=1999$: the condition of Section 6 is $d(n^2)\gt1999$, that is $d(n^2)\ge2001$, and this is what the code tests with target = 2 * 1000 + 1 and divisor_count >= target. Java's main sets the long variable best to Long.MAX_VALUE, calls dfs(0, 63, 1, 1, 2001) and prints best; the numbers are long and the divisor counts int. Because the first candidate, $510510$, appears on the first path, every later product is a number below $510510$ times a prime of at most $47$; so no product reaches $47\cdot510510$, far below $2^{63}$, and the largest one Java actually forms is $2162160$.

Complexity and Verification

Every multiplication by a prime either leads to a new node or is stopped by the second rule, which ends that node's loop; so with $V$ nodes there are at most $V-1$ multiplications of the first kind and $V$ of the second, and the time is proportional to the number of nodes. The memory is the current path of the recursion, at most $16$ calls: the root and one call for each of the $15$ primes. For the threshold $1000$ the search visits $187$ nodes and makes $186+107=293$ multiplications by a prime, as it stops $107$ times. On a desktop computer the C++ program, checkpoints included, runs in a few milliseconds, and the Python and Java versions in about $14$ and $26$ milliseconds; an empty program takes almost as long in each of the three languages, so nearly all of this time is start-up. We know no simple formula for the number of nodes, but every node is a different number of the shape of Section 8, and such numbers are rare: $246$ up to $510510$, which bounds the search for $1000$, and $47616$ below $2^{64}$, which bounds the C++ search for every threshold.

Unless --skip-checkpoints is given, the C++ program runs its two checkpoints before it searches, and all three programs print the same number. For this article the mathematics was also checked by separate computations that are not part of the published solutions:

a brute force over $x$ agrees with $s(n)=(d(n^2)+1)/2$ for every $n\le5000$, and a search over all $x\le y\le n^2+n$, without the restriction $n\lt x\le2n$, finds the same solutions for every $n\le20$;
the divisor formula agrees with a direct count for every $m\le3000$, and the counts for $n\le30$ agree with the OEIS entry; ^[3]
a scan of every $n\le510510$, factored with a smallest-prime-factor sieve, gives for every threshold from $1$ to $1000$ the same least $n$ as the search; for $1000$ it is the number that the three programs print;
the C++ option --target agrees with the same search on the first $40$ primes for every threshold up to $200$, for several larger ones up to $10^5$, and for four million.

Footnotes and References

The references below supply the background facts used above; the factor pairs, the closed form, the shape of the answer and the correctness of the search are derived in this article.

Project Euler 108 — Diophantine Reciprocals I. The official statement: x, y and n are positive integers with 1/x + 1/y = 1/n; for n = 4 there are exactly three distinct solutions, and the task asks for the least n for which the number of distinct solutions exceeds one thousand. It calls the problem an easier version of Problem 110.
Project Euler 110 — Diophantine Reciprocals II. The same equation; n = 1260, with 113 distinct solutions, is stated to be the least n with more than one hundred, and the task asks for the least n with more than four million.
OEIS A018892 — Number of ways to write 1/n as a sum of exactly 2 unit fractions. The entry begins 1, 2, 2, 3, 2, 5, 2, 4, 3, 5, 2, 8 and gives the formula a(n) = (tau(n²) + 1)/2, where tau counts divisors.
Wikipedia — Divisor function. The number of divisors, d(n) or σ₀(n), is the product of the numbers aᵢ + 1 over the prime powers pᵢ^aᵢ of n; for a square n the divisor √n has no different partner, so σ₀(n) is odd.
Wikipedia — Fundamental theorem of arithmetic. Every integer greater than 1 is a prime or a product of primes, unique up to the order of the factors; equivalently, a product of prime powers in exactly one way.
Wikipedia — Primorial. The primorial pₙ# is the product of the first n primes; the values begin 1, 2, 6, 30, 210, 2310, 30030, 510510, 9699690.
Wikipedia — Sieve of Eratosthenes. An ancient algorithm that finds all primes up to a given limit by marking the multiples of each prime, in O(n log log n) operations.
Wikipedia — Highly composite number. Such a number has more divisors than every smaller positive integer, and Ramanujan wrote a paper on them in 1915; its prime factors are the first k primes, its exponents are non-increasing, and so it is a product of primorials.
OEIS A025487 — Least integer of each prime signature; also products of primorial numbers. The terms are the numbers 2^k₁ · 3^k₂ · … · pₙ^kₙ with k₁ ≥ k₂ ≥ … ≥ kₙ, in increasing order: 1, 2, 4, 6, 8, 12, 16, 24, 30, …
Wikipedia — Depth-first search. An algorithm for traversing or searching a tree or graph: it starts at the root and follows each branch as deep as it goes before it backtracks.
Wikipedia — Branch and bound. A method for optimisation problems that explores a tree of subproblems and discards a branch if it cannot produce a better solution than the best one found so far.
C++ — Euler108.cpp. The linked immutable revision contains count_solutions, dfs_min_n, solve, run_checkpoints and the options --target= and --skip-checkpoints. It is the implementation record for this article.

Problem 108 source code

C++

#include <array>
#include <cstdint>
#include <iostream>
#include <limits>
#include <string>

namespace {

using u64 = std::uint64_t;
using u128 = unsigned __int128;

struct Options {
    int target = 1000;
    bool run_checkpoints = true;
};

struct SearchState {
    u128 best_n = std::numeric_limits<u64>::max();
    int need_divisors = 0;
};

bool parse_int_after_prefix(const std::string& arg, const std::string& prefix, int& value) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }

    int parsed = 0;
    for (char c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }
        parsed = parsed * 10 + static_cast<int>(c - '0');
    }
    value = parsed;
    return true;
}

bool parse_arguments(int argc, char** argv, Options& options) {
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);

        if (arg == "--skip-checkpoints") {
            options.run_checkpoints = false;
            continue;
        }
        if (parse_int_after_prefix(arg, "--target=", options.target)) {
            continue;
        }

        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }

    return options.target >= 1;
}

u64 count_solutions(const u64 n) {
    u64 x = n;
    u64 divisors_n2 = 1;

    for (u64 p = 2; p * p <= x; ++p) {
        if (x % p != 0) {
            continue;
        }
        int exp = 0;
        while (x % p == 0) {
            x /= p;
            ++exp;
        }
        divisors_n2 *= static_cast<u64>(2 * exp + 1);
    }
    if (x > 1) {
        divisors_n2 *= 3;
    }

    return (divisors_n2 + 1) / 2;
}

void dfs_min_n(const std::array<int, 15>& primes,
               int idx,
               int max_exp,
               u128 current_n,
               int divisor_count,
               SearchState& state) {
    if (divisor_count >= state.need_divisors) {
        if (current_n < state.best_n) {
            state.best_n = current_n;
        }
        return;
    }

    if (idx >= static_cast<int>(primes.size())) {
        return;
    }

    u128 n = current_n;
    for (int exp = 1; exp <= max_exp; ++exp) {
        n *= static_cast<u128>(primes[static_cast<std::size_t>(idx)]);
        if (n >= state.best_n) {
            break;
        }

        const int new_divisor_count = divisor_count * (2 * exp + 1);
        dfs_min_n(primes, idx + 1, exp, n, new_divisor_count, state);
    }
}

u64 solve(const int target) {
    static constexpr std::array<int, 15> primes = {
        2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47,
    };

    SearchState state;
    state.need_divisors = 2 * target + 1;
    dfs_min_n(primes, 0, 63, 1, 1, state);

    return static_cast<u64>(state.best_n);
}

bool run_checkpoints() {
    if (count_solutions(4) != 3) {
        std::cerr << "Checkpoint failed for n=4 solution count" << '\n';
        return false;
    }
    if (solve(4) != 6) {
        std::cerr << "Checkpoint failed for target=4" << '\n';
        return false;
    }
    return true;
}

}  // namespace

int main(int argc, char** argv) {
    Options options;
    if (!parse_arguments(argc, argv, options)) {
        return 1;
    }
    if (options.run_checkpoints && !run_checkpoints()) {
        return 2;
    }

    std::cout << solve(options.target) << '\n';
    return 0;
}

Python

# Problem 108: Diophantine reciprocals I
# Find the least n such that 1/x + 1/y = 1/n has more than 1000 distinct solutions.

def solve():
    # Number of solutions = (d(n²) + 1) / 2 where d = number of divisors
    # We need smallest n with d(n²) > 2*1000-1 = 2001
    # DFS over exponents in prime factorization
    primes = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47]
    target = 2 * 1000 + 1
    best = float('inf')
    
    def dfs(idx, max_exp, n, divisor_count):
        nonlocal best
        if divisor_count >= target:
            best = min(best, n)
            return
        if idx >= len(primes):
            return
        val = n
        for e in range(1, max_exp + 1):
            val *= primes[idx]
            if val >= best:
                break
            dfs(idx + 1, e, val, divisor_count * (2 * e + 1))
    
    dfs(0, 63, 1, 1)
    print(best)

solve()

Java

public class Euler108 {
    static int[] primes = { 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47 };
    static long best;

    static void dfs(int idx, int maxExp, long n, int divCount, int target) {
        if (divCount >= target) {
            best = Math.min(best, n);
            return;
        }
        if (idx >= primes.length)
            return;
        long val = n;
        for (int e = 1; e <= maxExp; e++) {
            val *= primes[idx];
            if (val >= best)
                break;
            dfs(idx + 1, e, val, divCount * (2 * e + 1), target);
        }
    }

    public static void main(String[] args) {
        best = Long.MAX_VALUE;
        dfs(0, 63, 1, 1, 2001);
        System.out.println(best);
    }
}