Problem 110: Diophantine Reciprocals II

Project Euler Problem 110 Solution

EulerSolve provides an optimized solution for Project Euler Problem 110, Diophantine Reciprocals II, 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 $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....

Detailed mathematical approach

Problem Summary

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

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

For $n=1260$ there are $113$ distinct solutions, and $1260$ is the least $n$ with more than one hundred of them. We must find the least $n$ for which the number of distinct solutions exceeds four million. ^[1] Two solutions that differ only in the order of $x$ and $y$ count as the same solution, as the example of $1260$ confirms (Section 3), so we may always assume $x\le y$. The statement calls this a much more difficult version of Problem 108, which asks the same question for more than one thousand solutions. The equation itself is known as the optic equation: two resistors of $x$ and $y$ ohms connected in parallel act like a single resistor of $n$ ohms, and for a thin lens with focal length $n$, an object at distance $x$ from the lens has its image at distance $y$. ^[2]

The whole solution in one sentence: the solutions for a given $n$ correspond exactly to the ways of splitting $n^2$ into a product of two factors, so their number depends only on the exponents in the prime factorisation of $n$, and the least $n$ with enough solutions is found by a short search that gives the larger exponents to the smaller primes and abandons every branch that can no longer beat the best number found so far. The sections below build this up from tiny cases, so that every step can be checked by hand.

Mathematical Approach

We first count the solutions for one fixed $n$, then turn the question into a question about prime exponents, and finally organise a search for the best exponents that is complete and still short.

1. Small values of $n$

For a single $n$ the solutions can be found by hand. Since $1/y$ is positive, $1/x$ is smaller than $1/n$, so $x\gt n$. Since $x\le y$, the fraction $1/x$ is the larger of the two, so it is at least half of $1/n$, and $x\le2n$. So we only have to try $x=n+1,\dots,2n$ and check whether $y=nx/(x-n)$ is a whole number. It comes from $1/y=1/n-1/x=(x-n)/(nx)$, and it is at least $x$ because $x-n\le n$, so each $x$ with a whole $y$ gives exactly one distinct solution. For $n=6$:

$x=7$ gives $y=42/1=42$, $x=8$ gives $y=48/2=24$ and $x=9$ gives $y=54/3=18$;
$x=10$ gives $y=60/4=15$, $x=11$ gives $66/5$, which is not a whole number, and $x=12$ gives $y=72/6=12$.

So $n=6$ has five solutions: $(7,42)$, $(8,24)$, $(9,18)$, $(10,15)$ and $(12,12)$. For instance, $\tfrac1{10}+\tfrac1{15}$${}=\tfrac3{30}+\tfrac2{30}=\tfrac16$. In the same way $n=1$ has the single solution $(2,2)$, $n=2$ has $(3,6)$ and $(4,4)$, and $n=4$ has $(5,20)$, $(6,12)$ and $(8,8)$. The numbers of solutions for $n=1,2,\dots,12$ are

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

The size of $n$ is clearly not what matters. The prime $11$ has only $2$ solutions, like every other prime, while $6$ has $5$, and so do $10$, $14$ and $15$, the other products of two different primes up to $15$. The squares of primes $4$ and $9$ have $3$ solutions each, and $12=2^2\cdot3$ has $8$. The count seems to depend only on the shape of the prime factorisation of $n$. The next three sections show why, and give the exact rule.

This direct count also solves the example of the statement: checking $n=1,2,\dots,1260$ in turn takes $1+2+\dots+1260=794\,430$ trials of $x$. It cannot go much further. Already the least $n$ with more than $10^5$ solutions is $53\,542\,288\,800$, as Section 11 shows, so the same scan would need about $1.4\cdot10^{21}$ trials there, and more than four million solutions need a larger $n$ still.

2. Turning the equation into a product

Multiplying both sides of the equation by $nxy$ clears the fractions:

$$ny+nx=xy.$$

Moving everything to one side gives $xy-nx-ny=0$, and adding $n^2$ to both sides completes a product:

$$(x-n)(y-n)=xy-nx-ny+n^2=n^2.$$

We saw that $x\gt n$, and in the same way $y\gt n$, so $u=x-n$ and $v=y-n$ are positive integers with $uv=n^2$. Conversely, if $u$ and $v$ are positive integers with $uv=n^2$, then $x=n+u$ and $y=n+v$ form a solution, because

$$\frac1{n+u}+\frac1{n+v}=\frac{2n+u+v}{n^2+n(u+v)+uv},$$

