0

It is clearly possible to compute sqrt(x*x + y*y) without using the hypot(x,y) function. Why is this function part of the standard?

Is this related to some internal (hardware based) optimizations?

nowox
  • 25,978
  • 39
  • 143
  • 293
  • 3
    All the functions in the `math.h` library can be programmed "by hand" with more or less difficulty, that is not a reason for it to not be part of the standard. `sqrt` can also be computed https://stackoverflow.com/a/60905298/6865932. – anastaciu Aug 10 '20 at 14:39
  • 1
    Good discussion [here](https://en.m.wikipedia.org/wiki/Hypot), with implementation technique and some references. – rici Aug 10 '20 at 14:45
  • 2
    In additional to avoiding undue overflows or underflows, `hypot(x,y)` may produce a more accurate result than doing it manually, alhough that's a quality of implemenation detail rather than part of the standard. For example, the IEEE-754 version of `hypot` in GNU Glibc claims the error in the return value is less than 1 ulps (units in the last place). – Ian Abbott Aug 10 '20 at 14:45
  • @rici Thanks for the link, yet ["Finally, the multiplication by |x| cannot underflow"](https://en.m.wikipedia.org/wiki/Hypot) fails in the corner case of `x,y` both subnormal and the product is inexact. the result can [underflow](https://stackoverflow.com/q/42277132/2410359). – chux - Reinstate Monica Aug 10 '20 at 15:05
  • @chux: AIUI, `|x|` is multiplied by the square root of a number between 1 and 2 (regardless of how that number was computed). That square root must therefore be between 1 and sqrt(2), and is not a subnormal. It's possible that |x| is subnormal, but how can multiplying a subnormal by a number between 1 and sqrt(2) underflow? The result is as precise as the subnormal was. Perhaps I don't understand the precise meaning of underflow. In any case, I doubt the wikipedians involved in that page are reading this thread, and I'm not one of them. Perhaps you should try the talk page. – rici Aug 10 '20 at 15:55
  • @rici 1) "how "can multiplying a subnormal by a number between 1 and sqrt(2) underflow? --> If the product is a subnormal and imprecise - in this case, `x` and `y` would both be subnormals. 2) I'll look into a wiki edit. – chux - Reinstate Monica Aug 10 '20 at 16:20
  • @chux: `y/x` is indeed problematic if both are subnormals, I agree. But the square root is still going to result in a number between 1 and sqrt(2). And the sentence you are talking about is about "the multiplication by `|x|`", which I take to mean to specifically be the multiplication of `|x|` by the square root, and not the computation of the square root (which has already happened at that point). So, `|x|` might be a subnormal. If it is, both are, becuase `|x|>=|y|`. The multiplication is `|x|` by a normal between 1 and 1.41 (approx.). Can that underflow? How? – rici Aug 10 '20 at 16:40
  • @chux: I only ask out of curiosity :-) I take no responsibility for the WP page, but I read it and it all seemed reasonable enough. Not that I'm a specialist in the area, or anything. (I also think that there are no problems if you happen to have 80-bit long doubles on your hardware, but I could be wrong about that, too.) – rici Aug 10 '20 at 16:43
  • @rici "Can that underflow? How?" --> if the product is an imprecise subnormal. "underflow" does not mean the product is zero, just less than the smallest _normal_ value and imprecise. Since `x` is a subnormal in this case, the product can also be-subnormal, even if it is larger than `x`. – chux - Reinstate Monica Aug 10 '20 at 16:48

3 Answers3

6

The hypot function performs this calculation while avoiding overflows.

Section 7.12.7.3p2 of the C standard regarding the hypot function states:

The hypot functions compute the square root of the sum of the squares of x and y, without undue overflow or underflow. A range error may occur.

dbush
  • 205,898
  • 23
  • 218
  • 273
  • Unless I'm missing something, [all glibc does](https://github.com/bminor/glibc/blob/92b9636/sysdeps/i386/fpu/e_hypot.S#L48_) is `sqrt(x*x + y*Y)`. – OrangeDog Aug 10 '20 at 15:19
  • Having a `double` overflow is quite hard isn't ? Does it mean this function is slower than doing the `sqrt` directly? – nowox Aug 10 '20 at 16:57
5

In addition to hypot(x,y) avoiding undo overflow (@bush) and potentially providing a more accurate answer (@Ian Abbot) than sqrt(x*x + y*y), a subtle corner case is defined.

If std library adheres to ANSI/IEEE 754, __STDC_IEC_559__ defined, then

hypot(±∞; y) returns +∞, even if y is a NaN. C17dr § F.10.4.3

sqrt(x*x + y*y) would not.

There might be a cost: Why hypot() function is so slow?. It is a quality of implementation concern.

Is this related to some internal (hardware based) optimizations?

@njuffa discusses this.

chux - Reinstate Monica
  • 143,097
  • 13
  • 135
  • 256
2

The value of the sqrt(x*x + y*y) doesn't have to be equal to the value of the hypot(x, y). Example for demonstration:

#include <stdio.h>
#include <math.h>
#include <float.h>

int main (void)
{
    double x = DBL_MAX/2, y = x;
    printf("x = %g, y = %g\n", x, y);
    printf("sqrt(x*x + y*y) = %g\n", sqrt(x * x + y * y));
    printf("hypot(x*x + y*y) = %g\n", hypot(x, y));
}

The hypot(x,y) gives the (approximately) correct result even if the arithmetical value of x*x + y*y is greater than DBL_MAX.

M. Nejat Aydin
  • 9,597
  • 1
  • 7
  • 17
  • Note that the compiler, compilation options, and the execution environment can affect the result of this code. One possible difference is if `sqrt(x*x + y*y)` is computed per the semantics of the abstract machine or if extended precision is allowable for intermediate results. In the first case, `x*x` can overflow, but might not if extended precision is used. See https://gcc.gnu.org/wiki/FloatingPointMath for some possibilities. – Andrew Henle Jan 19 '21 at 22:08