# Finding Pythagorean triples

I was writing a small routine to find Pythagorean triples. Getting an x^2 + y^2 list was easy, but testing for sqrt of that gives me a Float rather than an Int, and floats are subject to binary error so I might not have an actual perfect square, or so it seems to me. Is there a way around this or am I thinking wrong?

As long as the number in question doesnât have any digits after the comma (i.e. they look like xyz.0) and the integer in question is <= 2^53 - 1 (the largest integer representable by a Float64 - yes, that is a JavaScript documentation, as JavaScript by default only has Float64 numbers for all number types), there is no error in the calculation. In other words, integers in [-(2^53 - 1), (2^53 - 1)] are representable by Float64 perfectly fine.

How large are your numbers going to be?

3 Likes

you want isqrt, ie the integer square root. You can check that isqrt(x^2+y^2)^2=x^2+y^2

5 Likes

Not That large. I only have an older computer and Iâd be sitting a long time taking a double loop for every integer up that far.

All these functions I donât know. Someday I have to go through the entire function list.

The problem isnât integer representability or some fuzzy âbinary errorâ â itâs just that floating point values round to the nearest representable value. So, for floating point square root to be âsafelyâ an integer, the question isnât just if the root is within [-(2^{53} - 1), (2^{53} - 1)]âŚ the danger is if square root number could be so close to a perfect square that it rounds to an integer.

In other words, could there possibly be an integer x such that (x+\epsilon)^2 is an integer with \text{abs}(\epsilon) < \frac{\text{eps}(x)}{2}? This would require 2x\epsilon \ge 1, and that could happen for any x \ge 2^{26}. Or squares on the order of 2^{52}.

Would be interesting to try to find such a square. I wouldnât be surprised if they exist.

7 Likes

You may already be aware of this, but just in case, you can generate all Pythagorean triples more efficiently by using Euclidâs formula:
a=k.(m^2-n^2)
b=2kmn
c=k.(m^2+n^2)
where to ensure you produce each triplet just once you need m, n, and k to be positive integers with m > n , and with m and n coprime and not both odd.

2 Likes

Sure enough, it exists:

julia> function find_next_imperfect_square(n)
while true
root = ân
if isinteger(root) && Int(root)^2 != n
return n
end
n += 1
end
end
find_next_imperfect_square (generic function with 1 method)

julia> n = find_next_imperfect_square(2^52)
4503599627370497

julia> root = ân
6.7108864e7

julia> isinteger(root)
true

julia> Int(root)^2
4503599627370496

julia> Int(root)^2 == n
false


So thereâs the true bound: isintegerâsqrt will work for all n < 4503599627370497. Didnât take much time at all to find. Know what that number is?

3 Likes

Using 1 as the y value is an easy way to convert these things into false Pythagorean triples. So the next thing to check would be if squaring the result gives x^2 + y^2 back. In this case it doesnât:

julia> sqrt(2^52 + 1^2)^2 == float(2^52 + 1^2)
false


However:

julia> isinteger(sqrt(2^54 + 1^2))
true

julia> sqrt(2^54 + 1^2)^2 == float(2^54 + 1^2)
true


Fortunately, in Julia mixed type comparisons donât work by converting both values to float and then comparing them, and this comparison is false because they donât represent the same value:

julia> sqrt(2^54 + 1^2)^2 == 2^54 + 1^2
false


So that would work because Juliaâs mixed type equality comparisons are very carefulâthis would fail in C, C++ or most languages to be honest.

On the other hand, doing it with purely integer arithmetic using isqrt is safe and likely to be faster:

function is_triple(x::Integer, y::Integer)
zÂ˛ = x^2 + y^2
zÂ˛ == isqrt(zÂ˛)^2
end

4 Likes

Good idea, although not on my old refurb computer from Amazon. Maybe when I get access to a supercomputer.