Project Euler Problem 9 is a tough one. After reading the Pythagorean Triple Wikipedia page, I implemented a NumPy solution with Euclid’s Formula.

1. Create m and n arrays
The Euclid’s Formula defines indices m and n. We will create arrays to hold these indices.
1 2 | m = numpy.arange(33) n = numpy.arange(33) |
2. Calculate a, b and c
The second step is to calculate a, b and c of the Pythagorean triple using Euclid’s formula. Use the outer function to get the cartesian products, difference and sums we require.
1 2 3 | a = numpy.subtract.outer(m ** 2, n ** 2) b = 2 * numpy.multiply.outer(m, n) c = numpy.add.outer(m ** 2, n ** 2) |
3. Find the index
At this point we have a number of arrays containing a, b and c values. However, we still need to find the values that conform to the problem’s condition. Find the index of those values with the NumPy where function.
1 | idx = numpy.where((a + b + c) == 1000) |
4. Check solution
Check the solution with the numpy.testing module.
1 | numpy.testing.assert_equal(a[idx]**2 + b[idx]**2, c[idx]**2) |
Below is the complete code or get it from Github
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 | import numpy import numpy.testing #A Pythagorean triplet is a set of three natural numbers, a < b < c, for which, #a ** 2 + b ** 2 = c ** 2 # #For example, 32 + 42 = 9 + 16 = 25 = 52. # #There exists exactly one Pythagorean triplet for which a + b + c = 1000. #Find the product abc. #1. Create m and n arrays m = numpy.arange(33) n = numpy.arange(33) #2. Calculate a, b and c a = numpy.subtract.outer(m ** 2, n ** 2) b = 2 * numpy.multiply.outer(m, n) c = numpy.add.outer(m ** 2, n ** 2) #3. Find the index idx = numpy.where((a + b + c) == 1000) #4. Check solution numpy.testing.assert_equal(a[idx]**2 + b[idx]**2, c[idx]**2) print a[idx] * b[idx] * c[idx] |
If you liked this post and are interested in NumPy check out NumPy Beginner’s Guide by yours truly.



