## Pythagorean Triples

### October 26, 2012

Today’s exercise feels like a Project Euler problem:

A pythagorean triple consists of three positive integers a, b and c with a < b < c such that a2 + b2 = c2. For example, the three numbers 3, 4 and 5 form a pythagorean triple because 32 + 42 = 9 + 16 = 25 = 52.

The perimeter of a pythagorean triple is the sum of the three numbers that make up the pythagorean triple. For example, the perimeter of the 3, 4, 5 pythagorean triple is 3 + 4 + 5 = 12. There are 17 pythagorean triples with perimeter not exceeding 100. Ordered by ascending perimeter, they are: 3, 4, 5; 6, 8, 10; 5, 12, 13; 9, 12, 15; 8, 15, 17; 12, 16, 20; 7, 24, 25; 15, 20, 25; 10, 24, 26; 20, 21, 29; 18, 24, 30; 16, 30, 34; 21, 28, 35; 12, 35, 37; 15, 36, 39; 9, 40, 41; and 24, 32, 40.

How many pythagorean triples exist with perimeter not exceeding one million?

Your task is to write a program to compute the count of pythagorean triples with perimeter not exceeding one million. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.

Pages: 1 2

### 16 Responses to “Pythagorean Triples”

1. Andras said

Typo: 24, 30, 40 is wrong. Correct: 24, 32, 40.

2. programmingpraxis said

Fixed. Thanks.

3. […] Programming Praxis mentioned that the newest challenge sounded like a Project Euler problem, they were’t wrong. Basically, the idea is to count the […]

4. JP said

My Racket version here: Pythagorean Triples

It took a bit longer than I expected as the brute force solution would take entirely too long to run and I couldn’t quite hammer out all of the bugs using Euclid’s formula. Finally I stumbled upon Hall’s method (also used in one of the solutions here) and got a solution working.

It turns out though that generators in Racket are somewhat more expensive than I’d thought (based on a comment from programmingpraxis). So there’s now also a version that doesn’t use them that runs about 200x faster. Yay for premature optimization actually slowing things down!

5. Paul said

A Python version. It runs in about 0.15 seconds.

def pyth(N):
n = 0
Q = [(3,4,5)]
pop = Q.pop
while Q:
a, b, c = pop()
sumabc = a + b + c
if sumabc <= N:
a2, b2, c2, c3 = 2 * a, 2 * b, 2 * c, 3 * c
p1 = (-a+b2+c2, -a2+b+c2, -a2+b2+c3)
p2 = ( a+b2+c2,  a2+b+c2,  a2+b2+c3)
p3 = ( a-b2+c2,  a2-b+c2,  a2-b2+c3)
Q += [p1, p2, p3]
n += N // sumabc
return n
N = 1000000
print "{0:8d} {1:10d}".format(N, pyth(N))

6. cosmin said

Another Python solution, but not that faster.

7. cosmin said

I didn’t paste the right thing. Sorry for that :)

8. cosmin said
from math import floor, sqrt
from fractions import gcd

def count_pythagorean_triples(MAXP):
count = 0
for m in xrange(1, int(floor(sqrt(MAXP/2))) + 1):
for n in xrange (1, m):
if gcd(m, n) != 1 or gcd(gcd(m*m - n*n, 2*m*n), m*m + n*n) != 1: continue
count += int(floor(MAXP/(2*m*(m + n))))
return count

print count_pythagorean_triples(1000000)

9. Paul said

Cosmin, the second condition on line 8 is not necessary. Also the loop in line 7 has only to loop over the even ints if m is odd and vice versa (you can change to:
for n in xrange(m&1 + 1, m, 2)). Then your code runs a lot faster (on my PC 0.17 sec) and produces the same result.

10. cosmin said

Paul, you are right about the step=2 trick. The algorithm would run 2x faster. Initially I thought that the second condition is redundant and completely unnecessary, but I found a simple counterexample. m=3, n=1 generates (8, 6, 10) which is not a primitive triple, because it can be obtained from (3, 4, 5) (m=2 and n=1) and in this way some of the solutions will be counted multiple times.

11. Paul said

