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.

About these ads

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

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

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Follow

Get every new post delivered to your Inbox.

Join 635 other followers

%d bloggers like this: