Bug 27460 - hypot(__DBL_DENORM_MIN__, -0.) incorrectly raises FE_INEXACT and FE_UNDERFLOW
Summary: hypot(__DBL_DENORM_MIN__, -0.) incorrectly raises FE_INEXACT and FE_UNDERFLOW
Status: RESOLVED NOTABUG
Alias: None
Product: glibc
Classification: Unclassified
Component: math (show other bugs)
Version: 2.31
: P2 normal
Target Milestone: ---
Assignee: Not yet assigned to anyone
URL:
Keywords:
Depends on:
Blocks:
 
Reported: 2021-02-23 16:07 UTC by Matthias Kretz (Vir)
Modified: 2021-04-08 10:34 UTC (History)
2 users (show)

See Also:
Host:
Target: x86_64-pc-linux-gnu
Build:
Last reconfirmed:


Attachments

Note You need to log in before you can comment on or make changes to this bug.
Description Matthias Kretz (Vir) 2021-02-23 16:07:29 UTC
Testcase (https://godbolt.org/z/b4Ko4Y):

#include <fenv.h>
#include <math.h>

int main()
{
    volatile double a = __DBL_DENORM_MIN__;
    volatile double c = -0.;
    feclearexcept(FE_ALL_EXCEPT);
    //c = hypot(a, c); // this argument order raises no exceptions
    c = hypot(c, a);
    return fetestexcept(FE_ALL_EXCEPT);
}

This program raises FE_INEXACT and FE_UNDERFLOW. Since one argument to hypot is 0 the result is exact: the absolute value of the other argument. If you switch the argument order the program raises no floating-point exceptions.
Comment 1 Adhemerval Zanella 2021-02-23 20:13:25 UTC
Afaik since we define __STDC_IEC_559__ POSIX allows the implementation to return an error range [1].  This is in fact what the fix for BZ#18803 does, where it explicit raise underflow exception when the result is tiny and inexact.

I think we are in fact missing a underflow check for the case where second argument is '0':

diff --git a/sysdeps/ieee754/dbl-64/e_hypot.c b/sysdeps/ieee754/dbl-64/e_hypot.c
index 9ec4c1ced0..bcd172130e 100644
--- a/sysdeps/ieee754/dbl-64/e_hypot.c
+++ b/sysdeps/ieee754/dbl-64/e_hypot.c
@@ -100,7 +100,10 @@ __ieee754_hypot (double x, double y)
          uint32_t low;
          GET_LOW_WORD (low, b);
          if ((hb | low) == 0)
-           return a;
+           {
+             math_check_force_underflow (x);
+             return a;
+           }
          t1 = 0;
          SET_HIGH_WORD (t1, 0x7fd00000);       /* t1=2^1022 */
          b *= t1;


[1] https://pubs.opengroup.org/onlinepubs/9699919799/functions/hypot.html
Comment 2 Matthias Kretz (Vir) 2021-02-23 20:20:01 UTC
From C17:

> F.10.4.3
>  The hypot functions
> — hypot(x, y), hypot(y, x), and hypot(x, −y) are equivalent.
> — hypot(x, ±0) is equivalent to fabs(x).
> — hypot(±∞, y) returns +∞, even if y is a NaN.

Combining the first two bullets leads to: hypot(±0, x) is equivalent to fabs(x).
Comment 3 Adhemerval Zanella 2021-02-23 20:37:15 UTC
C11 already states these rules and I am not sure if the rules are transitive as you indicate. And I would not assume the rules are transitive becuase for 'pow' (F.10.4.4 The pow functions) the standard states all possible argument combinations.

Not assuming they are transitive means that my previous suggestion is incorrect and current code already follow the "hypot(x, ±0) is equivalent to fabs(x)" rule since hypot (__DBL_DENORM_MIN__, 0) returns __DBL_DENORM_MIN__.
Comment 4 Matthias Kretz (Vir) 2021-02-23 20:46:35 UTC
I had C17 handy...

My test case shows the current code violates "hypot(x, y), hypot(y, x), and hypot(x, −y) are equivalent". I don't think transitivity is the correct term here. All these items must be true for hypot to be conforming to Annex F.
Comment 5 jsm-csl@polyomino.org.uk 2021-02-23 20:54:24 UTC
I think the rules in Annex F allowing spurious "inexact" and "underflow", 
"unless stated otherwise", apply here (glibc limits spurious underflow to 
cases where the returned result would be possible for an underflowing 
computation).  (I consider the rules about functions bound to IEC 60559 
operations, and about functions with qNaN argument returning qNaN with no 
exceptions, to be "stating otherwise", but not the rules stating the value 
of a result for a non-NaN argument.)
Comment 6 Matthias Kretz (Vir) 2021-02-23 21:06:29 UTC
F10/9 Whether or when library functions raise an undeserved "underflow" floating-point exception is unspecified.379) Otherwise, as implied by F.8.6, the <math.h> functions do not raise spurious floating-point exceptions (detectable by the user), other than the "inexact" floating-point exception.

Footnote 379) It is intended that undeserved "underflow" and "inexact" floating-point exceptions are raised only if avoiding them would be too costly.

Additionally I believe F.10.4.3 states otherwise. I.e. hypot(denorm_min, 0) is equivalent to fabs(denorm_min) and hypot(0, denorm_min) is equivalent to hypot(denorm_min, 0).

But I'm new to interpreting floating-point exceptions legalese. I was only trying to verify my hypot SIMD implementation against yours and stumbled over this inconsistency.
Comment 7 Adhemerval Zanella 2021-02-24 13:32:11 UTC
Joseph, do you considere this a real issue for glibc (meaning that we should keep this bug opened) or should we close it?
Comment 8 jsm-csl@polyomino.org.uk 2021-02-24 22:48:56 UTC
I don't think there is a bug for glibc.  Maybe there is an issue with the 
C standard wording that it should say that hypot *returns* a value given 
by fabs rather than *is equivalent* to a given call to fabs (saying a 
function not bound to an IEEE operation is equivalent to one bound to such 
an operation is very questionable).
Comment 9 Paul Zimmermann 2021-04-08 04:38:24 UTC
see https://mailman.oakapple.net/pipermail/cfp-interest/2021-April/001966.html: full support for the inexact exception is not required for hypot (and mathematical functions in general)
Comment 10 Adhemerval Zanella 2021-04-08 10:34:10 UTC
Closed based on both Joseph and Paul comments.