Project Euler 39

April 12th 2008 06:49 pm

If p is the perimeter of a right angle triangle with integral length sides, {a,b,c}, there are exactly three solutions for p = 120.
{20,48,52}, {24,45,51}, {30,40,50}
For which value of p < 1000, is the number of solutions maximised?

WARNING: CONTAINS MATHEMATICS

You may remember from school



(Difference of 2 squares)
let

and

so


and substituting these into
gives

Set and to get rid of that root and you have


Primitive triples are ones in lowest terms (if you know (3,4,5) you can get (6,8,10) etc by multiplying each side by some integer)
If p and q are coprime and different parity you just get primitive ones (i did not bother and just used set() to remove dupes from the list)

Here is my code

<br />
from __future__ import division<br />
from collections import defaultdict<br />
from time import time</p>
<p>start = time()<br />
counter = defaultdict(list)</p>
<p>def triples():<br />
    for q in range(1, 35):<br />
        #35 is greater than sqrt(1000)<br />
        for p in range(q+1, 35):<br />
            a = p**2 - q**2<br />
            b = 2 * p * q<br />
            c = p**2 + q**2<br />
            for m in range(1,int(1000/c)):<br />
                #generate multiples till c is above 1000<br />
                #triples with perimeter > 1000 are filtered later<br />
                yield (m*min(a,b),m*max(a,b), m*c)</p>
<p>for triple in triples():<br />
    perimeter = sum(triple)<br />
    if perimeter < 1001:<br />
        counter[perimeter] += [triple]</p>
<p>print set((counter[120]))</p>
<p>maxlen = 1<br />
for i in counter:<br />
    length = len(set(counter[i]))<br />
    if length > maxlen:<br />
        maxlen = length<br />
        print i<br />
end = time()</p>
<p>print "took ", end-start</p>
<p>

Which outputs:
set([(24, 45, 51), (20, 48, 52), (30, 40, 50)])
520
528
630
660
720
840
took 0.0140960216522

I tried brute forcing it today and knocked up the following

<br />
from collections import defaultdict<br />
counter = defaultdict(int)<br />
from time import time</p>
<p>def run():<br />
	start = time()<br />
	for perimeter in range(1,1001):<br />
		#print perimeter<br />
		for a in range(1,perimeter):<br />
			for b in range(1,perimeter-a):<br />
				c = perimeter - a - b<br />
				if a**2 + b**2 == c**2:<br />
					counter[a+b+c] += 1<br />
	mosttriples = 0<br />
	for i in counter:<br />
		if counter[i] > mosttriples:<br />
			mosttriples = counter[i]<br />
			print i<br />
	end = time()<br />
	print "TOOK", end-start</p>
<p>print "without psyco"<br />
run()</p>
<p>print "with"<br />
import psyco<br />
psyco.full()<br />
run()<br />

Which gives
without psyco
516
520
528
630
660
720
840
TOOK 134.663383007
with
516
520
528
630
660
720
840
TOOK 21.0988409519

So psyco makes it 6 times faster but a smarter algorithm is 2000 times faster still.

Posted by tom under euler & Mathematics & Python |

2 Responses to “Project Euler 39”

  1. ejs responded on 14 Apr 2008 at 11:18 am #

    While I don’t have the code I used I remember solving it far faster using brute force. After 10 min playing this is what I think I had (hoping the indentation works):

    import math
    from eulerhelpers import timed

    @timed
    def solve39(cap=1000):
    result = [0] * (cap+1)
    for b in xrange(1, cap/2):
    for c in xrange(1, min(b, cap/2-b)+1):
    a = math.sqrt(b**2 + c**2)
    if a+b+c > cap:
    break
    if not a%1:
    result[int(a)+b+c] += 1
    if not b%(cap/20):
    print b,
    yield
    print max((i, e) for e, i in enumerate(result))

    On my computer, without psyco, it completes in just over 0.1 seconds or about 7 times longer than your more mathematical code. I stopped there thinking it was enough.

  2. tom responded on 14 Apr 2008 at 1:32 pm #

    http://paste.bradleygill.com/index.php?paste_id=238

Trackback URI | Comments RSS

Leave a Reply