# Project Euler 108: Solving the Diophantine equation 1/x + 1/y = 1/n.

Let us jump right into Problem 108 of Project Euler which reads

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

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

For n = 4 there are exactly three distinct solutions:

$\displaystyle\frac{1}{5}+\frac{1}{20} =\frac{1}{n}$

$\displaystyle\frac{1}{6}+\frac{1}{12} =\frac{1}{n}$

$\displaystyle\frac{1}{8}+\frac{1}{8} =\frac{1}{n}$

What is the least value of n for which the number of distinct solutions exceeds one-thousand?

NOTE: This problem is an easier version of problem 110; it is strongly advised that you solve this one first.

We could possibly brute force, but my guess is that it would take a rather long time to do so. So instead I started fiddling around with the equation to get some insight into the problem.

Starting with the original

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

we can see that both x and y needs to be greater than n since they both need to be non-negative. This means that we can express the equation as

$\displaystyle\frac{1}{n+r}+\frac{1}{n+s} =\frac{1}{n}$

where $r,s \in \mathbb{N}$, meaning that they are integers strictly larger than zero.

Putting this on a common denominator gives us

$\displaystyle\frac{n+s+n+r}{(n+r)(n+s)} =\frac{1}{n} \Leftrightarrow$

$\displaystyle n(n+s+n+r) =(n+r)(n+s) \Leftrightarrow$

$\displaystyle n(n+s+n+r) =n^2 + nr + ns + sr \Leftrightarrow$

$\displaystyle n^2 = sr \Leftrightarrow$

So both s and r are divisors of n^2, which also means that all divisor pairs of $n^2$ will give us a solution to the problem. However, only half of the divisor pairs will give a unique solution since you can simply swap r and s and get a solution as well.

So if d(n) is a function giving us the number of divisors of n, we are looking for the first n such that d(n^2) / 2 > 1000.  All numbers have an even number of divisors unless $\sqrt{n}$ is a divisor of n. In that case the number will have an odd number of divisors. I wont prove this, but I will leave it for you to figure out. In our case we have $\sqrt{n^2} = n$ which is also a divisor, therefore we are looking for a solution to which have (d(n^2) + 1) / 2 > 1000, where the division is an integer division.

## Finding divisors

The huge task in this exercise is to actually find the smallest number with many divisors. In Problem 12 we were asked to find the number of divisors of a number. Where we established that we could use the prime factorization to find the number of divisors such that if we have the number N we have

$\displaystyle N = p_1^{a_1} \cdot p_2^{a_2} \cdot p_3^{a_3} \dots = \prod_{n=1}^K p_n^{a_n}$

and then the number of divisors will be

$\displaystyle d(N) =(a_1+1) \cdot (a_2 +1) \cdot (a_3+1)\dots= \prod_{n=1}^K (a_n+1)$

However, we can cheat a bit so instead of factoring $n^2$ we will factor $n$ and the number of factors of  $n^2$ will then be

$\displaystyle d(N) =(2a_1+1) \cdot (2a_2 +1) \cdot (2a_3+1)\dots= \prod_{n=1}^K (2a_n+1)$

since we have that

$\displaystyle N^2 =( p_1^{a_1} \cdot p_2^{a_2} \cdot p_3^{a_3}\dots)^2 = (\prod_{n=1}^K p_n^{a_n})^2 =\prod_{n=1}^K p_n^{2a_n}$

This will help us significantly since it is much easier to factor a smaller number.

The algorithm I use for prime factorization requires a list of prime numbers, so I started probing around to figure out what I needed. I found that (2357111317)^2 has 2187 divisors, which is more than the 1999 we need in order to solve the problem, so we can set the upperlimit of the prime list to 17, and if a number have a prime factor larger than 17, that is not the number we are looking for anyway.

## Coding the solution

I have changed the divisor counting function as described above, which resulted in  the C# code

private long NoDSquared(long number, long[] primelist) {
long nod = 1;
long exponent;
long remain = number;

for (int i = 0; i < primelist.Length; i++) {
// In case there is a remainder this is a prime factor as well
// The exponent of that factor is 1
if (primelist&#91;i&#93; * primelist&#91;i&#93; > number) {
return nod * 2;
}

exponent = 1;
while (remain % primelist[i] == 0) {
exponent+=2;
remain = remain / primelist[i];
}
nod *= exponent;

//If there is no remainder, return the count
if (remain == 1) {
return nod;
}
}
return nod;
}


