A 1000 Pythagorean triplets – Problem 9

Written by on 15 November 2010

Topics: Project Euler

Today’s problem in Project Euler, is related to probably THE most well known theorem of all times.  Pythagoras theorem stating that Geometric interpretation of Pythagoras

The square of the hypotenuse is the sum of squares of the two other sides.

It can also be stated in a more geometrical way as

In any right triangle, the area of the square whose side is the hypotenuse (the side opposite the right angle) is equal to the sum of the areas of the squares
whose sides are the two legs (the two sides that meet at a right angle

Which can also be shown in a graphical sense, on the figure to the right, where the blue area equals the orange area.

But enough with the background info, the problem reads

A Pythagorean triplet is a set of three natural numbers, a < b < c, for which,

a2 + b2 = c2

For example, 32 + 42 = 9 + 16 = 25 = 52.

There exists exactly one Pythagorean triplet for which a + b + c = 1000.

Find the product abc.

There is an infinite number of Pythagorean triplets, so I am quite glad for the last condition stating that the circumference of the triangle has to be 1000. As a side note Fermat conjectured, and wrote that he had a proof in the margin of his text book that an + bn = cn has no integer solutions for n > 2. It has only recently been proven that he was right.

A great book on that topic is Fermat’s Last Theorem By Simon Singh. But that is just a side track.

There are two approaches for finding a solution in my post, there is the brute force method, which is widely described on the net, and works brilliantly for small numbers, and then there is a number theoretical approach which I myself had to study. So providing the latter is a test for me to see if I understood it to a level where I can explain it to you. And don’t hesitate to ask for clarification or point out errors in my deductions.

In my further deduction I will denote the circumference of the triangle s, which in the problem is fixed to 1000.

Brute Force

For the brute force algorithm we would need to loop over all the possible values of a and b, calculate c and check if that is a solution. We could be ignorant and just let a and b loop from 1-1000, or we could use the facts that a < b < c, and thus exploit that  a < s/3,  and a < b < s/2.

The code can be written as

int a = 0, b = 0, c = 0;
int s = 1000;
bool found = false;
for (a = 1; a < s / 3; a++) {
    for (b = a; b < s / 2; b++) {
        c = s - a - b;

        if (a * a + b * b == c * c) {
            found = true;
            break;
        }
    }

    if (found) {
        break;
    }
}

and the result

The Pythagorean triple is 200, 375, 425, and the sum is 1000
The product is 31875000
Solution took 0 ms

Number Theoretical Approach

As I mentioned in the beginning of the post, this approach took me a while to figure out, so I have made my explanation rather elaborate to make it understandable for myself.  The rest of this post will be more mathematical rigorous than any of my previous posts. I hope you enjoy reading it as much as I did writing it.

Prerequisites and Definitions

More than 2000 years ago Euclid established that Pythagorean Triplets can be generated by the formula.

Given any positive integers m and n where m > n > 0, the integers

a =m2-n2

b =2m*n

c =m2+n2

forms a Pythagorean triplet.  You can convince your self that (m2-n2)2 + (2m*n)2 = (m2+n2)2 by expanding both sides.  There are quite a number of other formulas for generating triplets, but this is the one I will use here. Depending on the choices of m and n, we might have to interchange a and b to adhere to the constraint a < b.

Furthermore we need a definition of primitive.

Definition: Primitive

A Pythagorean Triplet is called primitive if a, b and c are coprime. That means that the greatest common divisor (gcd) of the set {a,b,c} is one.

A Pythagorean triplet is primitive if and only if  m and n are coprimes and exactly one of them is even.

Proof

Showing that it is a necessary condition is rather easy. Assume that gcd(n,m)= d > 1, then d2 is a common divisor for {a,b,c} (you can convince your self about it by inserting m = d*p and n = d*q), and if both n and m are odd then a,b and c are even, and thus 2 is a common divisor. Thus these conditions are necessary.

Showing that they are sufficient is a bit more tricky.

Assume that exactly one of  m and n is even and and triplet is not primitive. Since the triplet is not primitive a prime p is divisible with a, b and c (according to the fundamental theorem of arithmetic as referred to in  problem 3 such a prime exists).

From the condition that exactly one of m and n is odd, it follows that a and c is also odd, and so is p then. Any odd prime that that divides b must also divide m and/or n, without loss of generality assume that it divides m. Thus it must also divide m2. Since the prime is also a divisor of c = m2 + n2, then it follows that it must divide n2 and n as well. (this follows from notions and facts part of the divisor article on wikipedia) Therefore p is a divisor of m and n, and therefore m and n are not coprimes. Thus it is sufficient that the conditions are fulfilled.

Finding a new unique representation

Every primitive Pythagorean triplet can be represented by a unique pair of coprimes m and n where exactly one of them is even. However, not all Pythagorean triplets are generated this way. Approaching it from the other end. From every Pythagorean triplet a primitive triplet can be generated by dividing by the greatest common divisor, thus every Pythagorean triplet can be represented by the unique set {d,m,n} by the given formula

a =d(m2-n2)

b =2d*m*n

c =d(m2+n2)

with m > n > 0, m and n being coprimes and exactly one of m and n even. d is the greatest common divisor of a, b and c.

Solving the problem

Now we have a new unique representation of any triplet, but how does that help us solving the problem?

Using the conditions

a + b+ c = s = d(m2-n2) + 2d*m*n + d(m2+n2) = 2m(m+n)d

We can see that m must be a divisor to s/2, and

and we need to find a k = m+n such that k is a divisor to s/(2m) and

m <k < 2m,

k < s/(2m),

k is odd and gcd(m,k) = 1.

In our problem it means that m < 23, thus decreasing the first number we need to find by a factor 10. The second number k, we need to find has an upper bound depending on m but is in the same vicinity. Once we have found m and k, the rest can be calculated easily. So if we can find a cheap way to calculate the greatest common divisor, we can make an efficient algorithm for finding the triplet.

Algorithm for greatest common divisor

On wikipedia there is a few different methods for calculating the greatest common divisor. I have chosen Euclid’s  algorithm, which is fairly fast and easy to implement. The code can be written as

public int gcd(int a, int b) {
    int y = 0;
    int x = 0;

    if (a > b) {
        x = a;
        y = b;
    } else {
        x = b;
        y = a;
    }

    while (x % y != 0) {
        int temp = x;
        x = y;
        y = temp % x;
    }
    return y;
}

Finding the Pythagorean triplet

We now have all parts we need to construct an algorithm, and a straight forward approach is

int a=0, b=0, c=0;
int s = 1000;
int m = 0, k = 0, n = 0, d = 0;
bool found = false;

int mlimit = (int) Math.Sqrt(s / 2);
for (m = 2; m <= mlimit; m++) {
    if ((s / 2) % m == 0) { // m found
        if (m % 2 == 0) { // ensure that we find an odd number for k
            k = m + 1;
        } else {
            k = m + 2;
        }
        while (k < 2 * m && k <= s / (2 * m)) {
            if (s / (2 * m) % k == 0 && gcd(k, m) == 1) {
                d = s / 2 / (k * m);
                n = k - m;
                a = d*(m * m - n * n);
                b = 2 * d * n * m;
                c = d * (m * m + n * n);
                found = true;
                break;
            }
            k += 2;
        }
    }
    if (found) {
        break;
    }
}

The result looks like

The pythagorean triple is 375, 200, 425, and the sum is 1000
The product is 31875000
Solution took 0 ms

Wrapping up

I am aware that this has been a rather long post. But I needed to explain my self, to really understand the underlying theory of the number theoretical approach. I hope you could follow it and learn something useful , otherwise ask questions and let me know if I can improve parts of the explanation.

As usual you can download the source code

18 Comments For This Post I'd Love to Hear Yours!

  1. Phil says:

    Could you explain the sufficency part of your proof for “A Pythagorean triplet is primitive if and only if m and n are coprimes and exactly one of them is even” a little more.

    If we simplify the orginal statement using the following substitutions
    A = Is a Pythagorean triplet
    B = m,n coprime
    C = exactly one of m,n is even

    the original statement breaks down to
    A B ^ C

    and in the first part of the proof you have shown
    A => B ^ C
    by showing
    !B U !C => !A
    (please correct me if I am wrong)

    I assume in the second part of the proof you are showing that
    B ^ C => A
    but I cannot follow this.

    Please can you further explain this part of the proof, if possible using the substitutions (A,B,C) that I have used.

  2. Phil says:

    my previous post should read

    A B ^ C

    where means if and only if

  3. Phil says:

    Arghhh. The characters are not showing in the post.

    Let’s just say that part should read

    A if and only if B ^ C

  4. Kristian says:

    Hi Phil

    I am not sure how much better I can explain the sufficiency part. Especially since it is a long time ago.

    What I do though is that I assume that we have B and C is fulfilled but we have found a non primtive triplet, meaning ~A. I then show that this assumption will lead to a contradiction. Since I end up showing that if we assume ~A and B then we will end up with a condition that infers ~C

    proof by contradiction is better explained in this post

    http://www.mathblog.dk/proof-by-contradiction/

    Does this help?

  5. Phil says:

    Hi Kristian,

    That did help, I enjoyed reading the post and found the Book of Proof you link to really helpful.

    I think I can now convince myself of the proof for B^C=>A:

    Consider the following truth table

    B C A B^C=>A

    T T T T
    T T F F
    T F T T
    T F F T
    F T T T
    F T F T
    F F T T
    F F F T

    Suppose that B^C=>A is false, from the above we can see that this must mean that B and C are true and A is false.
    However it can also be shown that ~A^C=>~B Therefore B is true and false, which is absurd, therefore what we originally supposed
    is wrong therefore B^C=>A is true.

    Thanks,

  6. Steve says:

    Understand this is not REALLY part of the problem, but wondered if you could prove uniqueness without bruteforce. I like the details of your proof, nice job! Thanks

  7. Enrique Rivera says:

    Java implementation:
    //——————————————————
    z:
    for (int i = 1;; i++) {
    for (int j = i + 1; j < 1000; j++) {
    double sqrt = Math.sqrt((i * i) + (j * j));
    if (sqrt % ((int)sqrt) == 0) {
    if (i + j + Math.sqrt((i * i) + (j * j)) == 1000) {
    System.out.println(i + " + " + j + " + " + ((int)Math.sqrt((i * i) + (j * j))) + " = " + 1000);
    break z;
    }
    }
    }
    }
    //——————————————————

  8. Samir says:

    First off, I wanted to say great blog. I found it while researching some algorithms for the Project Euler challenges I’m working on as I learn a new programming language.

    In looking at your recommendation in implementing the Brute Force approach, I understand how we can optimize the program by noting that a<s/3 (the smallest side has to be less than the average length of all of the sides) but how did you come about that b<s/2?

    Thanks!

  9. Kristian says:

    If a is really really small then b is almost as long as c, since they are almost parallel. Therefore b must be smaller than half the total circumference. Does that make sense?

  10. Samir says:

    Hi Kristian,

    I understand what you are saying but how can we be sure that a will be that small?

    Thanks again.

  11. Kristian says:

    We can’t be sure a is small. But this is a worst case. If a is larger than 1 then b will be even smaller. And thus this is just creating an upper bound.

  12. Jean-Marie Hachey says:

    Output of the algos and the sum of the triplet.
    Considering the following triplets:
    99, 440, 451;
    165, 396, 429;
    180, 385, 425;
    264, 315, 411.
    The sum of each of these triplets is 990.
    In these cases, the algo gives only the triplet with the smallest leg value. (Here 99).

  13. Jean-Marie Hachey says:

    Table 1
    Pythagorean triples with sums 12, 168, 990 and 1000.
    Re: Project Euler, Problem 9:
    A 1000 Pythagorean triplets

    http://img15.hostingpics.net/pics/332741PE9tableone.jpg

    ___

    Sources:

    1) Project Euler 9
    Kristian’s algorithm
    http://www.mathblog.dk/files/euler/Problem9.cs
    2) Microsoft Visual C# 2010 Express
    3) Microsoft Office Excel (2007)
    4) Pythagorean triples [including multiples] up to 2100
    http://www.tsm-resources.com/alists/trip.html

  14. No says:

    That’s a nice program, another approach would be using Herons formula and the triangle inequality.
    That is:
    ab/2=sqrt(500(500-a)(500-b)(500-sqrt(a^2+b^2)))
    ->
    a^2b^2/4=500(500-a)(500-b)(500-sqrt(a^2+b^2))
    This means that ab is divisible by 100.
    Using that, a+b>500 and a<250, b<(1000-a)/2, b<500, and then checking if a^2+b^2=c^2 should be quick.

  15. Nhan Do says:

    Start by solving the equation:
    a^2 + b^2 = c^2
    => (a+b)^2 – 2ab = c^2
    => (1000 – c)^2 = 2ab = c^2
    => 500000 = ab + 1000c (1)

    Hence, ab must divisible by 1000, let ab = 1000*k (k is an int)
    => b= 1000/a

    Substitute to (1):
    => 500000 = 1000*k + 1000 ( 1000 – a – b)
    => 500000 = 1000*k + 1000 ( 1000 – a – 1000*k /a)
    => a^2 – (k +500)a + 1000k = 0

    Using the quadratic formula for a:
    a= ( (k+500) + sqrt(D) ) /2

    All we have to see if sqrt(D) is an integer.
    Because D = k*k – 3000*k + 500*500 and D >0.
    We can see that k can only goes from 0 to 85 OR >2915
    but as a*b 2915

    So we only need to check for k from 1 to 85:
    Please see the code below:

    #include
    #include

    using namespace std;

    int main() {

    double delta;
    int k;
    for ( k = 1 ; k < 85 ; k++) {
    delta = sqrt( (double)(k*k – 3000*k + 500*500 ) );

    if ( delta – (int)delta < 0.00001)
    break;
    }

    int a = ( (k + 500) + delta ) /2 ;
    int b = 1000*k / a;
    int c = 1000 – (a + b);

    cout << "Product: " << a*b*c << endl;

    return 0;
    }

  16. Nhan Do says:

    Start by solving the equation:
    a^2 + b^2 = c^2
    => (a+b)^2 – 2ab = c^2
    => (1000 – c)^2 = 2ab = c^2
    => 500000 = ab + 1000c (1)

    Hence, ab must divisible by 1000, let ab = 1000*k (k is an int)
    => b= 1000/a

    Substitute to (1):
    => 500000 = 1000*k + 1000 ( 1000 – a – b)
    => 500000 = 1000*k + 1000 ( 1000 – a – 1000*k /a)
    => a^2 – (k +500)a + 1000k = 0

    Using the quadratic formula for a:
    a= ( (k+500) + sqrt(D) ) /2

    All we have to see if sqrt(D) is an integer.
    Because D = k*k – 3000*k + 500*500 and D >0.
    We can see that k can only goes from 0 to 85 OR >2915
    but as a*b2915

    So we only need to check for k from 1 to 85:
    Please see the code below:

    #include
    #include

    using namespace std;

    int main() {

    double delta;
    int k;
    for ( k = 1 ; k < 85 ; k++) {
    delta = sqrt( (double)(k*k – 3000*k + 500*500 ) );

    if ( delta – (int)delta < 0.00001)
    break;
    }

    int a = ( (k + 500) + delta ) /2 ;
    int b = 1000*k / a;
    int c = 1000 – (a + b);

    cout << "Product: " << a*b*c << endl;

    return 0;
    }

  17. Nhan Do says:

    I don’t know why the line got cut:
    “but as a*b 2915″

  18. Nhan Do says:

    Here is the missing line:
    “but as a*b is smaller than 1000^2, remember that a*b= 1000*k, k cannot greater than 2915″

Trackbacks For This Post

  1. Project Euler - Problem 9 Solution

Leave a Comment Here's Your Chance to Be Heard!

You can use these html tags: <a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <strike> <strong>

You can use short tags like: [code language="language"][/code]

The code tag supports the following languages: as3, bash, coldfusion, c, cpp, csharp, css, delphi, diff,erlang, groovy, javascript, java, javafx, perl, php, plain, powershell, python, ruby, scala, sql, tex, vb, xml

You can use [latex]your formula[/latex] if you want to format something using latex

Notify me of comments via e-mail. You can also subscribe without commenting.