Date: Sun, 07 Feb 2021 16:17:17 +0000 From: bugzilla-noreply@freebsd.org To: bugs@FreeBSD.org Subject: [Bug 253313] lib/msun: hypotl(3) mishandles subnormal numbers Message-ID: <bug-253313-227-NPLfMBYbut@https.bugs.freebsd.org/bugzilla/> In-Reply-To: <bug-253313-227@https.bugs.freebsd.org/bugzilla/> References: <bug-253313-227@https.bugs.freebsd.org/bugzilla/>
next in thread | previous in thread | raw e-mail | index | archive | help
https://bugs.freebsd.org/bugzilla/show_bug.cgi?id=3D253313 --- Comment #2 from Dimitry Andric <dim@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)=3D=3D0) return a; - t1=3D0; - SET_HIGH_WORD(t1,ESW(MAX_EXP-2)); /* t1=3D2^(MAX_EXP-= 2) */ + /* t1=3D2^(MAX_EXP-2) (or maybe just t1=3D0x1p+16382 ?) */ + INSERT_LDBL80_WORDS(t1,ESW(MAX_EXP-2),0x8000000000000000); b *=3D t1; a *=3D t1; k -=3D 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 {=20=20=20=20=20=20=20=20=20=20=20=20=20=20=20=20=20= =20=20=20=20=20=20=20=20=20=20=20=20=20=20=20=20=20=20=20=20=20=20=20=20=20= =20=20=20=20=20=20=20=20=20=20=20=20=20 \ #define INSERT_LDBL80_WORDS(d,ix0,ix1) \ do { \ - union IEEEl2bits iw_u; \ + volatile union IEEEl2bits iw_u; \ iw_u.xbits.expsign =3D (ix0); \ iw_u.xbits.man =3D (ix1); \ (d) =3D iw_u.e; \ Giving as output: Via scaling and sqrtl x=3D2.853684e-4932 y=3D1.444012e-4922 a=3D1.444012e-4922 x=3D0x1.b2933cafa0bb7p-16383 y=3D0x1.fffffffffffffp-16351 a=3D0x1.000000006ca4c72cp-16350 Via hypotl x=3D2.853684e-4932 y=3D1.444012e-4922 z=3D1.444012e-4922 x=3D0x1.b2933cafa0bb7p-16383 y=3D0x1.fffffffffffffp-16351 z=3D0x1.fffffffffffffp-16351 This consists of two parts: the first being that the "t1=3D0; 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. Thi= s 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=3Dglibc.git;a=3Dblob;f=3Dsysdeps/ieee754/ldb= l-96/e_hypotl.c): 85 if(__builtin_expect(eb < 0x20bf, 0)) { /* b < 2**-8000 */ 86 if(eb =3D=3D 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)=3D=3D0) return a; 91 SET_LDOUBLE_WORDS(t1, 0x7ffd, 0x80000000, 0); /* t1=3D2^16382 */ 92 b *=3D t1; 93 a *=3D t1; 94 k -=3D 16382; The second part is the addition of volatile to the INSERT_LDBL80_WORDS macr= o. 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.) --=20 You are receiving this mail because: You are the assignee for the bug.=
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?bug-253313-227-NPLfMBYbut>