Cosmin. If you make sure, that only odd-even coprimes are used, then you have no problem. With odd-even coprimes only primitive triples are formed. Therefore it is important to loop only over even n if m is odd and vice versa. If m,n are odd-odd coprimes you will be double counting and you need the second condition on line 8.

12. Graham said

In Shen (without type checking, since I couldn’t quite get it to work):

\* Counts all triples with perimeter less than n, given a method of generating
primitive Pythagorean triples
*\

(define neg
N -> (- 0 N))

(define hall
N -> (hall-1 N [3 4 5]) where (and (natural? N) (positive? N))
_ -> [[]])

(define hall-1
N [A B C] -> (if (< N (+ A B C))
[]
(append
[(if (< A B) [A B C] [B A C])]
(hall-1 N [(+ A (neg B) (neg B) C C)
(+ A A (neg B) C C)
(+ A A (neg B) (neg B) C C C)])
(hall-1 N [(+ A B B C C)
(+ A A B C C)
(+ A A B B C C C)])
(hall-1 N [(+ (neg A) B B C C)
(+ (neg A) (neg A) B C C)
(+ (neg A) (neg A) B B C C C)]))))

(define f
Pyth N -> (sum (map (/. P (div N P)) (map sum (Pyth N)))))

Basically a translation of the given solution as I try to learn Shen.

13. Graham said

Apologies for botching posting the code.

14. ardnew said

another typo: “21, 27, 35” should be “21, 28, 35”

15. programmingpraxis said

Fixed. Thanks again.

16. ardnew said

This is basically the same approach that Paul took. Performs the count for P < 1000000 in just under a second on my machine.

First command line argument defines P. If a second argument exists, all triples are generated and printed to STDOUT (significantly slower).

use strict;
use warnings;
use List::Util qw( sum );

# matrices discovered by F.J.M. Barning
# (source: http://oai.cwi.nl/oai/asset/7151/7151A.pdf)
our \$A = [[ 1, -2, 2], [ 2, -1, 2], [ 2, -2, 3]];
our \$B = [[ 1, 2, 2], [ 2, 1, 2], [ 2, 2, 3]];
our \$C = [[ -1, 2, 2], [ -2, 1, 2], [ -2, 2, 3]];

# returns three primitive triples from a given input primitive triple
sub mprod
{
my (\$a, \$b, \$c) = @_;

return map
{[
\$a * \$\$_[0][0] + \$b * \$\$_[0][1] + \$c * \$\$_[0][2],
\$a * \$\$_[1][0] + \$b * \$\$_[1][1] + \$c * \$\$_[1][2],
\$a * \$\$_[2][0] + \$b * \$\$_[2][1] + \$c * \$\$_[2][2]
]}
\$A, \$B, \$C;
}

# prints the multiples of a primitive triple, including itself, until
# their elements' sum is greater than n
sub pprim
{
my (\$a, \$b, \$c, \$n) = @_;
my (\$x, \$y, \$z, \$k) = (\$a, \$b, \$c, 1);

while (\$x + \$y + \$z <= \$n)
{
printf "%8d %8d %8d\$/", \$x, \$y, \$z;
++\$k;
\$x = \$a * \$k;
\$y = \$b * \$k;
\$z = \$c * \$k;
}

return \$k – 1;
}

# counts all Pythagorean triples less than N. if parameter \$m is
# non-zero, then all primitive and multiples are printed to STDOUT
sub generate
{
my \$n = shift;
my \$m = shift;
my @r = [@_];

my (\$s, \$c) = (0, 0);

while (@r)
{
my \$v = pop @r;
my \$p = sum @\$v;
if (\$p < \$n)
{
\$c += pprim(@\$v, \$n) if \$m;
\$s += int(\$n / \$p);
push @r, mprod(@\$v);
}
}

printf "\$/generated triples: %10d".
"\$/sum of quotients: %10d\$/\$/", \$c, \$s;
}

die "usage:\$/\tperl \$0 <upper limit> [0|1 (print triples)]\$/\$/"
unless @ARGV > 0;

generate(\$ARGV[0], @ARGV > 1, (3, 4, 5));