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.
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 algorithm 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 individual 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.