## Pythagorean Triples

### October 26, 2012

Today’s exercise feels like a Project Euler problem:

A pythagorean triple consists of three positive integers

a,bandcwitha<b<csuch thata^{2}+b^{2}=c^{2}. For example, the three numbers 3, 4 and 5 form a pythagorean triple because 3^{2}+ 4^{2}= 9 + 16 = 25 = 5^{2}.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.

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