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.
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.
tom responded on 14 Apr 2008 at 1:32 pm #
http://paste.bradleygill.com/index.php?paste_id=238