Project Euler 110: Find an efficient algorithm to analyse the number of solutions of the equation 1/x + 1/y = 1/n.

Project Euler 110: Find an efficient algorithm to analyse the number of solutions of the equation 1/x + 1/y = 1/n.

Written by on 18 August 2012

Topics: Project Euler

Problem 110 of Project Euler is a continuation of Problem 108 except that we need a solution approach which is much more efficient. The problem reads

In the following equation xy, and n are positive integers.

It can be verified that when n = 1260 there are 113 distinct solutions and this is the least value of n for which the total number of distinct solutions exceeds one hundred.

What is the least value of n for which the number of distinct solutions exceeds four million?

NOTE: This problem is a much more difficult version of problem 108 and as it is well beyond the limitations of a brute force approach it requires a clever implementation.

I wont go through all the basics of this which is covered in the solution of Problem 108, so go read that before going on, I will wait….

Okay done? So what happens if we just run the same code? Well first of all we would get in problems with the capacity of the variables, but we could just change it to BigInteger and then run it. I tried, and gave up after a while. So I changed the solution approach.

We know that we can use the prime factorization of a number to calculate the total number of divisors. So why not build numbers from prime factors? We only need a few rules in order to pursue this approach.

The first observation we need is that if we have a prime factorization sorted such that  the this is a candidate only if .

The reason is that if we have we have two exponents where then we can make a number which is smaller by interchanging the exponents. If we write this out as a product consisting of numbers, interchanging the exponents will still give us factors  of these will be smaller when interchanging the exponents. Therefore lowering the product. One example is that we have . Then we can make this smaller number by interchanging the exponents such that we get  while we have the same number of divisors.

The second thing we need is to realise that we need at most 14 distinct prime factors since 3^14 = 4782969.

The third thing I did was realizing that we don’t need to loop through all possibilities, if we have the exponents for all the prime factors except the exponent for 2, we can just calculate the size of it based on the limit.

I did this calculation this way

public void Twos(long limit) {
    exponents[0] = 0;
    long divisors = 1;

    for (int i = 0; i < exponents.Length; i++) {
        divisors *= 2 * exponents[i] + 1;
    }
exponents[0] = (int)(limit / divisors - 1) / 2;

while (divisors*(2*exponents[0]+1) < limit)
exponents[0]++;
}

The last while loop I did because there are several ways that we can get some rounding errors, so I found this way the easiest to implement. I am sure that you can analyse the algorithm and find out exactly how to do it.

once we have a set of exponents we can calculate the resulting number as

public BigInteger Number(long[] primes){
    BigInteger number = 1;

    for(int i = 0; i < exponents.Length;i++){
        number *= BigInteger.Pow(primes[i], exponents[i]);
    }
    return number;
}

The only thing missing now is a complicated while loop to increment the exponents to give us the right candidates. Every time we have a candidate number where the number of twos are lower than the number of threes we need to increment the exponent of a larger prime factor and then set all smaller prime factors’s exponents to the same number. This is done with the following code

long[] primes = ESieve(2,45);
exponents = new int[primes.Length];

//Find the absolute max we need
BigInteger result = 1;
for(int i = 0; i < primes.Length; i++){
    result *= primes[i];
}

long limit = 2*4000000-1;
int counter = 1;

while (true) {
    Twos(limit);

    if (exponents[0] < exponents[1]) {
        //invalid combination
        counter++;
    } else {
        //The number of twos are larger than
        // number of threes, so good solution
        BigInteger number = Number(primes);
        if (number < result)
            result = number;
        counter = 1;
    }
    if (counter >= exponents.Length)
        break;

    exponents[counter]++;
    SetAllSmallerExponents(counter);
}

This gives us

The first n with more than 4000000 solutions is 9350130049860600

which runs in

Solution took 305,4984 ms

Optimizing the code

Working with BigInteger is always a performance killer, so I thought I would hunt for a bit of performance there. When we calculate the number we can break on two instances. If the exponents start becoming 0 we don’t really need to continue but can break the for loop.

The other instance is that if the number becomes larger than the current best result. So I changed the Number algorithm to

public BigInteger Number(long[] primes, BigInteger result){
    BigInteger number = 1;

    for(int i = 0; i < exponents.Length;i++){
        if (exponents[i] == 0) break;
        number *= BigInteger.Pow(primes[i], exponents[i]);

        if(result < number)
            return result+1;
    }

    return number;
}

This gives us

The first n with more than 4000000 solutions is 9350130049860600

 

Solution took 105,6026 ms

which is a decent speed up. I am uncertain if we can do more to speed this up.

Wrapping up

We managed to change the solution approach from Problem 108 to be fast enough to solve it for a larger number. If we apply this solution to the problem in 108, we solve it in less than 1ms. Which is a nice speedup.

The code for this problem can be found here.

1 Comment For This Post I'd Love to Hear Yours!

  1. David says:

    You can speed up the code even more by not using BigInteger at all. I coded up your solution in FORTH and it runs instantaneously by doing everything with 64 bit math.

    The trick is in Number(), stop the computation as soon as you exceed 14 primorial, which is the upper bound you determine at the start of the main loop. Since you don’t care about values greater than 14#, there’s no reason to perform the full computation.

    I don’t think you can use exponentiation by squaring, but it’s not important since Number() would do at most 54 multiplications (log2(14#) = 53)

Leave a Comment Here's Your Chance to Be Heard!

You can use these html tags: <a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <strike> <strong>

You can use short tags like: [code language="language"][/code]

The code tag supports the following languages: as3, bash, coldfusion, c, cpp, csharp, css, delphi, diff,erlang, groovy, javascript, java, javafx, perl, php, plain, powershell, python, ruby, scala, sql, tex, vb, xml

You can use [latex]your formula[/latex] if you want to format something using latex

Notify me of comments via e-mail. You can also subscribe without commenting.