and when $uv=n^2$ the denominator equals $2n^2+n(u+v)$, which is $n\,(2n+u+v)$, so the fraction is $1/n$. The two directions undo each other, so the solutions $(x,y)$ correspond exactly to the ways of writing $n^2$ as a product $u\cdot v$ of two positive factors. For $n=6$ the factorisations $36=1\cdot36=2\cdot18=3\cdot12$${}=4\cdot9=6\cdot6$ give back the five solutions of Section 1: add $6$ to both factors. The solutions can also be described all at once as $x=cr(r+s)$, $y=cs(r+s)$ and $n=crs$, where $c$, $r$ and $s$ are positive integers and $r$ and $s$ are coprime, that is, have no common divisor larger than $1$; ^[2] for example, $(8,24)$ comes from $c=2$, $r=1$ and $s=3$. For counting, however, the factorisations of $n^2$ are much more convenient.

3. Counting the factor pairs

Every divisor $u$ of $n^2$ determines its partner $v=n^2/u$, so the number of solutions in which the order matters, with $(8,24)$ and $(24,8)$ counted separately, is $d(n^2)$, the number of divisors of $n^2$. For distinct solutions we keep only $x\le y$, that is $u\le v$. Since $uv=n^2$, this means $u^2\le n^2$, so $u\le n$. The divisors smaller than $n$ and the divisors larger than $n$ are matched in pairs $u\leftrightarrow n^2/u$, and $u=n$ is the only divisor that is its own partner. Hence $d(n^2)$ is always odd, and the number of distinct solutions is

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

For $n=6$ the divisors of $36$ are $1,2,3,4,6,9,12,18,36$: four of them are below $6$, four are above, and one is $6$ itself, so $N(6)=(9+1)/2=5$. The example of the statement settles what "distinct" means: $1260^2$ has $225$ divisors (Section 4 shows how to count them), so there are $225$ solutions in which the order matters, but $(225+1)/2=113$ distinct ones, as stated. The numbers $N(1),N(2),N(3),\dots$ form entry A018892 of the On-Line Encyclopedia of Integer Sequences (OEIS), the number of ways to write $1/n$ as a sum of two unit fractions, and the entry gives the same formula. ^[3]

4. Counting divisors from the prime factorisation

Every whole number larger than $1$ is a product of primes in exactly one way, apart from the order of the factors; this is the fundamental theorem of arithmetic. ^[4] Take $m=p_1^{e_1}p_2^{e_2}\cdots p_k^{e_k}$ with different primes $p_i$. A divisor of $m$ can contain no other prime, and it cannot contain $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$; conversely, every such product divides $m$. By the uniqueness of the factorisation, different choices of the exponents $f_i$ give different numbers. There are $e_i+1$ choices for each $f_i$, so

$$d(m)=(e_1+1)(e_2+1)\cdots(e_k+1).$$

This is the standard formula for the divisor function. ^[5] For $36=2^2\cdot3^2$ the divisors $2^f3^g$ with $f,g\in\{0,1,2\}$ fill a square of three rows and three columns: $1,2,4$ in the first row, $3,6,12$ in the second and $9,18,36$ in the third, $3\cdot3=9$ divisors in all.

Now write $n=p_1^{a_1}\cdots p_k^{a_k}$. Squaring doubles every exponent, $n^2=p_1^{2a_1}\cdots p_k^{2a_k}$, so

$$D(n)=d(n^2)=(2a_1+1)(2a_2+1)\cdots(2a_k+1),$$

a product of odd numbers, which shows once more that it is odd, and the number of distinct solutions is $N(n)=(D(n)+1)/2$. The list of Section 1 follows at once. The number $1$ has no prime factors at all, and a product without any factors, the empty product, is $1$ by convention; so $D(1)=1$ and $N(1)=1$, the single solution $(2,2)$. This is also why the search of Section 9 starts from the number $1$ with $D=1$. A prime has $D=3$ and $N=2$; the square of a prime has $D=5$ and $N=3$; a product of two different primes has $D=3\cdot3=9$ and $N=5$; $8=2^3$ has $D=7$ and $N=4$; and $12=2^2\cdot3$ has $D=5\cdot3=15$ and $N=8$. For $1260=2^2\cdot3^2\cdot5\cdot7$ the rule gives $D=5\cdot5\cdot3\cdot3=225$ and $N=113$, the value of the statement.

5. The question in terms of exponents

