Project Euler Problem #108

June 5, 2010

Read the details of the problem here

Summary

Solving the Diophantine equation 1/x + 1/y = 1/n.

Solution

For a while I’ve been seeing that people have been finding this blog through searches for solutions to Project Euler 110 so I thought I’d give it a go so that I could actually satisfy those readers. However, when I went to solve that I found that this one was strongly recommended to be done first. So, here I am!

As Groovy is generally very slow in brute-force approaches if there’s a lot of maths processing going on, as I suspect there would likely be for such a solution as this (using nested loops for values of x and y and so expecting an O(n2) performance), it’s worth trying to analyse the problem to reduce the problem scope down somewhat.

Try and find y in terms of x and n in order to save having to run two loops.

Given: 1/x + 1/y = 1/n

    n/x + n/y = 1
    n + nx/y = x
    ny + nx = xy
    nx = xy – ny
    nx/y = x – n 
    1/y = (x-n)/nx 
    y = nx/(x-n)

Substitute: d = x – n
Giving: y = nx/d

As: d = x – n then x = d + n
Substitute back into y = nx/d 
Gives: y = n(d+n)/d
 
    y = (nd + n2)/d
    y = n + n2/d

This means that, for y to be an integer, then d (i.e. xn) must be a divisor of n2. We’ve also got obvious bounds on the value of x which are evident from the question itself and from the formulae above; n < x  ≤ 2n.

Plugging this back into the exemplar where n = 4. The divisors of 16 (n2) being 1, 2, 4, 8 & 16.

divisor (d) x = (d + n) y = (n + n2/d) n < x ≤ 2n
1 5 20 True
2 6 12 True
4 8 8 True
8 12 6 False
16 20 5 False

So from the 5 solutions to the equation, only 3 are permissible given the constraints on the value of x.

Note the symmetry for the values of x and y around the value of n. A little bit of experimentation shows that the number of solutions for any particular n is simply: (divisors(n2) + 1)  / 2.

So now there’s no need to iterate over x and y looking for solutions, it’s just a matter of counting the number of divisors for the square of any value of n looking for at least 2000 divisors i.e. > (1000*2)-1 – unfortunately as the formula uses n2 those numbers get really large, really fast, which doesn’t bode well for a Grooovy factorisation solution. We’re still looking at O(n2) complexity.

As per various other Project Euler questions, the number of divisors is obviously going to be high when the square of the number is a product of a large number of small primes. So that seems to be a good place to start.

As learned in earlier solutions, the number of divisors of an integer are the product of the powers of each prime factor plus one. So, for instance, looking at 16 which is 24 – the power is 4 so there are 5 divisors. Using 24 as another example, 24 is 23 * 31 – so there are (3+1) *(1+1) = 8 divisors ( 1, 2, 3, 4, 6, 8, 12, 24 ).

Given that we’re trying to find a square number (i.e. n2), and that (ab)2 = a2 * b2, it means that the powers of the factors are going to be even and so therefore at least 2, giving minimum element multipliers of 2+1 = 3.

Assuming that all the factors are just squares of the primes what would be lowest number be? 36 =  729,  37 = 2187. So we’re looking for an answer for n2 somewhere between 22 * 32 * 52 * 72 * 112 * 132 = 300302  and 22 * 32 * 52 * 72 * 112 * 132 * 172 = 5105102 with the upper bound being inclusive as a potential solution.

This is a substantially smaller search space than before, and I’ll make the assumption that it’s going to be a multiple of the lower bound, as this is where maxima are likely to be seen, which means that it’s (probably) only going to be necessary to check 16 values of n at most. Odds are it’s going to be higher powers of the smaller primes in the solution but these assumptions make a brute-force approach viable.

The program is therefore quite straightforward, with the divisor method using prime factorisation (as above) to calculate the number of divisors.

def divisors(num, primes) {

  def (n, count, idx) = [ num, 1, 0 ]

  while (n > 1) {
    def base = primes[idx]
    while (n % base != 0) {
      base = primes[++idx]
    }

    def exp = 1
    while ((n > 0) && (n % base == 0)) {
      n = n.intdiv(base)
      exp++
    }
    idx++
    count *= exp
  }
  count
}

def primes = (2..17).findAll { it.toBigInteger().isProbablePrime(100) }
def inc = primes[0..<6].inject(1) { p, v -> p * v }
def n = inc + inc
def (answer, LIMIT) = [ 0, 1999 ]

while (answer == 0) {
  if (divisors(n ** 2, primes) > LIMIT) answer = n
  n += inc
}

This ran in 0.10 seconds. Well within the 15s permitted time. I get the feeling that I may have got lucky with this one though – my start number and increment were somewhat arbitrarily chosen. I’ll have to see if this logic stands up to scrutiny for Problem 110 which looks a lot harder with at least 4 million solutions being sought.

Note: It took me a LOT longer to work this out than it did to write about it!

Conclusion

One of those engage brain before typing questions. Although I used a few Groovy listop features and took advantage of the terser syntax, it wasn’t really a Groovy solution. After the analysis to reduce the search scope performance was entirely adequate.

Advertisements

6 Responses to “Project Euler Problem #108”


  1. […] why this is one of the most popular search terms through which people find this blog. Having done Problem 108 which was recommended to do before attempting this one, I was hoping that I could just repurpose […]


  2. […] blog to see how I could code this up more elegantly, but after googling around for a bit I found a solution which puts mine to shame. His (or her) solution runs an order of magnitude faster, despite using a […]

  3. Anders Says:

    “A little bit of experimentation shows that…”

    Dude, not cool. This holds true, as is easily realized, but mathematical truths are not verified by checking three cases, five, or a thousand. Either we do math or we can climb right back up in the trees.

    Otherwise, good blog. 🙂

    • keyzero Says:

      True. I do realise this. I am just a hack after instant gratification. Got stung badly on trying to scale this to Problem 110 too.


  4. […] Problem #108 and Problem #110 also tackled the domain of products of powers of primes and the solution to the latter seemed a good place to start as the one to the former ended up as being a bit of a lucky hack. As per that solution, taking the approach of assuming that the answer was to be found somewhere at n * product(prime factors) and driving off that meant that factorisation of large numbers wasn’t needed so Groovy would remain a viable language to do this in. […]


  5. […] KeyZero Conversation: Comments on Computing […]


Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

%d bloggers like this: