## Pythagorean Triples

### October 26, 2012

I have long been fascinated by how much the ancient mathematicians were able to accomplish with their limited resources, and today’s exercise offers an opportunity to look at one of those accomplishments.

Our solution today uses an algorithm that was codified by Euclid (Elements, Book X, Proposition 29) about 300 B.C. but was probably know to Pythagoras two hundred years before that. Take any two numbers m and n, with n > n, with m and n relatively prime and with one of them even and the other odd. Then (m2n2, 2 m n, m2 + n2) is a primitive pythagorean triple (a, b, c) in which a and b are relatively prime. Since every m, n pair generates a different triple, you can generate all the primitive triples by iterating through all combinations of m and n:

```(define (euclid limit)   (define (trip m n)     (let ((m2 (* m m)) (n2 (* n n)))       (let ((a (- m2 n2)) (b (* 2 m n)) (c (+ m2 n2)))         (if (< a b) (list a b c) (list b a c)))))   (let m-loop ((m 1) (trips (list)))     (if (< limit (* m m)) trips       (let n-loop ((n (+ (modulo m 2) 1)) (trips trips))         (if (<= m n) (m-loop (+ m 1) trips)           (let* ((t (trip m n)) (p (apply + t)))             (if (< limit p)                 (m-loop (+ m 1) trips)                 (n-loop (+ n 2) (if (= (gcd (car t) (cadr t)) 1)                                     (cons t trips) trips)))))))))```

There are seven primitive pythagorean triples with perimeter less than a hundred:

```> (euclid 100) ((12 35 37) (9 40 41) (20 21 29) (7 24 25) (8 15 17) (5 12 13) (3 4 5))```

A more modern method comes to us from a British math professor named Hall. He defines a genealogy of primitive pythagorean triples: if a, b, c is a primitive pythagorean triple, so too are:

a – 2b + 2c, 2ab + 2c, 2a – 2b + 3c

a + 2b + 2c, 2a + b + 2c, 2a + 2b + 3c

-a + 2b + 2c, =2a + b + 2c, -2a + 2b + 3c

The code is straight forward, and quite pretty:

```(define (hall n)   (let loop ((a 3) (b 4) (c 5))     (if (< n (+ a b c)) (list)       (append         (list (if (< a b) (list a b c) (list b a c)))         (loop (+ a (- b) (- b) c c)               (+ a a (- b) c c)               (+ a a (- b) (- b) c c c))         (loop (+ a b b c c)               (+ a a b c c)               (+ a a b b c c c))         (loop (+ (- a) b b c c)               (+ (- a) (- a) b c c)               (+ (- a) (- a) b b c c c))))))```

Once again, there are seven pythagorean triples with perimeter less than a hundred:

```> (hall 100) ((3 4 5) (5 12 13) (7 24 25) (9 40 41) (20 21 29) (8 15 17) (12 35 37))```

Whatever method we use for generating primitive pythagorean triples, it is easy to count all the triples with perimeter less than n; simply sum for all the triples the quotient of n and the triple’s perimeter:

```(define (f pyth n)   (apply +     (map (lambda (p) (quotient n p))          (map (lambda (xs) (apply + xs))               (pyth n)))))```

```> (f euclid 1000000) 808950```
```> (f hall 1000000) 808950```

You can run the program at http://programmingpraxis.codepad.org/pnV0GZKI.

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 * \$\$_ + \$b * \$\$_ + \$c * \$\$_,
\$a * \$\$_ + \$b * \$\$_ + \$c * \$\$_,
\$a * \$\$_ + \$b * \$\$_ + \$c * \$\$_
]}
\$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, @ARGV > 1, (3, 4, 5));