Efficiently calculating all pythagorean triples knowing the hypoteneuse

Question

Given a hypoteneuse (c in the typical equation a*a + b*b = c*c), what is an efficient way to calculate all possible integer values of a and b, such that a < b?

Note: I've seen c be greater than 1e12, thus c*c is greater than long.MaxValue, from what I can tell, c*c does fit into a decimal, though.

Answer 1

All pythagorean triples (a,b,c) satisfy the property that, for some integers k,m and n,

a=k(m^2-n^2), b=2kmn, c=k(m^2 + n^2)

So start by factoring c. Then for every distinct factor k of c (i.e. for every distinct subset of the factors, multiplied together), find all m and n that satisfy c/k = (m^2 + n^2). Doing the latter will take significantly less time than the brute-force approach others have suggested (you only have to find squares that sum to c/k, instead of c^2), but will identify all triples (a,b,c). You also don't need to use bignums, because all of the intermediate results fit into a long.

I also suggest you check out the web page http://www.maths.surrey.ac.uk/hosted-sites/R.Knott/Pythag/pythag.html under the heading "A more general Pythagorean Triple Calculator" which contains an embedded calculator, written in javascript, that does precisely what you want.

Answer 2

There is a mathematical solution that finds a and b fast even for large values of c. Unfortunately, it is not that simple. I'm trying to give a short explanation of the algorithm. I hope it is not too confusing.

Since c is given and you are looking for a and b, you basically want to solve diophantine equations of the form

n=x^2+y^2,

where n is given. It doesn't help much that n = c*c is a square and thus I'm describing a solution for any n. If n were a prime number then you could use Fermat's theorem, to decide if your equation is solvable, and use, as Moron pointed out the Hermite-Serret algoritm to find the solutions if there are any.

To solve the case where n is not prime it is a good idea to use Gaussian integers. (Gaussian integers are just complex numbers with integral coefficients). In particular one notes that the norm of x+yi is

N(x+yi) = x^2+y^2.

Hence one has to find the Gaussian integers x+yi whose norm is n. Since the norm is is multiplicative it is sufficient to solve this equation for the factors of n and then to multiply the Gaussian integers of the individal equations. Let me give an example. We want to solve

65 = x^2 + y^2.

This is equivalent to find x,y such that

N(x+yi) = 65

over the Gaussian integers. We factor 65 = 5 * 13, then we use Fermat to note that both 5 and 13 can be represented as sum of two squares. Finally, we find either by using brute force of by using the Hermite-Serret algorithm

N(2+i) = N(1+2i) = ... = 5
N(2+3i) = N(3+2i) = ... = 13

Note, I've left out the Gaussion integers 2-i, -2+i, etc with negative coefficients. Those are of course solutions too.

Hence we can now multiply these solutions together to get

65 = 5*13 = N(2+i) * N(2+3i) = N((2+i) * (2+3i)) = N(1 + 8i)

and

65 = 5 * 13 = N(2+i) * N(3+2i) = N((2+i) * (3+2i)) = N(4 + 7i).

Hence, one gets the two solutions

1*1 + 8*8 = 65
4*4 + 7*7 = 65

The other combinations e.g. with negative coefficients need to be checked too. They don't give new solutions other than permutations and changed signs.

---

To compute all the solutions one might also add that there is no need to ever compute c*c. Finding the factors of c is all that is necessary. Also since a and b are both smaller than c, it will not happen that products of Gaussian integers are not representable with 64-bit integer coefficients. Hence, if one is careful, 64-bit integer are enough precision to solve the problem. Of course, it is always easier to just use a language like Python that does not have this kind of overflow problems.