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 n² represents a distinct solution to (1). Thus, if we can find the number of divisors that n² has, we’ll know that half that number (well, actually divisors(n²) + 1 / 2 because one of the divisor pairs – n * n = n² – 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 s_i 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:

// This is not optimized, and not the real sieve
// -- see Melissa E. O'Neill, The Genuine Sieve
// of Eratosthenes -- because it doesn't matter
// here and my brain hurts.

// Returns a list of primes up to a given maximum.
let primes max =
    let rec sieve list =
        match list with
        | [] -> []
        | hd :: rest ->
            hd :: (sieve <| List.filter (fun x -> x % hd <> 0) rest)
    sieve [ 2 .. max ]

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

// Returns the prime factorization of a number
// as a list of tuples of (prime, exponent). This
// only looks for prime factors up to a given
// maximum. Again this isn't optimized, because
// it doesn't matter for this size of problem.
// Definitely won't work for #110 though.
let primeFactors m max =
    let rec primeFactorsHelper (m, factorsSoFar) primesToTest =
        let singleFactorHelper m p factorsSoFar =
            let rec exponentHelper m p e =
                if m % p = 0 then exponentHelper (m/p) p (e+1)
                else (m, e)
            let (rem, exp) = exponentHelper m p 0
            rem, (p, exp) :: factorsSoFar
        match primesToTest with
        | [] -> factorsSoFar
        | hd :: rest when hd > m -> factorsSoFar
        | hd :: rest when m % hd = 0 ->
            primeFactorsHelper (singleFactorHelper m hd factorsSoFar) rest
        | hd :: rest -> primeFactorsHelper (m, factorsSoFar) rest
    primeFactorsHelper (m, []) (primes max)

Factoring n² still takes some time, though, as n² 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):

// Returns the number of solutions for a given n.
// This still only uses prime factors up to a
// given maximum.
let solutions m maxPrime =
    primeFactors m maxPrime
    |> List.fold (fun acc (p, e) -> acc * ((2*e) + 1)) 1

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

// Returns the solution, finally, I think.
let problem108answer =
    { 1 .. System.Int32.MaxValue }
    |> Seq.find (fun n -> (((solutions n 17) + 1) / 2) > 1000)

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.