## 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 (*m*^{2} − *n*^{2}, 2 *m* *n*, *m*^{2} + *n*^{2}) 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, 2a–b+ 2c, 2a– 2b+ 3ca + 2b + 2c, 2a + b + 2c, 2a + 2b + 3

c-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))

There are many other algorithms for generating pythagorean triples; search Google for “generating pythagorean triples” to learn more about this fascinating topic.

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

Now we are ready to answer the challenge question:

`> (f euclid 1000000)`

808950

`> (f hall 1000000)`

808950

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

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

Fixed. Thanks.

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

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!

A Python version. It runs in about 0.15 seconds.

Another Python solution, but not that faster.

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

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.

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.

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.

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

*\

(load "shen-libs/maths-lib.shen")

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

Apologies for botching posting the code.

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

Fixed. Thanks again.

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));