Solution

Answer: 9350130049860600

We seek the least positive integer $n$ such that the number of distinct positive integer solutions of

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

exceeds four million.

The key to the problem is transforming the Diophantine equation into a divisor-counting problem.

Step 1 — Transform the equation

Start with

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

Multiply through by $xyn$:

$$ny + nx = xy$$

Rearrange:

$$xy - nx - ny = 0$$

Add $n^2$ to both sides:

$$xy - nx - ny + n^2 = n^2$$

Factor:

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

So every positive divisor pair $(d,e)$ with

$$de=n^2$$

corresponds to a solution

$$x=n+d,\qquad y=n+e$$

Thus the number of ordered solutions equals the number of divisors of $n^2$:

$$\tau(n^2)$$

where $\tau$ is the divisor-counting function.

But the problem asks for distinct solutions, meaning $(x,y)$ and $(y,x)$ are considered the same. Since divisor pairs are symmetric, the number of distinct solutions is

$$\frac{\tau(n^2)+1}{2}$$

(the $+1$ accounts for the square root divisor when $d=e=n$).

We need

$$\frac{\tau(n^2)+1}{2} > 4{,}000{,}000$$

Therefore

$$\tau(n^2) > 7{,}999{,}999$$

So we need

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

Step 2 — Express $\tau(n^2)$ from prime exponents

Suppose

$$n=\prod p_i^{a_i}$$

Then

$$n^2=\prod p_i^{2a_i}$$

and therefore

$$\tau(n^2)=\prod (2a_i+1)$$

So the problem becomes:

Find the smallest integer

$$n=\prod p_i^{a_i}$$

such that

$$\prod (2a_i+1) > 8{,}000{,}000$$

Step 3 — Minimizing $n$

To minimize $n$:

use the smallest primes first,
and enforce nonincreasing exponents:

$$a_1\ge a_2\ge a_3\ge \cdots$$

Otherwise swapping exponents onto smaller primes would reduce $n$.

This becomes a recursive search / branch-and-bound optimization problem.

Step 4 — Python implementation

from math import prod

TARGET = 8_000_000

# Small primes are always optimal
primes = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41]

best = float('inf')

def search(idx, max_exp, current_n, divisor_count):
    """
    idx            : index into primes
    max_exp        : largest exponent allowed here
    current_n      : current value of n
    divisor_count  : current value of tau(n^2)
    """
    global best

    # If enough divisors have been achieved, update best
    if divisor_count > TARGET:
        best = min(best, current_n)
        return

    # No point continuing if already worse
    if current_n >= best:
        return

    p = primes[idx]
    value = current_n

    # Try exponents in decreasing-valid order
    for exp in range(1, max_exp + 1):
        value *= p

        if value >= best:
            break

        # multiplying divisor count by (2*exp + 1)
        new_divisors = divisor_count * (2 * exp + 1)

        search(idx + 1, exp, value, new_divisors)

search(0, 60, 1, 1)

print(best)

Step 5 — Code walkthrough

We recursively build

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

while tracking:

$$\tau(n^2)=\prod (2a_i+1)$$

Important details

Exponent ordering

search(idx + 1, exp, ...)

ensures later exponents cannot exceed earlier ones.

This guarantees minimality.

Branch-and-bound pruning

if current_n >= best:
    return

If we already exceeded the best candidate, stop exploring.

This is essential for performance.

Divisor-count update

If we add exponent $e$, then

$$(2e+1)$$

multiplies the divisor count.

Hence:

new_divisors = divisor_count * (2 * exp + 1)

Step 6 — Verify against the example

For

$$n = 1260 = 2^2\cdot3^2\cdot5\cdot7$$

we get

$$\tau(n^2) = (5)(5)(3)(3)=225$$

Therefore the number of distinct solutions is

$$\frac{225+1}{2}=113$$

which matches the problem statement exactly.

Step 7 — Final computation

Running the search yields:

$$n = 9350130049860600$$

Check: