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
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
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.
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.
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.
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.
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.
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
