3

hypot() function returns infinity for a = b = 1e154 on OSX. It returns a valid value on Ubuntu.

#include <stdio.h>
#include <math.h>
int main() {
  long double h1 = hypot(1e153,1e153);
  long double h2 = hypot(1e154,1e154);
  long double h3 = hypot(1e155,1e155);
  printf("h1 = %Le\n", h1);
  printf("h2 = %Le\n", h2);
  printf("h3 = %Le\n", h3);
  return 0;
}

Result:

h1 = 1.414214e+153
h2 = inf
h3 = 1.414214e+155

It's a rare case but it overflows way under the maximum value (1e308 something)... Is there any way to avoid it? It overflows only on OSX.

Jihyun
  • 883
  • 5
  • 17
  • 4
    Use [`hypotl()`](https://pubs.opengroup.org/onlinepubs/9699919799/functions/hypot.html)? – pmg Dec 03 '21 at 18:12
  • 2
    *way under the maximum value* Are you sure about that? What is `x*x` when `x = 1e155`? – Andrew Henle Dec 03 '21 at 18:16
  • thank you... I will try hypotl() but I am calling hypot() from another language where long double (128bit) is not officially supported.... I know when x is 1e155, x*x is 1e310 ... but hypot() is working fine up to 1e308... – Jihyun Dec 03 '21 at 18:23
  • 2
    Can't replicate with MSVC, even though it implements `long double` as `double`. Can't replicate using `double` either. Note that the link from @pmg says "These functions shall compute the value of the square root of x²+ y² *without undue overflow or underflow.* – Weather Vane Dec 03 '21 at 18:39
  • 1
    Evidently, IEEE-754 recommends a `hypot()` implementation that is resistant to unnecessary overflows. It looks like glibc and and the MS runtime provide such an implementation, but whatever math library the OP is using does not. An alternative would be to implement one's own. Wikipedia offers some [algebraic variants](https://en.wikipedia.org/wiki/Pythagorean_addition#Implementation_2) on the naive formula that are more resistant to overflow. – John Bollinger Dec 03 '21 at 18:44
  • 2
    Sure looks like a bug in the MacOS libc. (I've confirmed it on my Mac.) – Steve Summit Dec 03 '21 at 19:29
  • 1
    I filed a bug report and emailed the manager of Apple’s math group. – Eric Postpischil Dec 04 '21 at 00:00
  • 1
    For a workaround, you might write a C routine with another name that takes `double` arguments (so it is compatible with the other language you are using), calls `hypotl`, and returns a `double` result. If you need strict semantics regarding rounding or exceptions, update the question with the requirements, and people can suggest additional mitigations. – Eric Postpischil Dec 04 '21 at 00:02
  • thank you everyone. I used a hypot() function found in https://stackoverflow.com/questions/3764978/why-hypot-function-is-so-slow (A simple FMA-based double-precision implementation one) and used it in my code. It is working fine. thank you again – Jihyun Dec 04 '21 at 02:01

0 Answers0