Why is sqrt(x*x + y*y) != math.hypot(x, y) in Python 3.8?

Question

I'm trying to run some tests on the shiny new Python 3.8 and noticed an issue with math.hypot. From the docs:

For a two dimensional point (x, y), this is equivalent to computing the hypotenuse of a right triangle using the Pythagorean theorem, sqrt(x*x + y*y).

However, these are not equivalent in 3.8:

>>> from math import hypot, sqrt
>>> x, y = 95, 168
>>> sqrt(x*x + y*y), hypot(x, y), sqrt(x*x + y*y) == hypot(x, y)
(193.0, 193.00000000000003, False)
>>> sqrt(x*x + y*y).is_integer(), hypot(x, y).is_integer()
(True, False)

In 3.7 both ways produce exactly the same result ("193.0", which is considered an integer).

Answer 1

The function hypot offers another approximation of the mathematical expression √(x² + y²), just like the floating-point expression sqrt(x*x + y*y) is an approximation of this same mathematical expression.

The function hypot is recommended because it solves very noticeable defects that are present in the floating-point computation sqrt(x*x + y*y) with very large or small values. For instance, if x is only a bit larger than the square root of the maximum finite floating-point value, sqrt(x*x + y*y) always produces +inf because x*x produces +inf.

Compare:

>>> x, y = 95E200, 168E200
>>> sqrt(x*x + y*y), hypot(x, y)
(inf, 1.93e+202)
>>> z, t = 95E-200, 168E-200
>>> sqrt(z*z + t*t), hypot(z, t)
(0.0, 1.93e-198)

For these two (respectively very large and very small) pairs of inputs, hypot is doing fine, whereas sqrt(x*x + y*y) is catastrophically wrong.

---

When the naïve version sqrt(x*x + y*y) works reasonably well (when the values x and y are neither very large nor very small), it can be may be more or less accurate than the function hypot depending on the values of x and y. They can both be expected to produce a result that is a few ULPs away from the mathematical result. But since they are different approximations obtained by different methods, they may differ (in the worst case by twice “a few ULPs”).

One typical implementation for hypot(x, y) is first to swap x and y if necessary so that x has the largest magnitude, and then compute x * sqrt(1 + (y/x)*(y/x)). This solves the problem with x*x overflowing. As a side-effect, it means that even when there is no overflow, the result is slightly different from sqrt(x*x + y*y).

Note that it's normal that sqrt(x*x + y*y) is more precise when you apply it to small integers (as you do in your test): when x and y are small integers, x*x and y*y and their sum can be computed exactly as floating-point values. If this sum is the square of an integer, the floating-point function sqrt can only compute this integer. In short, in this scenario the computations, despite being floating-point, are exact from beginning to end. In contrast, the typical hypot implementation above starts by computing x/y (in your test, 95.0/168.0), and this result is not in general representable exactly as a floating-point value. The first step already incurs an approximation, and this approximation can result in the final result being wrong (as it is in your test)!

---

There is no standard algorithm for hypot: it is only expected to compute a good approximation of the mathematical expression √(x² + y²) while avoiding the overflow and underflow problems. This article shows different implementations, and points out that the popular implementation that I mentioned sacrifices accuracy to avoid overflow and underflow (but the article also provides a floating-point implementation for hypot that is more accurate than sqrt(x*x + y*y) even where sqrt(x*x + y*y) works).