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

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,a^{2}+ b^{2}= c^{2}

For example, 3^{2}+ 4^{2}= 9 + 16 = 25 = 5^{2}.

There exists exactly one Pythagorean triplet for whicha+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 a^{n} + b^{n} = c^{n} 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 =m^{2}-n^{2}

b =2m*n

c =m^{2}+n^{2}

forms a Pythagorean triplet. You can convince your self that (m^{2}-n^{2})^{2} + (2m*n)^{2} = (m^{2}+n^{2})^{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: PrimitiveA Pythagorean Triplet is called primitive if

a,bandcare 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 d^{2} 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 m^{2}. Since the prime is also a divisor of c = m^{2} + n^{2}, then it follows that it must divide n^{2} 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(m^{2}-n^{2})

b =2d*m*n

c =d(m^{2}+n^{2})

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(m^{2}-n^{2}) + 2d*m*n + d(m^{2}+n^{2}) = 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

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.

my previous post should read

A B ^ C

where means if and only if

Arghhh. The characters are not showing in the post.

Let’s just say that part should read

A if and only if B ^ C

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?

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,

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

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;

}

}

}

}

//——————————————————

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!

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?

Hi Kristian,

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

Thanks again.

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.

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

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

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.

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;

}

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;

}

I don’t know why the line got cut:

“but as a*b 2915”

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”

Why is k an odd number?

Why is k an odd number?

When I write my code and did not check k for odd number, I got duplicate set:

findtriplets(240):

d = 20; m = 2; n = 1; a = 60; b = 80; c = 100;

d = 10; m = 3; n = 1; a = 80; b = 60; c= 100;

d = 8; m = 3; n = 2; a = 40; b = 96; c = 104;

d = 4; m = 5; n = 1; a = 96; b = 40; c = 104;

….

This problem is highly addictive:) Here is another solution.

Examine each divisor of s / 2, d

Find m and n that

m > 0

n > 0

m + n is an odd number

m and n are coprime

m * (m + n) = d

m > n

If found,

a = (s / 2 / d) * (m^2 – n^2)

b = ( s / 2 / d) * (2 * m * n)

c = (s / 2 / d) * ( m^2 + n ^2)

The only thing I am unsure of is (in the brute force method = BFI brute force and ignorance); how do we know that a < s/3 and b < s/2? My solution (nested loops) using QB64 (no cracks please; I'm over 60 and I fell in love with qbasic & quickbasic decades ago). I did love quick C but; that was back in the early 90's. I really should find something like it to wake up some old memory cells…gsw

@GSW

Since we know that a + b + c = 1000, and that a < b < c. Then we can deduce that

a < 1000 / 3 because if 'a' were not less than 333.333, then there would not be enough left for 'b' and 'c' to be greater. If 'a' was equal to 333.333 then there would only be 666.666 left for 'b' and 'c' to be.

The same can be deduced for b < 1000 / 2. If 'b' were greater than 500, then that means 'c' would have to be less than 500 cause thats all thats left, and that would break the rule

a < b < c.

I don’t write code. I worked on this using pencil and paper:

perimeter=a+b+c=(m^2-n^2)+2mn+(m^2+n^2)=1000. This gave me

2m(m+n)=1000, which gave m(m+n)=500

This can be factored 1(1+499), 2(2+248), 4(4+121), 5(5+95),

10(10+40), 20(20+5). Only in the last case is m^2-n^2 positive. I tried this with “perimeter 990 triplets” (see Jean-Marie Hachey’s post dated 11/20/2013), but couldn’t solve for Jean-Marie’s triplets. None of these triplets are primitive, yet I don’t know why I was able to solve for perimeter 1000, but not for perimeters 990. Is this because there is only one Pythagorean triplet with perimeter 1000?

Hi! Nice blog!

Here are my 2 cents about brute force solution:

1.

“… or we could use the facts that a < b < c, and thus exploit that a < s/3, and a < b < s/2."– Why b < s/2, I supposeb <(s-a)/2

2. Nested loop can be terminated on a*a+b*b > c*c

3. Nested loop can be replaced with binary search:

start=a; end=

~~/2~~while start c*c: end=b

else: start=b+1

> and thus exploit that a < s/3, and a < b < s/2.

I'm so sorry, but can anybody explain me this proposition? I don't understand this.

thx

Here is an answer in java. It is a bad brute force example, but it works. It outputs a, b, and c plus the answer. Enjoy!

public class pythagoreanTriplets{

//Variables

private static int a = 1000;

private static int b = 1000;

private static int c = 1000;

private static boolean found = false;

public static void main(String args[]){

while (!found){

if(a + b + c == 1000 && a < b && b < c){

if((a*a) + (b*b) == c*c){

found = true;

System.out.println(a*b*c);

System.out.println("A: " + a + " B: " + b + " C:" + c);

}

}

c –;

if(c == 0){

c = 1000;

b –;

}

if(b == 0){

b = 1000;

a –;

}

}

}

}

Hi, I have a question about the final code/algorithm. In the outer for loop, can you explain why “(s/2) % m == 0” means that we have found m? Since s/2 = m(m+n)d, I understand that m will always be found before (m+n), so we don’t have to worry that we have actually found k instead of m, but how do we know for sure that we haven’t found d (which is gcd(a,b,c)) before m? Isn’t it possible that d is smaller than m?

Actually there’s an easier solution which solves the problem in

O(n).we have N, that is a+b+c;

Therefore a+b+c = N

c = N-a-b……………[1]

since a^2 + b^2 = c^2, (let’s call b^2 as b2 for simplicity)

b2 = c2 – a2

b2 = (N2 – (a+b)2) – a2 ………. squaring equation [1]

solving the equation further , you get

b = N^2 – 2Na / (2N – 2a)

and c = N – a – b

Java code :

public static Number doStuff(int num){

int a, b, c;

int lim1 = num/2;

int lim2 = num/3;

outer:

for(a=1; a<lim1; a++){

for(b=1; b<lim2; b++){

c = num - a - b;

if(isTriplet(a, b, c)){

return (new Solution().new Number(a,b,c));

}

}

}

`return (new Solution().new Number(-1, -1, -1));`

}

Nice analysis

Here’s my groovy code:

`def a = {n, m -> n * n - m * m}`

def b = {n, m -> 2 * n * m}

def c = {n, m -> n * n + m * m}

def pythagorean = { n, m -> [a(n, m), b(n, m), c(n, m)] }

`int limit = Math.sqrt(1000/2)`

`def n = (limit..1).findResult { n->`

def searchM = (1..n).find{ m -> pythagorean(n, m).sum() == 1000 }

searchM?[n,searchM]:null

}

`println pythagorean(n).inject{ x, y -> x * y }`