Read the details of the problem here

Summary

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

Solution

Well now I see 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 some of the code…and I could, to some extent, but (i) I couldn’t get it running in the time allowed as the assumptions I made regarding the start and the increments weren’t as straightforward to make for this problem so a more exhaustive search of the problem space was required and (ii) because I was doing a more exhaustive search I was hitting situations where the number couldn’t be factored into primes without creating a considerably larger prime list (or having many expensive isPrime() calls).

In other words, as I suspected, I just got lucky with the previous solution.

I did consider changing the divisors() function to deal with the latter case by simply failing & so ruling out the number if it couldn’t factor it with an “appropriate” number of smallish primes but it started to feel a little hackish. I backtracked a little on my reasoning for Problem 108 and got something working that is a little simpler, is more scalable across a larger search space, and is actually much neater in general though it does still rely on “educated” guesses for some of the bounds limiting values.

The solution to Problem 108 found the answer at 22 * 32 * 5 * 7 * 11 * 13, with the factors of the square of the answer being 24 * 34 * 52 * 72 * 112 * 132. The number of divisors was therefore 5*5*3*3*3*3 = 2025.

With the same logic as before we’d be looking for numbers whose square had at least 8 million divisors. Assuming that just squared factors were used that would need a number with up to around 14-15 prime factors as 314 = 4,782,979, 315 = 14,348,907.

I dithered around with these numbers for a while, trying different bases and increments as per my solution to Problem 108 but, as noted above, didn’t get anywhere.

I did some more reading of the forum for the previous problem.

Some people had solved it on paper just by “guessing” the values of the powers of the primes such that the number of solutions would exceed the limit and then working out the actual number from there, checking for closest proximity. This seemed eminently sensible and do-able in code. There was no need to factorize a number; the number of primes being inspected & the allowable powers could be easily tweaked to suit.

Based on the numbers above, the highest prime factor to be considered (assuming that the answer is likely going to be made up of higher powers of the lower primes) is 47. The range of powers I kicked off with would be from 0 (zero – meaning that that prime wasn’t a factor in the answer) to 5, which was fairly arbitrary in that I just doubled what the maximum was for Problem 108 and added 1 for luck…

def ( LIMIT, POWER_RANGE ) = [ 4e6, (0..5) ]
def primes =
    (2..47).findAll { it.toBigInteger().isProbablePrime(100) }

def gen_powers(range, maxlen, yield) {
  gen_powers(range, maxlen, [], yield)
}

def gen_powers(range, maxlen, pows, yield) {
  if (pows.size() == maxlen)
    yield(pows)
  else if (pows.size() == 0)
    range.each { gen_powers(range, maxlen, [it], yield) }
  else
    range.findAll { it >= pows[0] }
         .each { gen_powers(range, maxlen, [it]+pows, yield) }
}

def answer = Long.MAX_VALUE

gen_powers(POWER_RANGE, primes.size()) { powers ->
  if ((1 + powers.inject(1) { p, v -> p * (2 * v + 1) }).intdiv(2) >= LIMIT) {
    def n = 1G
    powers.eachWithIndex { v, i -> n *= primes[i] ** v }
    answer = [ answer, n ].min()
  }
}

This gave the correct answer in 1.13 seconds, so well within my 15 seconds cut off. After having determined what the answer was, and so what the real bounds were, I narrowed the search range as much as possible, and it ran inside 0.06 seconds.

Note: Having gained access into the forum, I wasn’t too surprised to see that some people had done this one with pencil & paper during an idle five minutes…

Conclusion

Groovy was good to do this in, providing listops to allow simple one-liners to be used to keep the code concise and the ability to easily pass a closure into a generator. Performance was adequate for the scale of solution.

It would be nice if the number promotion worked transparently though. I was left scratching my head again when n kept overflowing a Long that I was originally using before I realised I had to specify a BigInteger (the G suffix) instead. I guess it’s a trade-off between performance and ease-of-use but could prove to be a time-bomb in inadequately tested code.

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.

Project Euler Problem #102

November 20, 2009

Read the details of the problem here

Summary

For how many triangles in the text file does the interior contain the origin?

Solution

This question had me wondering for a while as to what the best approach to take might be. In the end I went for a simple approach that just relied on the fact that the interior angles from any point within a triangle to the vertexes would add up to 360° or 2 * Pi radians. To calculate the angles at the origin I relied on some things I’d learned when I used to write computer games as a hobby.

The dot product of two 2D vectors, U and V, is equivalent to the following two formulae:

U•V = U0V0 + U1V1

U•V = |U| * |V| * cos θ

Where |U| is the magnitude of the vector U, i.e. sqrt( U0² + U1² )

So this gives: cos θ = ( U0V0 + U1V1 ) / ( |U| * |V| )

The coded solution is then quite trivial and just calculates the angle at the origin subtended by vectors out to the pairs of vertices for each line segment.

def answer = 0

new File("triangles.txt").eachLine {
    def (x0, y0, x1, y1, x2, y2) = it.split(",")*.toInteger()
    def ( a, b, c ) = [ [ x0, y0 ], [ x1, y1 ], [ x2, y2 ] ]
    def theta = 0.0

    [ [ a, b ], [ a, c ], [ b, c ] ].each {
        def ( u, v ) = it
        def dot = u[0] * v[0] + u[1] * v[1]
        def ul =  Math.sqrt(u[0] * u[0] + u[1] * u[1])
        def vl =  Math.sqrt(v[0] * v[0] + v[1] * v[1])
        theta += Math.acos(dot / (ul * vl))
    }
    
    if (Math.abs(theta - 2 * Math.PI) < 1E-7) answer++
}

It ran in 0.96 seconds and wasn’t difficult to write once I’d decided what approach to take! Note the importance of using an epsilon value (1E-7) when comparing floating-point numbers as you can’t rely on an exact equality.

There are probably quicker ways to do this (in fact, I know there are) but this approach was nice & straightforward to code off the top of my head and ran in an acceptable time.

Conclusion

Groovy makes reading a text file and parsing it into appropriate data structures quite simple. It didn’t really make any difference to the rest of the solution apart from perhaps the nice syntax on the <list of sides>.each{} loop as it was mostly just straightforward maths.