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.
Programming Praxis 
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.
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))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));