The number of distinct solutions exceeds four million exactly when $(D+1)/2\gt4\,000\,000$, that is, when $D\gt7\,999\,999$. Since $D$ is odd, it cannot be $8\,000\,000$, so the condition is

$$D(n)=\prod_i\,(2a_i+1)\ge8\,000\,001.$$

In general, more than $T$ solutions means $D(n)\ge2T+1$; for the example of the statement, more than one hundred solutions means $D(n)\ge201$. We write $K$ for this bound. The variables $x$ and $y$ have disappeared: we look for the least number $n=\prod p_i^{a_i}$ whose exponents make the product $\prod(2a_i+1)$ at least $K$.

It helps to think of $n$ as a price and of $D$ as what the price buys. Raising the exponent of a prime $p$ from $a$ to $a+1$ multiplies the price by $p$ and multiplies $D$ by $(2a+3)/(2a+1)$: the first factor $p$ multiplies $D$ by $3$, the second one only by $5/3$, the third one by $7/5$, and so on. Small primes are cheap, fresh primes give the biggest gain, and the whole difficulty lies in balancing the two.

6. Small primes first, large exponents first

Two exchange arguments show what the least $n$ with $D(n)\ge K$ must look like, for any bound $K$.

Its primes are the first primes $2,3,5,\dots$, without gaps. If a prime $q$ divided $n$ while a smaller prime $p$ did not, we could replace every factor $q$ by $p$. The exponents stay the same, so $D$ stays the same, but the number becomes smaller, which is impossible for the least $n$.
Its exponents never increase: $a_1\ge a_2\ge\dots\ge a_k\ge1$, where $a_1$ is the exponent of $2$, $a_2$ that of $3$, and so on. If a smaller prime $p$ had a smaller exponent $a$ than a larger prime $q$ with exponent $b$, swapping the two exponents would keep $D$ and multiply the number by $(p/q)^{b-a}$, which is less than $1$.

For example, $18=2\cdot3^2$ and $12=2^2\cdot3$ both have $D=15$, so $18$ can never be the least number for any bound: $12$ does the same job for less. The same two arguments are used for highly composite numbers, the numbers with more divisors than every smaller number, and they lead to the same shape. ^[6]

For a prime $p$, the primorial $p\#$ is the product of all primes up to $p$: $2\#=2$, $3\#=6$, $5\#=30$, $7\#=210$, and so on. A number of our shape is a product of primorials: $2^{a_1}3^{a_2}5^{a_3}\cdots$ contains $2\#$ exactly $a_1-a_2$ times, $3\#$ exactly $a_2-a_3$ times, and so on. For example, $1260=2^2\cdot3^2\cdot5\cdot7=6\cdot210$. These numbers are also the least numbers of each prime signature, that is, of each list of exponents, and some authors call them Hardy–Ramanujan integers. They begin $1,2,4,6,8,12,16,24,30,\dots$ ^[7] They are rare: there are only $289$ of them up to $10^6$, $4357$ up to $10^{12}$ and $32\,749$ up to $10^{18}$.

7. How many primes, how large the exponents

Simple bounds make the search finite. The primorial $47\#$, the product of the first $15$ primes, has all exponents equal to $1$, so

