Incorrect overflow in std::hypot

[efriedma-quic]

Consider:

#include <cmath>
float f = std::hypot(1e20f,1e20f,1e20f);

With libc++, this produces +inf. Other libraries appear to handle this correctly. (The C standard suggests hypot() should return a result "without undue overflow or underflow"... although that only applies to the two-argument form).

[PaulXiCao]

When fixing this issue one could also look at the symmetrical problem for underflows, e.g.

#include <cmath>
float f = std::hypot(1e-30f, 1e-30f, 1e-30f);

This yields 0 for libc++ and the correct non-zero result for libstdc++. See this godbolt which shows clang using both standard libraries.

[PaulXiCao]

The current implementation uses a literal translation of the mathematical definition, e.g.

float hypot(float a, float b, float c) { return std::sqrt(a*a + b*b + c*c); }

I would propose a version which scales the arguments. Similar to the following:

#include <algorithm> // std::max
float hypot(float a, float b, float c) { 
  // ... todo: handle infs, NaNs
  const float M = std::max({std::abs(a), std::abs(b), std::abs(c)});
  if (M == 0)
    return 0;
  else
    return M * std::sqrt( std::pow(a/M, 2) + std::pow(b/M, 2) + std::pow(c/M, 2) );
}

Note, more thought needs to be put into handling input validation for infinities and/or NaNs.

Futhermore, one could take a look at new possibilities for underflowing of intermediate terms, e.g. hypot(1e-15f, 1e15f, 1e15f) for which a*a did not underflow, but std::pow(a/M, 2) does.

[efriedma-quic]

See also https://github.com/llvm/llvm-project/blob/main/libc/src/math/generic/hypotf.cpp

e.g. hypot(1e-15f, 1e15f, 1e15f) for which a*a did not underflow, but std::pow(a/M, 2) does.

The underflow doesn't really matter in this case: hypot(0, 1e15f, 1e15f) produces the same result. It's only relevant if all the operands are small.

[efriedma-quic]

See also https://github.com/llvm/llvm-project/blob/main/libc/src/math/generic/hypotf.cpp

See also https://github.com/llvm/llvm-project/blob/main/libc/src/__support/FPUtil/Hypot.h

[lntue]

   return M * std::sqrt( std::pow(a/M, 2) + std::pow(b/M, 2) + std::pow(c/M, 2) );

This scaling is expensive and incurs extra rounding, also I think it's better to not call std::pow. To prevent overflow and underflow, we can use some fixed scaling factor that is exact to perform:

if (M > 0x1.0p511) {
  // a*a + b*b + c*c might overflow double precision
  a *= 0x1.0p-600;  // scale down
  b *= 0x1.0p-600;  // scale down
  c *= 0x1.0p-600;  // scale down
  M = 0x1.0p600; // scale up factor
} else if (M < 0x1.0p-511) {
  // a*a, b*b, c*c might underflow double precision
  a *= 0x1.0p600;  // scale up
  b *= 0x1.0p600;  // scale up
  c *= 0x1.0p600;  // scale up
  M = 0x1.0p-600; // scale down factor
} else {
  M = 1.0;
}
return M * sqrt(a*a + b*b + c*c);

[PaulXiCao]

I will try to create a patch similar to your proposed version. For that purpose I will try to give an explanation for the "magic-numbers".

A double precision IEEE 754 has a range for exponents of [-1022, 1023]. (Biased 11-bit exponent.) The numbers a,b,c squared should not exceed those bounds. This explains 0x1.0p+-511.

The scaling factor needs to be between +-511 and the exponent bounds. This explains 0x1.0p600.

Similar arguments can be made for floats and long double.