Finding Pythagorean triples

New to Julia
1

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?

2

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
3

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

5 Likes
4

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. :grinning:

5

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

6

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
7

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
8

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? :slight_smile:

3 Likes
9

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
10

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