and then the main loop for the code is

long n = 1;
long result = 0;
long limit = 1000;

while (true) {
if ((NoDSquared(n, primes) + 1) / 2 > limit) {
result = n;
break;
}
n++;
}


Nothing more to say than we should run the code which on my computer gives me the result

The first n with more than 1000 solutions is 180180


and the time

Solution took 31,0487 ms


## Wrapping up

While this most certainly solved problem 108 within the time limit, it is far from fast enough for Problem 110, so we will be revisiting this problem soon enough.

You can find the complete source code here as always.

### Posted by Kristian

Eman

Hello, nice solution I solved this on the third or fourth try. I used this approach: I found the the number we are looking for must be very composite (lot of factors) so I made some test, some tries and finally found solution.

P.S this is very stupid factorization metod but it was from my beginings

#include
#include
#include

using namespace std;

const int MAX_SOLUTIONS = 1000;
typedef unsigned long long Long;
int getSolutions(Long n);

int main()
{

time_t start = clock();

int solutions = 0;
Long n=2*3*5*7*11*13;
Long x = 2*3*5*7*11*13;

while(solutions <= MAX_SOLUTIONS)
{

solutions = getSolutions(n);
//cout<<n<<'\t'<<solutions<<endl;
//cin.get();
n+=x;
}

cout<< "Project Euler 108:\t"<<n - x<<"\t"<<solutions<<endl;
cou<<"PROGRAM LAST:\t"<<(clock() - start)"<< ms."<<endl;

}

int getSolutions(Long  n)
{

Long maxX = 2*n;
register int solutions=0;
Long differ;
long double y;
Long x;

for(x = n+1; x<=maxX; x++)
{
differ = x-n;

y = (x*n)/(long double)differ;
//if(y < x)
// return solutions;
if((Long)y == y)
{
solutions++;
// cout<<x<<"\t"<<y<<"\t"<<differ<<endl;

}

}

return solutions;

}


robert

Your pages show very well with the MathML. I wonder whether you ca let me know what sort of software that you are using to get the Math symbols right. Thanks.

Kristian

Hi Robert

I use wordpress as you probably have figured out and to render the math I use the optimized latex plugin

Ralf

How would you solve Diophantine equation a*c1+b*c2+c*c3+d*c4+…+i*c9=n?

Manas Srivastava

Consider the equation: 1/x + 1/y= 1/6
Acc. to you,

6=2^1 * 3^1
So, No. of distinct possible are= (1+1) * (1+1)
= 4

But following solution are possible: (7,42) (8,24) (9,18) (10,15) (12,12)
These are 5 solutions.
Please clear my doubt and correct me where I am wrong.

Individ
QuasiChameleon

@Manas Srivastava

According to Kristian, the number of distinct possibilities is floor((d(n^2) + 1) / 2), which for 6 becomes floor([(2*a2 + 1)*(2*a3 + 1) + 1]/2) = 5.

Note that a2 and a3 are both 1.

BartV

2187 is simply 3^7 , the number of divisors of any number that’s a product of 7 different prime factors, following the formula you quoted.

I think you were well on the way to solving this problem without any code whatsoever.

2*3*5*7*11*13*17 > 2^4 *3*5*7*11*13*17 > 2^2 * 3^2 * 5*7*11*13*17
3^7 == 3^5 * (2*4+1) > 3^4 * (2*2+1) * (2*2+1) > 1999

_tr_#

when i find analysis’s Kristan i try it in hackkerank but i only earn 50 d 🙂 , i used miller rabin but still wrong with testcacse
5,7,8,9,11.
https://www.hackerrank.com/contests/projecteuler/challenges/euler108/problem?h_r=profile
I think this analysis’s Kristan “n^2=rs ” have relate to analysis’s pytagor
r
s=[(r+s)/2]^2 – [(r-s)/2]^2=n^2
-> (2*n)^2+(r-s)^2=(r+s)^2 , and i do not know it help everything lighting up ?