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

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.

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 574 other followers

%d bloggers like this: