23
\$\begingroup\$

I wrote this code for Project Euler problem 75, which asks how many integers ≤ 1500000 exist, where the integer is a perimeter length that can be divided into three integer-length sides of a right triangle in one unique way.

I was curious if anyone knew how to improve its speed. It runs fine, but I'm just looking to improve my own coding know-how.

from functools import reduce
import math

primes=set([2,3,5,7,11,13,17,19,23])
def isPrime(n):
    n=abs(n)
    if n in primes:
       return True
    if n<2:
       return False
    if n==2:
       return True
    if n%2==0:
       return False
    for x in range(3, int(n**0.5)+1, 2):
        if n % x == 0:
            return False
    primes.add(n)
    return True

def aFactors(n):
    if isPrime(n):
        return [1,n]
    return set(reduce(list.__add__,([i,n//i] for i in range(1,int(math.sqrt(n))+1) if n%i==0)))


count=0
number=12
while number<=1500000:
    p=number/2
    f=aFactors(p)
    triangles=[]
    pairs=[(i, int(p/i)) for i in f]
    add=triangles.append
    for i in pairs:
        mList=aFactors(i[0])
        pairsOfmc=[(k,int(i[0]/k)) for k in mList]
        for j in pairsOfmc:
            add((2*i[1]*i[0]-i[1]*i[1]*j[0],2*i[0]*i[1]-2*i[0]*j[1],2*i[0]*j[1]+i[1]*i[1]*j[0]-2*i[0]*i[1]))
    r=0
    while r<len(triangles):
        if any(triangles[r][i]<=0 for i in range(len(triangles[r]))):
            del triangles[r]
        else:
            l=list(triangles[r])
            l.sort()
            triangles[r]=tuple(l)
            r+=1
    trianglesFinal=list(set(triangles))
    for i in trianglesFinal:
        print(number, i)
    if len(trianglesFinal)==1:
        count+=1
    number+=2
print(count)

Please note that I am not looking for a different calculating method (I am sure there is one, but, for me, Project Euler is about finding your own methods. Using yours would, to me, be cheating). However, any faster functions, ways to combine blocks of code, simplified tests, or the like (such as not checking ever other number, but every xth number, etc) would be greatly appreciated!

\$\endgroup\$
  • \$\begingroup\$ Probably not what you are looking for but in my opinion if you are worried about run time, then python is not the right language to use. It is hard to fine tune and optimize in general. It is design for rapid prototyping more then for efficient code. \$\endgroup\$ – gbtimmon Apr 23 '13 at 20:06
49
+50
\$\begingroup\$

For the record, this is Project Euler problem 75.

1. Admonition

You write:

Please note that I am not looking for a different calculating method (I am sure there is one, but, for me, Project Euler is about finding your own methods. Using yours would, to me, be cheating.)

Abandon this attitude! A key part of programming — just as important as coding — is developing a wide repertoire of algorithms and techniques that you can apply to problems you face. Although for your own intellectual satisfaction it's great to discover algorithms on your own, you should always go on to compare your solution to the best solutions discovered by others, so that you can improve your repertoire for the next problem.

In particular, when you want to speed up a program, it's counterproductive to say, "I don't want to implement a different algorithm, I just want to speed up the code I wrote". The biggest speedups come from finding better algorithms.

You say that your program "runs fine" but that doesn't seem to be true. Project Euler says:

Each problem has been designed according to a "one-minute rule", which means that although it may take several hours to design a successful algorithm with more difficult problems, an efficient implementation will allow a solution to be obtained on a modestly powered computer in less than one minute.

I tried running your program, but after ten minutes it still had not produced an answer (and my laptop was running hot), so I killed it.

So in section 2 below I'll be using a smaller test size, 100,000 instead of 1,500,000, to keep the runtimes manageable.

2. Piecewise optimization

  1. Your program is hard to test from the interactive interpreter because it has code at the top level. Better to put the main program inside a function so that you can call it from interpreter and then add an if __name__ == '__main__': section so that you can run it as a script if you want to.

    I've used the function name problem75, and it takes an argument, which is the largest value of the length of wire in the problem. Now:

    >>> from timeit import timeit
    >>> timeit(lambda:problem75(10**5),number=1)
    12 (3, 4, 5)
    24 (6, 8, 10)
    30 (5, 12, 13)
    [... much output deleted ...]
    33.87288845295552
    
  2. The print statement near the end of the main loop is unnecessary. This saves about 3 seconds, bringing the time down to 31.0 seconds.

  3. Your function aFactors makes an initial call to isPrime, taking \$ O(\sqrt n) \$, in order to avoid a loop that also takes \$ O(\sqrt n) \$. The test costs as much as it saves, so it's not worth it. Remove the call to isPrime and rewrite the function like like this:

    def factors(n):
        """Return the set of factors of n."""
        result = set([1, n])
        add = result.add
        for i in range(2, int(math.sqrt(n)) + 1):
            if n % i == 0:
                add(i)
                add(n // i)
        return result
    

    Here I've picked a better name for the function, written a docstring explaining what it does, and cached the value of result.add so that it does not need to be looked up each time arounnd the loop. This saves about 4 seconds; time now 27.7 seconds.

  4. Since the only thing you do with the set of factors is to turn it into pairs which you then iterate over, why not generate the pairs, like this:

    def product_pairs(n):
        """Generate pairs (i, j) such that i * j = n."""
        yield 1, n
        yield n, 1
        for i in range(2, int(math.sqrt(n)) + 1):
            if n % i == 0:
                j = n // i
                yield i, j
                yield j, i
    

    and then process them like this:

    triangles = []
    append = triangles.append
    for i in product_pairs(number // 2):
        for j in product_pairs(i[0]):
            append((2*i[1]*i[0]-i[1]*i[1]*j[0],2*i[0]*i[1]-2*i[0]*j[1],2*i[0]*j[1]+i[1]*i[1]*j[0]-2*i[0]*i[1]))
    

    This avoids having to construct an intermediate list in memory. This saves about 3 seconds; time now 24.3 seconds.

  5. Instead of writing for i in product_pairs(...): and then looking up i[0] and i[1], split up the pair into two variables when you assign it in the loop, like this:

    for i, j in product_pairs(number // 2):
        for k, l in product_pairs(i):
            append((2*j*i - j*j*k, 2*i*j - 2*i*l, 2*i*l + j*j*k - 2*i*j))
    

    This avoids the sequence lookups. This saves about 3 seconds; time now 21.2 seconds.

  6. You construct a list of triangles, some of which have sides with negative or zero length. You then go through the list and del the triangles which are invalid. But del on a list is a potentially expensive operation: it has to copy the remainder of the list to keep it contiguous. See the TimeComplexity page on the Python wiki for the cost of basic operations on Python's built-in data structures, where you can see that "delete item" operation on a list takes \$O(n)\$.

    So don't add these invalid triangles to the list in the first place. In fact, don't bother with the list at all, just construct the set of triangles directly:

    triangles = set()
    add = triangles.add
    for i, j in product_pairs(number // 2):
        for k, l in product_pairs(i):
            triangle = 2*j*i - j*j*k, 2*i*j - 2*i*l, 2*i*l + j*j*k - 2*i*j
            if all(side > 0 for side in triangle):
                add(tuple(sorted(triangle)))
    

    This saves about 6 seconds; time now 15.1 seconds.

  7. You can avoid some repeated multiplications by naming the various products of i, j, k, and l:

    for i, j in product_pairs(number // 2):
        for k, l in product_pairs(i):
            p, q, r = 2 * j * i, j * j * k, 2 * i * l
            triangle = p - q, p - r, r + q - p
            if all(side > 0 for side in triangle):
                add(tuple(sorted(triangle)))
    

    This saves about 1 second; time now 13.9 seconds.

And I'm afraid that's as far as I got by applying piecewise speedups to your code. About 60% speedup in total. The revised program now gets the solution to the full Project Euler problem in about eight minutes on my laptop:

>>> timeit(lambda:main(1500000),number=1)
[... solution deleted ...]
467.3059961795807

Perhaps some other participant here can make more piecewise progress, but for dramatic speedups, you really need ...

3. A better algorithm

(Avert your eyes if you must.)

The basic insight is that instead of iterating over the perimeter and finding the Pythagorean triples with that perimeter, we can iterate over all Pythagorean triples in some more convenient order, and keep count of how many triangles are found at each perimeter length.

We can further speed things up by only generating the primitive Pythagorean triples, and then multiplying the primitive triples by successive numbers 1, 2, 3, ... to generate the remaining Pythagorean triples in the required range. In fact, we only need to multiply their perimeters.

Euclid’s formula can be used to generate all primitive Pythagorean triples. Given coprime positive integers \$m\$ and \$n\$, with \$m > n\$, and exactly one of \$m\$ and \$n\$ even, $$ \eqalign{ a &= m^2 − n^2 \\ b &= 2mn \\ c &= m^2 + n^2} $$ is a primitive Pythagorean triple.

The remaining subtlety is that Euclid’s formula doesn’t generate triples in order by length of the perimeter, so there needs to be some way to determine when to stop. My strategy in the code below is to iterate over ascending values of \$m\$. The perimeter of the triple is \$a + b + c = 2m^2 + 2mn\$, which is at least \$2m(m + 1)\$, since \$n ≥ 1\$. So \$2m(m + 1)\$ is a lower bound on the perimeter for which all triples have been generated so far. When this exceeds the limit, we can stop the search: there are no more triples to be found in the required range.

from collections import Counter
from fractions import gcd
from itertools import count

def coprime(m, n):
    """Return True iff m and n are coprime."""
    return gcd(m, n) == 1

def all_primitive_triples():
    """Generate all primitive Pythagorean triples, together with a lower
    bound on the perimeter for which all triples have been generated
    so far.

    """
    for m in count(1):
        for n in range(1, m):
            if (m + n) % 2 and coprime(m, n):
                a = m * m - n * n
                b = 2 * m * n
                yield a, b, m * m + n * n, 2 * m * (m + 1)

def problem75(limit=1500000):
    """Return the number of values of L <= limit such that there is
    exactly one integer-sided right-angled triangle with perimeter
    L.

    """
    triangles = Counter()
    for a, b, c, q in all_primitive_triples():
        if q > limit:
            break
        p = a + b + c
        for i in range(p, limit + 1, p):
            triangles[i] += 1
    return sum(n == 1 for n in triangles.values())

This solves the Project Euler problem in a more reasonable amount of time:

>>> timeit(problem75, number=1)
2.164217948913574
\$\endgroup\$
  • 3
    \$\begingroup\$ On the subject of more piecewise progress, the double for-loop to split a number into 3 factors can be replaced with a method which returns product triples. I got a further 30% speedup (from 9 secs to 6 secs) this way. \$\endgroup\$ – Peter Taylor Apr 23 '13 at 22:58
  • \$\begingroup\$ Pure awesomeness, codereview really needs a +2 \$\endgroup\$ – konijn Apr 23 '13 at 23:27
  • 4
    \$\begingroup\$ Half tempted to edit this post and add all manner of bolding to The biggest speedups come from finding better algorithms. \$\endgroup\$ – Ben Jackson Apr 24 '13 at 0:33
  • \$\begingroup\$ I found the linear algebra approach to generateing all primative pythagorean tripples more intuitive. \$\endgroup\$ – recursion.ninja Feb 18 '14 at 14:31
  • \$\begingroup\$ @awashburn: With the linear algebra approach, it's not so easy to see how to bound the perimeter for which all triples have been generated so far, which is a crucial part of solving this problem. \$\endgroup\$ – Gareth Rees Feb 18 '14 at 14:35
3
\$\begingroup\$

It is not necessary to construct all the triangles to find out if there is exactly one: you can stop when you've found two.

Some little things to clean up the code that should give a minor speedup too:

Instead of indexing (triangles[r][i]<=0 for i in range(len(triangles[r]))), iterate directly:

(x <= 0 for x in triangles[r])

reduce(list.__add__, ... looks ugly to me and creates some intermediate lists you can avoid by:

(j for i in range(1,int(math.sqrt(n))+1) if n%i==0 for j in (i,n//i))

Instead of the outer while loop you can use:

for number in range(12, 1500001, 2): 

(This is for Python 3, in Python 2 xrange should be used)

\$\endgroup\$
  • \$\begingroup\$ I tried out the "stop when you've found two triangles" approach, but found that it actually made things slower. (The test doesn't succeed often enough to pay for its own cost.) \$\endgroup\$ – Gareth Rees Apr 24 '13 at 8:36
  • \$\begingroup\$ @GarethRees Is this in reference to your best effort? I get almost 40 % time reduction using this trick, when my starting point is the OP's code plus an optimization similar to your number 6. \$\endgroup\$ – Janne Karila Apr 24 '13 at 9:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.