How to optimize this Python script?

Question

There are numerous solutions (a, b, c) for an integer right triangle with perimeter p, and for all these solutions a+b+c == p and the Pythagorean theorem also applies. I am writing a Python script to calculate the maximum number of solutions possible for a triangle with perimeter <= 1000.

My script is correct, however it takes forever to run. I'm sure it would take more than 30 minutes even with my i7 processor, so I need to optimize it. Can someone help me? (This is a problem on Project Euler, in case anyone is wondering)

def solutions(p):
    result = []

    for a in range(1, p + 1):
        for b in range(1, p - a + 1):
            for c in range(1, p - a - b + 1):
                if a + b + c == p and a < b and b < c:
                    d = a ** 2
                    e = b ** 2
                    f = c ** 2

                    if (d + e == f) or (e + f == d) or (f + d == e):
                        result.append((a, b, c))
    return len(result)


max_p = 0
max_solutions = 0

for p in range(3, 1001):
    print("Processing %d" % p)
    s = solutions(p)

    if s > max_solutions:
        max_solutions = s
        max_p = p

print("%d has %d solutions" % (max_p, max_solutions))

Answer 1

This isn't your code optimized, but my own code (which I used for this problem). I started off by doing some algebra to make the program very easy and not have to iterate 1000^3 times (1-1000 for a, then 1-1000 for b for every value of a, and 1-1000 for c for every value of b.

# Project Euler 9
'''
Algebra behind Final Method:
a + b + c = 1000 | a^2 + b^2 = c^2
c = 1000 - (a + b) # Solving for C
a^2 + b^2 = (1000 - (a + b))^2 # Substituting value for C
a^2 + b^2 = 1000000 - 2000a - 2000b + (a + b)^2 # simplifying
a^2 + b^2 = 1000000 - 2000a - 2000b + a^2 + b^2 + 2ab # simplifying again
0 = 1000000 - 2000a - 2000b + 2ab # new formula
2000a - 2ab = 1000000 - 2000b # isolating A
1000a - ab = 500000 - 1000b # divide by 2 to simplify
a(1000 - b) = 500000 - 1000b # factor out A
a = (500000 - 1000b) / (1000 - b) # solve for A
'''
def pE9():
   from math import sqrt
   a, b, c = 1, 1, 1
   while True:
      b += 1
      a = (500000 - 1000 * b) / (1000 - b)
      c = sqrt(a**2 + b**2)
      if a + b + c == 1000 and a**2 + b**2 == c**2:
         break
   print int(a * b * c)

from timeit import timeit
print timeit(pE9, number = 1)

Using number = 1 so just test it for one run.

Output:

>>> 
# Answer hidden
0.0142664994414

Answer 2

Here is my rewrite of your program.

First, I precompute all the squared values. Not only does that avoid the multiplication, but it means that Python won't be constantly creating and garbage-collecting objects for all the square values.

Next, I got rid of the loop for the third side of the triangle. Once you have chosen values for a and b there is only one possible value that fits the criteria a + b + c == 1000 so this just tests that one. This turns the problem from approximately O(n^3) to approximately O(n^2), a vast improvement.

Then I tried running it. On my four-year-old computer it finished in about 46 seconds, so I stopped there and here you go.

I did a Google search and found discussion of this problem; if the discussion I saw was correct, then this program prints the correct answer.

upper_bound = 1000

sqr = [i**2 for i in range(upper_bound+1)]

def solutions(p):
    result = []

    for a in range(1, p - 1):
        for b in range(1, p - a):
            c = p - (a + b)
            if a < b < c:
                d = sqr[a]
                e = sqr[b]
                f = sqr[c]

                if (d + e == f) or (e + f == d) or (f + d == e):
                    result.append((a, b, c))
    return len(result)


max_p = 0
max_solutions = 0

for p in range(3, upper_bound+1):
    print("Processing %d" % p)
    s = solutions(p)

    if s > max_solutions:
        max_solutions = s
        max_p = p

print("%d has %d solutions" % (max_p, max_solutions))

EDIT: Here's a somewhat faster version I was playing around with. It incorporates a suggestion from @gnibbler in the comments.

upper_bound = 1000

sqr = [i**2 for i in range(upper_bound+1)]

def solution(p):
    count = 0
    for a in range(1, p - 1):
        for b in range(a, p - a):
            c = p - (a + b)
            d = sqr[a]
            e = sqr[b]
            f = sqr[c]

            if (d + e == f):
                count += 1
    return count

c, p = max((solution(p), p) for p in range(3, upper_bound+1))
print("%d has %d solutions" % (p, c))

On my computer this version takes 31 seconds instead of 46.

The tricky business with using max() doesn't really make it noticeably faster. I tried it without pre-computing the squares and it was very much slower, so slow I didn't want to wait for an exact time.

Answer 3

Got it. It just relied on setting a^2+b^2=c^2 and then substituting p - a - b = c

 1 from math import pow
  2 
  3 def see_if_right_triangle(p):
        solutions = 0
  4     # Accepts the perimeter as input
  5     for a in range(1, p): 
  6         for b in range(1, p):
  7             if 2*p*b + 2*p*a - pow(p, 2) == 2*a*b:
  8                solutions += 1
        print "The perimeter {p} has {sol} number of solutions".format(p=p, sol=solutions)
 10                        
 11 
 12 for p in range(3, 1001):
 13     see_if_right_triangle(p)

I think this can be optimized more... especially if you figure out some math to narrow down the ranges you are going to accept for a and b

Answer 4

A better one:

def solution(n):
    count = 0
    for c in range(n // 3 + 1, n // 2):
        for a in range(1, n // 3):
            b = n - a - c
            if b <= 0:
                continue
            if a >= b:
                continue
            if a * a + b * b != c * c:
                continue
            count += 1
    return count