$$D(47\#)=3^{15}=14\,348\,907\ge8\,000\,001.$$

The answer is therefore at most

$$47\#=614\,889\,782\,588\,491\,410.$$

The product $43\#$ of the first $14$ primes would not do, since $D(43\#)=3^{14}=4\,782\,969$ is too small, but larger exponents can replace some of the primes. A number with $16$ or more different prime factors is at least the product of the first $16$ primes, ^[8]

$$53\#=32\,589\,158\,477\,190\,044\,730,$$

which is larger than $47\#$. Hence the answer has at most $15$ different prime factors, and by Section 6 they are among $2,3,\dots,47$. Its exponent of $2$ is at most $59$, because $2^{60}=1\,152\,921\,504\,606\,846\,976$ is already larger than $47\#$, while $2^{59}=576\,460\,752\,303\,423\,488$ is smaller; by Section 6 every other exponent is at most that of $2$. These bounds justify the choices made in the programs: they use the $16$ primes from $2$ to $53$ and allow exponents up to $63$, a little more than necessary.

8. Why a simple rule is not enough

It is tempting to build $n$ step by step, each time multiplying by the prime that buys the most for its price. On a logarithmic scale, where prices and gains add up, this is the prime with the largest ratio $\ln\frac{2a+3}{2a+1}\big/\ln p$, where $a$ is the current exponent of $p$. For more than one hundred solutions, $D\ge201$, this greedy rule multiplies by $2,3,2,5,7,2,3$ in this order and passes through

$$2,\ 6,\ 12,\ 60,\ 420,\ 840,\ 2520$$

with $D=3,9,15,45,135,189,315$. It stops at $2520=2^3\cdot3^2\cdot5\cdot7$, twice the true answer $1260=2^2\cdot3^2\cdot5\cdot7$. Its last step jumps from $189$ far beyond $201$, while $1260$ reaches $225$, just enough. The rule fails even for more than two solutions: from $2$ it moves to $6$, because a fresh prime $3$ triples $D$ for the price $3$, and ends with $5$ solutions, while the least answer is $4=2^2$, with $3$ solutions.

The best numbers are not even built on each other. The least $n$ with more than $50$ solutions is $360=2^3\cdot3^2\cdot5$, with $53$ solutions, but the least $n$ with more than $100$ solutions, $1260$, is not a multiple of $360$: the exponent of $2$ went down from $3$ to $2$ when the target grew. So the answer for four million need not extend the answer for a smaller target, and we cannot rely on building it that way. We need a search that is complete, so that no candidate of Section 6 can be missed, and that rules out whole groups of candidates at once.

9. The search tree

Every candidate of Section 6 is reached from $1$ by choosing the exponent of $2$, then the exponent of $3$, which is not larger, then the exponent of $5$, and so on. This gives a tree. Its root is $1$, the empty list of exponents, and the children of the list $(a_1,\dots,a_j)$ are the lists $(a_1,\dots,a_j,e)$ with $1\le e\le a_j$; the children of the root are the lists $(e)$ with any $e\ge1$. Going to a child multiplies the number by $p_{j+1}^{\,e}$, where $p_{j+1}$ is the next prime, and multiplies $D$ by $2e+1$. Every candidate appears exactly once. A depth-first search walks through this tree: it follows one branch down as far as it goes, then returns to the most recent node that still has an untried child, and it tries the children in the order $e=1,2,3,\dots$.

The search keeps the best number found so far. Before anything has been found, it holds a starting value, which only has to be larger than the least number we are looking for: infinity, which is larger than every candidate, always works, and so does any number above a known number with $D\ge K$; for four million, any number above $47\#$ (Section 7). The search uses three rules:

At a node with $D\ge K$, record the number if it is smaller than the best one, and do not go deeper: every node below it is a multiple of this number and therefore larger.
Do not enter a child whose number is at least the best one: everything below it is at least as large.
When rule 2 rejects the child with exponent $e$, do not try larger exponents of the same prime either, because the number with exponent $e+1$ is $p$ times larger still.

For more than two solutions, $K=5$, the whole search visits four nodes, the root included, and uses all three rules. From the root $1$ it enters $2$, with $D=3$, and then $6$, with $D=9$, which rule 1 records. The exponent of $3$ may not exceed that of $2$, so the search returns to the root and enters $4$, with $D=5$, which is recorded too, since $4\lt6$. Rule 2 rejects the next child, $8$, which is not smaller than $4$, and rule 3 ends the search. The result is $4$, the answer that the greedy rule of Section 8 missed.

This method is called branch and bound: the bound is the best number found so far, and a branch is dropped as soon as it cannot beat it. ^[9] Because the search tries the exponent $1$ first everywhere, its very first path goes through $2$, $6$, $30$, $210$, $\dots$, the primorials, and, unless the starting value is smaller, it records the first primorial with $D\ge K$. For four million this is $47\#$, so the search finds the bound of Section 7 by itself and then improves on it.

The rules never lose the answer. Let $n^*$ be the least number with $D(n^*)\ge K$; by Sections 6 and 7 its list of exponents is a node of the tree. Every node $m$ on the path from the root to $n^*$ is a divisor of $n^*$. If such an $m$ other than $n^*$ had $D(m)\ge K$, it would be a smaller number with enough solutions, which is impossible; so rule 1 never ends the walk before $n^*$. Rules 2 and 3 reject only numbers that are at least as large as the current best value. That value is either the starting value, which is larger than $n^*$, or a recorded number with $D\ge K$, which is at least $n^*$. A node on the path is at most $n^*$, so it can be rejected only if it is $n^*$ itself and $n^*$ has already been recorded. Either way the search ends with the best number equal to $n^*$.

10. Worked example: more than one hundred solutions

Here $K=201$. We write each node as its number followed by its value of $D$ in brackets.

The first path takes the exponent $1$ everywhere: $2$ $(3)$, $6$ $(9)$, $30$ $(27)$, $210$ $(81)$ and $2310=2\cdot3\cdot5\cdot7\cdot11$ $(243)$. Since $243\ge201$, the best number is now $2310$.
Exponent $2$ for the prime $2$: $4$ $(5)$, $12$ $(15)$, $60$ $(45)$, $420$ $(135)$. The next step would give $420\cdot11=4620$, which is not smaller than $2310$, so rule 2 rejects it. Back at $4$, the exponent $2$ for the prime $3$ gives $36$ $(25)$, then $180$ $(75)$ and $1260$ $(225)$, and the best number becomes $1260$. The sibling of $180$, $900=2^2\cdot3^2\cdot5^2$ $(125)$, is still smaller than $1260$, but its first child $900\cdot7=6300$ is rejected by rule 2, and rule 3 skips the second one, $900\cdot49$.
Exponent $3$ for the prime $2$: the path $8$, $24$, $120$, $840$ reaches only $D=189$, and $840\cdot11=9240$ is rejected. Below $72=2^3\cdot3^2$, the node $360$ $(105)$ has the child $2520$, the result of the greedy rule, which is rejected, and $216=2^3\cdot3^3$ leads to $1080$ $(147)$, whose child $7560$ is rejected.
Exponent $4$ for the prime $2$: the path $16$, $48$, $240$ ends with the rejected child $1680=2^4\cdot3\cdot5\cdot7$. This number does have $D=9\cdot3\cdot3\cdot3=243$, enough solutions, but it is larger than $1260$. The other nodes entered in this branch are $144$, $720$ and $432$.
Exponents $5$ to $10$ for the prime $2$ produce only short branches, with the nodes $32$, $96$, $480$, $288$, $864$; $64$, $192$, $960$, $576$; $128$, $384$, $1152$; $256$, $768$; $512$; and $1024$. Every node entered there has $D$ below $201$, and the children that would reach $D\ge201$, such as $3360=2^5\cdot3\cdot5\cdot7$ with $D=297$, are at least $1260$ and are rejected. Finally $2^{11}=2048$ is at least $1260$, and rule 3 ends the whole search.

Together with the root, these are $44$ calls of the search function. A number is recorded twice, first $2310$ and then $1260$, and the result is $1260$ with $113$ solutions, as the statement says.

11. Least values for other targets

Running the same search for other targets shows how the answer grows. The least $n$ with more than $T$ solutions is:

for $T=1$, $n=2$, with $2$ solutions; for $T=2$, $n=4=2^2$, with $3$; for $T=3$ and $T=4$, $n=6$, with $5$;
for $T=5,6,7$, $n=12=2^2\cdot3$, with $8$; for $T=8,9,10$, $n=24=2^3\cdot3$, with $11$;
for $T=20$, $n=60=2^2\cdot3\cdot5$, with $23$; for $T=50$, $n=360$, with $53$; for $T=100$, $n=1260$, with $113$;
for $T=10^4$, $n=116\,396\,280$, the product of $2^3$, $3^2$ and the primes from $5$ to $19$, with $12\,758$ solutions;
for $T=10^5$, $n=53\,542\,288\,800$, the product of $2^5$, $3^2$, $5^2$ and the primes from $7$ to $23$, with $100\,238$ solutions.

The shapes change in irregular ways, as Section 8 predicted: along this list the exponent of $2$ runs through $1,2,1,2,3,2,3,2,3,5$. The number of nodes the search visits grows slowly: $44$ for $T=100$, $835$ for $T=10^4$, $2665$ for $T=10^5$ and $8576$ for $T=10^6$.

The exchange arguments of Section 6 use only that $D$ does not change when exponents are moved from one prime to another. Besides this, the search of Section 9 needs only the factor by which a new prime with exponent $e$ multiplies $D$, here $2e+1$, and bounds like those of Section 7, which are worked out from that factor. With the factor $e+1$ instead, $D$ becomes the ordinary number of divisors, a product of $k$ different primes has $2^k$ divisors instead of $3^k$, and the same search finds the least number with at least $K$ divisors; for every $K\le288$, the most divisors of any number up to $2\cdot10^6$, it agrees with a direct scan of these numbers.

How the Code Works

The C++ program ^[10] keeps its numbers in u128, an alias for the compiler's unsigned 128-bit integer type unsigned __int128, ^[11] and its settings in the structure Options: the target, $4\,000\,000$ by default, and whether to run the checkpoints. parse_arguments accepts --skip-checkpoints and --target=T; parse_int_after_prefix reads the number after the prefix digit by digit, and any other argument, or a target below $1$, makes the program stop with exit code $1$. solve_u128(target) sets the required divisor count to $2T+1$ for the target $T$, the bound $K$ of Section 5, and starts the search at the root with the sixteen primes from $2$ to $53$ and the exponent limit $63$ of Section 7. The structure SearchState holds that count and the best number found so far, best_n, which starts at $2^{64}-1$, written as static_cast<u128>(~0ULL), the largest 64-bit value converted to 128 bits. This start is the program's own choice: GCC has no 128-bit integer constants when long long is narrower than 128 bits, ^[11] but the largest 128-bit value can be written as ~static_cast<u128>(0). The last paragraph of this section shows what the choice costs.

dfs_min_n is the depth-first search of Section 9. Besides the prime list and the state, its arguments are the position of the next prime in the list, the largest exponent that prime may get, which is the exponent of the previous prime, the number built so far and its value of $D$, called divisor_count. If divisor_count has reached the required count, the function records the number in best_n when it is smaller and returns (rule 1). If all sixteen primes are used, it returns. Otherwise it multiplies a copy of the number by the next prime once per round of a loop over the exponents $1,2,\dots$ up to the limit, so no power is computed from scratch. As soon as the product reaches best_n, it leaves the loop with break (rules 2 and 3); otherwise it calls itself for the next prime, with the current exponent as the new limit, which keeps the exponents non-increasing, and with divisor_count * (2 * exp + 1) as the new value of $D$. Every product is below $53\cdot2^{64}$, well within 128 bits; in the search for four million every product even stays below $2^{61}$. The divisor counts are plain int values, and in that search they never exceed $3^{15}=14\,348\,907$.

to_string_u128 turns the result into decimal digits by repeated division by $10$, because the standard output streams have no output operator for this type: writing such a value with std::cout does not compile with GCC or Clang. run_checkpoints, which runs before the main computation unless --skip-checkpoints is given, checks that the target $4$ gives $6$, in agreement with Section 1, and that the target $1000$ gives the known least $n$ with more than one thousand solutions, the question of Problem 108; if a check fails, the program stops with exit code $2$. main reads the options, runs the checkpoints and prints the least $n$ for the chosen target.

The Python version is a single function solve that contains the same search as a nested function dfs, which updates the variable best of the enclosing function through nonlocal. The bound $2\cdot4\,000\,000+1$ is written into the code, solve prints the result at the end, and there are no options or checkpoints. Python integers have no fixed size, and best starts as float('inf'), positive infinity. The test val >= best compares a whole number with a floating-point number, and Python makes such comparisons exactly, so infinity acts as a bound above every number. ^[12] The Java version keeps the primes and the best number in static fields and uses BigInteger for the numbers, with multiply and compareTo; ^[12] its best number starts at Long.MAX_VALUE, which is $2^{63}-1$, and the bound $8\,000\,001$ is written into main. In the search for four million, all three versions visit the same nodes in the same order; the C++ program runs its two checkpoint searches before that search.

The starting values $2^{64}-1$ in C++ and $2^{63}-1$ in Java are not true infinities, but for this problem that does not matter: the first number recorded, $47\#$, is below both. The option --target of the C++ program inherits the limit, though. A number below $2^{64}$ has at most $15$ different prime factors, because $53\#\gt2^{64}$, and its exponent of $2$ is below $64$. By Section 6 such a number has the same $D$ as some product of primorials that is not larger, and among the $47\,616$ products of primorials below $2^{64}$ the largest value of $D$ is $72\,335\,025$, at $16\,484\,279\,277\,904\,237\,800$${}=2^3\cdot3^3\cdot5^2\cdot7^2$${}\cdot11\cdot13\cdots43$. So below the starting value the sixteen primes and the exponent limit $63$ always suffice and the divisor counts fit into an int, and $36\,167\,512=(72\,335\,025-1)/2$ is exactly the last target whose least $n$ lies below it; up to this target the program is right. Beyond it the program fails without any warning. Up to the target $2^{30}-1$ the least $n$ is larger than the starting value, nothing is ever recorded, and the program prints $2^{64}-1$. Larger targets break the arithmetic of the usual 32-bit int: from $2^{30}$ on, the bound $2T+1$, computed as 2 * target + 1, overflows, and from $2^{31}$ on, so does the target itself while it is read. Signed integer overflow is undefined behaviour in C++, so for these targets the language guarantees no result at all. ^[13]

Complexity and Verification

The work is proportional to the number of nodes visited. A node may try several exponents, but every round of its loop either enters a child, which is a node of its own, or is the single rejected round that ends the loop, so there are fewer rounds than twice the number of nodes, and each round costs a few multiplications and comparisons of numbers with at most $128$ bits. For the target of the problem the search function is called $17\,411$ times, a loop over exponents is cut short by rule 2 or 3 at $7384$ of these calls, so the loops run $17\,410+7384=24\,794$ rounds, one for each call except the root and one for each rejection, and a number is recorded $7$ times: first $47\#$, then six smaller ones. The number of calls is bounded in advance: every call handles a different product of primorials, and none of them exceeds $47\#$, because the first path stops at $47\#$ and afterwards only numbers below the best one are entered. There are $30\,726$ products of primorials up to $47\#$, the bound of Section 7, so the three rules skip $13\,315$ of them. Listing all of them and keeping the least one with $D\ge8\,000\,001$ would also solve the problem, but the search does not need to know the bound beforehand. The memory is the recursion stack: at most $17$ nested calls, one for the root and one for each of the sixteen primes; in the search for four million at most $16$ are nested, on the path to $47\#$, which uses fifteen primes. On the computer used for these checks, each of the three programs finished in less than a tenth of a second, start-up included.

The method would also cope with much larger targets. Hardy and Ramanujan proved that the number of products of primorials up to $x$ is

$$\exp\Bigl(\bigl(\tfrac{2\pi}{\sqrt3}+o(1)\bigr)\sqrt{\frac{\log x}{\log\log x}}\Bigr),$$

where $o(1)$ stands for a quantity that tends to $0$ as $x$ grows. This grows more slowly than any power of $x$. ^[7] A search that starts at infinity visits $26\,807$ nodes for $T=10^7$; the C++ program visits $26\,806$, because its starting value $2^{64}-1$ is already below $53\#$, the first primorial with enough solutions there, so rule 2 rejects $53\#$. Its limit is set by its integer types, as described in the previous section, not by the running time.

The C++ program checks itself with its two checkpoints before every normal run, and all three programs print the same number. For this article the following independent checks were made; they are not part of the published solutions.

For every $n\le3000$, trying all $x$ from $n+1$ to $2n$ gives exactly $(d(n^2)+1)/2$ solutions, and every solution found satisfies $(x-n)(y-n)=n^2$. The same direct count confirms that $1260$ is the least $n$ with more than one hundred solutions and that it has $113$ of them.
A sieve, which like the sieve of Eratosthenes runs through the multiples of each prime, stored the smallest prime factor of every $n\le1.2\cdot10^8$; repeated division by it factorised each $n$ and gave $d(n^2)$, without any assumption about the shape of $n$. For every target $T$ from $1$ to $12\,757$, the least $n$ found this way agrees with the output of the C++ program run with --target=T; this covers every target of Section 11 except $T=10^5$.
All $30\,726$ products of primorials up to $47\#$ were listed in two different ways, by lists of exponents and as products of primorials, with the same result, and the least of them with $D\ge8\,000\,001$ is the number printed by the programs; the same list confirms the value for $T=10^5$. A search with the first $30$ primes and exponents up to $200$ gives the same numbers, and for $155$ further targets between $12\,758$ and $36\,167\,512$ it agrees with the C++ program. For the target $36\,167\,513$ the C++ program prints its starting value $2^{64}-1$, as described in the previous section.
An instrumented copy of the C++ program counted the calls and records given above, for $K=5$, $K=201$, four million and $T=10^7$, and an independent reimplementation in Python found the same numbers.

Footnotes and References

The references below supply the background facts used above. The factorisation of the equation, the counting formula, the bounds on primes and exponents, the correctness of the pruned search and all counts of nodes are derived or measured explicitly in this article.

Project Euler 110 — Diophantine Reciprocals II. The official statement: x, y and n are positive integers with 1/x + 1/y = 1/n; for n = 1260 there are 113 distinct solutions, and 1260 is the least n with more than one hundred of them; the task asks for the least n for which the number of distinct solutions exceeds four million. A note calls the problem a much more difficult version of Problem 108 that is well beyond a brute-force approach.
Wikipedia — Optic equation. The optic equation asks for positive integers a, b and c with 1/a + 1/b = 1/c; all its solutions are a = km(m + n), b = kn(m + n) and c = kmn with coprime m and n. The same equation describes two resistors in parallel and the thin lens formula, and c is half the harmonic mean of a and b.
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 a(n) = (τ(n²) + 1)/2, describes a(n) as the number of divisors of n² that are at most n, and states the product formula ((2a₁ + 1)(2a₂ + 1)⋯(2aₜ + 1) + 1)/2 for n = p₁^a₁ ⋯ pₜ^aₜ.
Wikipedia — Fundamental theorem of arithmetic. Every integer greater than 1 is either prime or can be written as a product of primes in a way that is unique up to the order of the factors; the article's example is 1200 = 2⁴ · 3 · 5².
Wikipedia — Divisor function. The function σ₀(n), also written d(n) or τ(n), counts the divisors of n; for n = p₁^a₁ ⋯ pᵣ^aᵣ it equals (a₁ + 1)(a₂ + 1)⋯(aᵣ + 1), and the divisor functions are multiplicative. In the article's example, 24 = 2³ · 3 has 8 divisors.
Wikipedia — Highly composite number. In the prime factorisation of a highly composite number the primes must be the first k primes, since otherwise one of them could be replaced by a smaller prime, and the exponents must be non-increasing, since otherwise exchanging two exponents would give a smaller number with the same number of divisors; such a number is a product of primorials. Ramanujan wrote a paper on these numbers in 1915.
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, 32, 36, 48, … The entry notes that Choie et al. call them Hardy–Ramanujan integers, and that Hardy and Ramanujan proved that there are exp((2π + o(1))/√3 · √(log x / log log x)) of them up to x.
Wikipedia — Primorial. The primorial pₙ# is the product of the first n primes, and n# is the product of the primes up to n; the article's table lists 47# = 614889782588491410 and 53# = 32589158477190044730.
Wikipedia — Branch and bound. Branch and bound splits the set of candidate solutions into a tree of smaller subproblems and discards a subproblem when a bound proves that it cannot contain a better solution than the best one found so far; keeping the open subproblems on a stack gives a depth-first search, which quickly produces complete solutions and therefore bounds.
C++ — Euler110.cpp. The linked immutable C++ revision contains dfs_min_n, solve_u128, run_checkpoints and the options --target and --skip-checkpoints. It is the implementation record for this article; the independent checks described above are not part of it.
GCC manual — 128-bit Integers. GCC provides __int128 and unsigned __int128 for targets with an integer mode wide enough to hold 128 bits, and it cannot express an integer constant of type __int128 on targets whose long long is narrower than 128 bits.
Python — Value comparisons: numbers of the built-in numeric types, such as int and float, can be compared across their types, and they compare mathematically correctly without loss of precision. Java — BigInteger: immutable arbitrary-precision integers, with multiply returning the product and compareTo returning −1, 0 or 1 when the number is less than, equal to or greater than the argument.
cppreference — Arithmetic operators. When a signed integer arithmetic operation overflows, that is, when its result does not fit into the result type, the behaviour is undefined; the page lists wrapping around, trapping and the compiler removing the operation among the possible effects.

Problem 110 source code

C++

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

namespace {

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

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

struct SearchState {
    u128 best_n = static_cast<u128>(~0ULL);
    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;
}

void dfs_min_n(const std::array<int, 16>& 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);
    }
}

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

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

std::string to_string_u128(u128 value) {
    if (value == 0) {
        return "0";
    }

    std::string s;
    while (value > 0) {
        const int digit = static_cast<int>(value % 10);
        s.push_back(static_cast<char>('0' + digit));
        value /= 10;
    }
    std::reverse(s.begin(), s.end());
    return s;
}

bool run_checkpoints() {
    if (to_string_u128(solve_u128(4)) != "6") {
        std::cerr << "Checkpoint failed for target=4" << '\n';
        return false;
    }
    if (to_string_u128(solve_u128(1000)) != "180180") {
        std::cerr << "Checkpoint failed for target=1000" << '\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 << to_string_u128(solve_u128(options.target)) << '\n';
    return 0;
}

Python

# Problem 110: Diophantine reciprocals II
# Find the least n such that 1/x + 1/y = 1/n has more than 4 million solutions.

def solve():
    primes = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53]
    target = 2 * 4000000 + 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

import java.math.BigInteger;

public class Euler110 {
    static int[] primes = { 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53 };
    static BigInteger best = BigInteger.valueOf(Long.MAX_VALUE);

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

    public static void main(String[] args) {
        dfs(0, 63, BigInteger.ONE, 1, 8000001);
        System.out.println(best);
    }
}