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

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.

$\displaystyle\frac{1}{x}+\frac{1}{y} =\frac{1}{n}$

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 $p_1^{a_1} \cdot p_2^{a_2} \cdot p_3^{a_3} \dots$ sorted such that $p_1 < p_2  the this is a candidate only if $a_1 \geq a_2 \geq a_3 \geq \dots$.

The reason is that if we have we have two exponents $a_i < a_j$ where $i < j$ then we can make a number which is smaller by interchanging the exponents. If we write this out as a product consisting of $a_ia_j$ numbers, interchanging the exponents will still give us $a_ia_j$ factors $a_i - a_j$ of these will be smaller when interchanging the exponents. Therefore lowering the product. One example is that we have $2^2 3^3 = 108$. Then we can make this smaller number by interchanging the exponents such that we get $2^3 3^2 =72$ 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.

### Posted by Kristian

David

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)

Erick

In the beginning, we calculate a result which is a product of the first 14 primes. I believe this should be 15 primes: 3^15 = 14348907, 3^14 = 4782969. And since we need at least 7999999 factors, 3^14 would not be sufficient.

This is important because I want to suggest using this result to limit the number of candidates for our brute force search. Because we know this result produces the required number of solutions, we don’t want any number bigger than this and we partially eliminate them using the following steps:

1) We calculate (int) Math.Log(result,2). This number is the highest power of 2 for which the resulting number would be below the result.

In this case it would be 59. The implications of this is as follows:
– If we do 2^60, we would reject it anyway, so there is no point in doing this
– If we do 2^59 only, we will get only 119 factors which isn’t sufficient.
– Therefore, we don’t need to search a combination with only powers of 2

This in itself is a small finding. So we repeat it with the next prime.

2) The highest power of 3 that we can use is 22. We get this by doing (int) Math.Log(result,6), with 6 being 2 * 3 given that lower prime numbers can’t have lower power.

From this, we get the following deduction:
– There is no point to do 3^23 or any higher than that
– While we can prove that 3 cannot go above 22, we cannot prove the same for 2. But we still have a limit for powers of 2 from number 1)
– If we use 2^59 * 3^22, we know that we would get 119 * 45 = 5355 factors which isn’t enough
– And we can deduce that if we increase the power of either 2 or 3, we can’t use the resulting number
The logic goes like this: if 2^60 * 1 is too big 2^60 * x will also be too big if x > 1.
Therefore, if we write x = 3^y where y >= 0, it doesn’t matter what y is, we can’t use any power of 2 bigger than 59.
– Therefore, we cannot use only a combination of 2 and 3.

3) If we then repeat this process, we would find that we need at least 5 prime numbers: 2,3,5,7,11. As a bonus, we can also deduce that we should start the brute force search with 11^2. This is a significant reduction in possibilities.

So I’ve modified your code to include this logic. Between the loop where result is calculated and the main while loop, I have:

long limit = 2*4000000-1;
int counter = 0;
double remainder = limit;
int log_base = 1;
int max_power = 0;
while (max_power < remainder) {
remainder = remainder / (max_power * 2 + 1);
log_base = log_base * (int) primes[counter];
max_power = (int) BigInteger.Log(result,log_base);
counter++;
}
counter–;
exponents[counter] = (int) Math.Ceiling((decimal)(remainder-1)/2);
SetAllSmallerExponents(counter);

I get a time of around 78-79 ms
(my machine is slower than yours)

Antonio Sanchez

Minor correction: we need at most 15 primes, since we really need 2*4000000-1 divisors of N^2.

The simplest approach is to start building up exponents: 2, 2*3, 2*3*4, …, 2^2*3, 2^2*3^2…

We can quit our search whenever our current N is greater than the product of the first 15 primes = 614889782588491410. We also know that the highest possible exponent is 60, since log2(614889782588491410) = 59.1.

I took a recursive approach, incrementing exponents in order and terminating when we either surpass the current best solution, or find the next best solution (since multiplying by any more factors will result in a worse solution than the one we just found).

Runs < 1 ms.

typedef unsigned long long ULL;
typedef ULL value_type; // or big-integer type

void recursive_build(ULL dmin, ULL d, vector<ULL>& primes, vector<size_t>& exponents, size_t k, value_type& x, value_type& xmin) {

value_type y = x;
for (size_t i=1; i<= exponents[k-1]; ++i) {
exponents[k] = i;
ULL dnext = d*(2*i+1);
y = y * primes[k];

// quit if bigger than current best
if (y > xmin) {
return;
}

// quit and save if we have new best
if (dnext > dmin) {
xmin = y;
return;
}
recursive_build(dmin, dnext, primes, exponents, k+1, y, xmin);
}
}

void euler(ULL N) {

// we need X such that X^2 has at least this many divisors
ULL dmin = 2*N-1;

// count how many unique primes until we get # divisors (provides a bound on X)
size_t nprimes = 1;
ULL p = 3;
while (p < dmin) {
p = p*3;
++nprimes;
}
vector<ULL> primes = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47};
// gen_primes<ULL>(nprimes); // or generate primes if more required

// current best estimate is first k primes
value_type xmin = 1;
for (size_t i=0; i<nprimes; ++i) {
xmin = xmin * primes[i];
}

// recursively build up number by slowly increasing exponents
vector<size_t> exponents(nprimes, 0);
value_type y = 1;
recursive_build(dmin, 1, primes, exponents, 0, y, xmin);

// solution
cout << xmin << std::endl;

}