[Bug 253313] lib/msun: hypotl(3) mishandles subnormal numbers

[bugzilla-noreply at freebsd.org]

https://bugs.freebsd.org/bugzilla/show_bug.cgi?id=253313

--- Comment #2 from Dimitry Andric <dim at FreeBSD.org> ---
The most minimal fix I could come up with is:

diff --git a/lib/msun/src/e_hypotl.c b/lib/msun/src/e_hypotl.c
index 9189b1fab54d..c66d2246c8e2 100644
--- a/lib/msun/src/e_hypotl.c
+++ b/lib/msun/src/e_hypotl.c
@@ -82,8 +82,8 @@ hypotl(long double x, long double y)
                man_t manh, manl;
                GET_LDBL_MAN(manh,manl,b);
                if((manh|manl)==0) return a;
-               t1=0;
-               SET_HIGH_WORD(t1,ESW(MAX_EXP-2));       /* t1=2^(MAX_EXP-2) */
+               /* t1=2^(MAX_EXP-2) (or maybe just t1=0x1p+16382 ?) */
+               INSERT_LDBL80_WORDS(t1,ESW(MAX_EXP-2),0x8000000000000000);
                b *= t1;
                a *= t1;
                k -= MAX_EXP-2;
diff --git a/lib/msun/src/math_private.h b/lib/msun/src/math_private.h
index b91b54cea689..200a734f1233 100644
--- a/lib/msun/src/math_private.h
+++ b/lib/msun/src/math_private.h
@@ -272,7 +272,7 @@ do {                                                       
        \

 #define        INSERT_LDBL80_WORDS(d,ix0,ix1)                          \
 do {                                                           \
-  union IEEEl2bits iw_u;                                       \
+  volatile union IEEEl2bits iw_u;                              \
   iw_u.xbits.expsign = (ix0);                                  \
   iw_u.xbits.man = (ix1);                                      \
   (d) = iw_u.e;                                                        \

Giving as output:

Via scaling and sqrtl
x=2.853684e-4932 y=1.444012e-4922 a=1.444012e-4922
x=0x1.b2933cafa0bb7p-16383 y=0x1.fffffffffffffp-16351
a=0x1.000000006ca4c72cp-16350

Via hypotl
x=2.853684e-4932 y=1.444012e-4922 z=1.444012e-4922
x=0x1.b2933cafa0bb7p-16383 y=0x1.fffffffffffffp-16351
z=0x1.fffffffffffffp-16351

This consists of two parts: the first being that the "t1=0;
SET_HIGH_WORD(t1,ESW(MAX_EXP-2));" statements result in a long double that
printf's as 0x1p+16382, but has the high part of its mantissa set to 0. This is
different from the 'real' 0x1p+16382 constant, which has high part of the
mantissa set to 0x80000000 instead. E.g. compare this with glibc's
implementation
(https://sourceware.org/git/?p=glibc.git;a=blob;f=sysdeps/ieee754/ldbl-96/e_hypotl.c):

    85          if(__builtin_expect(eb < 0x20bf, 0)) {  /* b < 2**-8000 */
    86              if(eb == 0) {       /* subnormal b or 0 */
    87                  uint32_t exp __attribute__ ((unused));
    88                  uint32_t high,low;
    89                  GET_LDOUBLE_WORDS(exp,high,low,b);
    90                  if((high|low)==0) return a;
    91                  SET_LDOUBLE_WORDS(t1, 0x7ffd, 0x80000000, 0); /*
t1=2^16382 */
    92                  b *= t1;
    93                  a *= t1;
    94                  k -= 16382;

The second part is the addition of volatile to the INSERT_LDBL80_WORDS macro.
This is basically a hack to force the compile to not optimize away the
undefined behavior or reading and writing from different union fields. It
should really be fixed more properly, and for all these macros, but that is
much more invasive.

(Note that the second part isn't needed for gcc, as it appears to avoid
optimizing around this type of union field access.)

-- 
You are receiving this mail because:
You are the assignee for the bug.