Project Euler: #108, Diophantine reciprocals

Thing practiced: math

Tools used: pen, paper, and MonoDevelop

Problem:

In the equation
$\frac{1}{x} + \frac{1}{y} = \frac{1}{n} \tag{1}$ where x, y, and n are positive integers, what is the least value of n where the number of distinct solutions (i.e. solutions where x, y does not recur as y, x) exceeds 1000?

Solution:

Since x, y, and n are positive integers, we know that x and y must each be greater than n. Thus, we can rewrite (1) as
$\frac{1}{n+a} + \frac{1}{n+b} = \frac{1}{n} \tag{2}$ where a and b are again positive integers.

Doing some algebra, we get:
$\frac{n+b}{(n+a)(n+b)} + \frac{n+a}{(n+a)(n+b)} = \frac{1}{n} \tag{3}$ $n(2n+a+b) = (n+a)(n+b) \tag{4}$ $n^2 = ab \tag{5}$

This means that each distinct pair of divisors of n2 represents a distinct solution to (1). Thus, if we can find the number of divisors that n2 has, we’ll know that half that number (well, actually divisors(n2) + 1 / 2 because one of the divisor pairs – n * n = n2 – won’t be doubled) is the number of solutions to (1).

Now we just need an algorithm to find the number of divisors of a given integer m. An easy way to do this would be to enumerate each of the numbers up to sqrt(m) and count each number that divides m evenly as a divisor. Since we’re going to use this algorithm many times, though, this brute force approach won’t work.

Instead, we’ll use the divisor function. This tells us that we can obtain the number of divisors of m from its prime factorization. Specifically, the number of divisors of m is the product over all distinct prime factors of the power of the prime factor, plus 1. (That is, if we can write
$m = \prod^r_{i=1} p_i^{s_i} \tag{6}$ then we have $d(m) = \prod^r_{i=1} (s_i + 1) \tag{7}$ .) Think of it this way: each distinct divisor of m is some distinct combination of m’s prime factors, and each prime factor can appear 0, 1, …, or si times, so multiplying through for each of the prime factors we get (7).)

Now we just need an algorithm to factor a number into primes. We can use a (pseudo-)sieve of Eratosthenes to generate a seed list of potential primes:

We use the list of primes to obtain the prime factorization:

Factoring n2 still takes some time, though, as n2 grows much more quickly than n. Luckily, we can factor n instead. From (6) we can see that $m^2 = (\prod^r_{i=1} p_i^{s_i})^2 = \prod^r_{i=1} p_i^{2{s_i}} \tag{8}$ so with the same combinatorics that gave us (7) we have
$d(m^2) = \prod^r_{i=1}(2{s_i} + 1) \tag{9}$

Now, finally, we can find the solution using (9):

And obtain the solution using an (empirically-determined-to-be) appropriate max prime seed:

This gives us a final answer of 180,180.

Meant to have a crack at #110 today also, but this problem took an embarassingly long time to solve.