From https://github.com/JuliaMath/openlibm/issues/224, via Steve Kargl: Test program: ====================================================================== #include <stdio.h> #include <math.h> #if defined(__i386__) #include <ieeefp.h> #endif int main() { long double x, y, z, a; #if defined(__i386__) fpsetprec(FP_PE); #endif x = 0x3.6526795f4176ep-16384L; y = 0x3.ffffffffffffep-16352L; z = hypotl(x, y); if (x > y) { a = y / x; a = x * sqrtl(1 + a); } else { a = x / y; a = y * sqrtl(a + 1); } printf("Via scaling and sqrtl\n"); printf("x=%Le y=%Le a=%Le\n", x, y, a); printf("x=%La y=%La a=%La\n", x, y, a); printf("\nVia hypotl\n"); printf("x=%Le y=%Le z=%Le\n", x, y, z); printf("x=%La y=%La z=%La\n", x, y, z); return 0; } ====================================================================== Result on FreeBSD 14.0-CURRENT #0 main-n244563-d6f9c5a6d2f8 amd64: % cc -O2 -fno-builtin z.c -o z-amd64 -lm && ./z-amd64 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=nan x=0x1.b2933cafa0bb7p-16383 y=0x1.fffffffffffffp-16351 z=nan On FreeBSD 14.0-CURRENT #0 main-n244563-d6f9c5a6d2f8 i386: % cc -O2 -fno-builtin z.c -o z-i386 -lm && ./z-i386 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=nan x=0x1.b2933cafa0bb7p-16383 y=0x1.fffffffffffffp-16351 z=nan The expectation is that the results should be roughly similar, but at least hypotl should not produce nan in this case. With gcc 10.2 and not using -fno-builtin, the resulting assembly doesn't use *any* calls to libm, and produces: 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 clang 11.0 apparently doesn't have a built-in hypotl, so the output isn't different.
Another interesting detail, if I recompile libm's e_hypotl.c with gcc 10.2 and link it into libm.so, the result from the earlier z-amd64 executable becomes: % LD_LIBRARY_PATH=/usr/obj/usr/src/amd64.amd64/lib/msun /share/dim/bugs/bug253313/z-amd64 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=0.000000e+00 x=0x1.b2933cafa0bb7p-16383 y=0x1.fffffffffffffp-16351 z=0x0p+0 so now hypotl doesn't result in nan, but zero isn't right either, I guess... Interestingly the glibc and msun hypotl code seem to have the same Sun Microsystems origin, but have diverged quite a bit. Hard to tell from a simple viewing what causes this difference.
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.)
A commit in branch main references this bug: URL: https://cgit.FreeBSD.org/src/commit/?id=d3338f3355a612cf385632291f46c5777bba8d18 commit d3338f3355a612cf385632291f46c5777bba8d18 Author: Dimitry Andric <dim@FreeBSD.org> AuthorDate: 2021-02-10 22:28:43 +0000 Commit: Dimitry Andric <dim@FreeBSD.org> CommitDate: 2021-02-10 22:28:43 +0000 Fix incorrect hypotl(3) result with subnormal numbers This adjusts the factor used to scale the subnormal numbers, so it becomes the right value after adjusting its exponent. Thanks to Steve Kargl for finding the most elegant fix. Also enable the hypot tests, and add a test case for this bug. PR: 253313 MFC after: 1 week contrib/netbsd-tests/lib/libm/t_hypot.c | 20 ++++++++++++++++++++ lib/msun/src/e_hypotl.c | 2 +- lib/msun/tests/Makefile | 1 + 3 files changed, 22 insertions(+), 1 deletion(-)
(In reply to commit-hook from comment #3) FYI, from the lists: On Wed, Feb 10, 2021, at 4:30 PM, Dimitry Andric wrote: > + volatile long double a = 0x1.b2933cafa0bb7p-16383L; > + volatile long double b = 0x1.fffffffffffffp-16351L; > + volatile long double e = 0x1.fffffffffffffp-16351L; These are compile errors on most platforms. I believe these constants are specific to the x86-centric 80-bit long double format, whereas all other platforms use either 64-bit or 128-bit long doubles.
A commit in branch main references this bug: URL: https://cgit.FreeBSD.org/src/commit/?id=25120662284466ecef976df8f86e97bafdedf991 commit 25120662284466ecef976df8f86e97bafdedf991 Author: Dimitry Andric <dim@FreeBSD.org> AuthorDate: 2021-02-11 11:01:10 +0000 Commit: Dimitry Andric <dim@FreeBSD.org> CommitDate: 2021-02-11 11:01:10 +0000 Fix lib/msun/test builds on platforms without 80-bit long doubles After d3338f3355a612cf385632291f46c5777bba8d18, the lib/msun test case 'hypotl_near_underflow' would fail to compile on platforms where long doubles weren't 80 bit, like on x86. Disable this particular test on such platforms for now. PR: 253313 MFC after: 1 week X-MFC-With: d3338f3355a612cf385632291f46c5777bba8d18 contrib/netbsd-tests/lib/libm/t_hypot.c | 4 ++++ 1 file changed, 4 insertions(+)
A commit in branch stable/13 references this bug: URL: https://cgit.FreeBSD.org/src/commit/?id=97232b0dc03146c0e6466a4886edaaca4f95c3a7 commit 97232b0dc03146c0e6466a4886edaaca4f95c3a7 Author: Dimitry Andric <dim@FreeBSD.org> AuthorDate: 2021-02-10 22:28:43 +0000 Commit: Dimitry Andric <dim@FreeBSD.org> CommitDate: 2021-02-18 18:08:27 +0000 Fix incorrect hypotl(3) result with subnormal numbers This adjusts the factor used to scale the subnormal numbers, so it becomes the right value after adjusting its exponent. Thanks to Steve Kargl for finding the most elegant fix. Also enable the hypot tests, and add a test case for this bug. PR: 253313 (cherry picked from commit d3338f3355a612cf385632291f46c5777bba8d18) Fix lib/msun/test builds on platforms without 80-bit long doubles After d3338f3355a612cf385632291f46c5777bba8d18, the lib/msun test case 'hypotl_near_underflow' would fail to compile on platforms where long doubles weren't 80 bit, like on x86. Disable this particular test on such platforms for now. PR: 253313 (cherry picked from commit 25120662284466ecef976df8f86e97bafdedf991) contrib/netbsd-tests/lib/libm/t_hypot.c | 24 ++++++++++++++++++++++++ lib/msun/src/e_hypotl.c | 2 +- lib/msun/tests/Makefile | 1 + 3 files changed, 26 insertions(+), 1 deletion(-)
A commit in branch stable/12 references this bug: URL: https://cgit.FreeBSD.org/src/commit/?id=be49cd3ec6a53aa6b0b5c5635202dd92bcaacfc6 commit be49cd3ec6a53aa6b0b5c5635202dd92bcaacfc6 Author: Dimitry Andric <dim@FreeBSD.org> AuthorDate: 2021-02-10 22:28:43 +0000 Commit: Dimitry Andric <dim@FreeBSD.org> CommitDate: 2021-02-18 18:10:32 +0000 Fix incorrect hypotl(3) result with subnormal numbers This adjusts the factor used to scale the subnormal numbers, so it becomes the right value after adjusting its exponent. Thanks to Steve Kargl for finding the most elegant fix. Also enable the hypot tests, and add a test case for this bug. PR: 253313 (cherry picked from commit d3338f3355a612cf385632291f46c5777bba8d18) Fix lib/msun/test builds on platforms without 80-bit long doubles After d3338f3355a612cf385632291f46c5777bba8d18, the lib/msun test case 'hypotl_near_underflow' would fail to compile on platforms where long doubles weren't 80 bit, like on x86. Disable this particular test on such platforms for now. PR: 253313 (cherry picked from commit 25120662284466ecef976df8f86e97bafdedf991) contrib/netbsd-tests/lib/libm/t_hypot.c | 24 ++++++++++++++++++++++++ lib/msun/src/e_hypotl.c | 2 +- lib/msun/tests/Makefile | 1 + 3 files changed, 26 insertions(+), 1 deletion(-)
A commit in branch stable/11 references this bug: URL: https://cgit.FreeBSD.org/src/commit/?id=6287b5978c5f3959ad7f9494681cf94d1243120e commit 6287b5978c5f3959ad7f9494681cf94d1243120e Author: Dimitry Andric <dim@FreeBSD.org> AuthorDate: 2021-02-10 22:28:43 +0000 Commit: Dimitry Andric <dim@FreeBSD.org> CommitDate: 2021-02-18 18:11:19 +0000 Fix incorrect hypotl(3) result with subnormal numbers This adjusts the factor used to scale the subnormal numbers, so it becomes the right value after adjusting its exponent. Thanks to Steve Kargl for finding the most elegant fix. Also enable the hypot tests, and add a test case for this bug. PR: 253313 (cherry picked from commit d3338f3355a612cf385632291f46c5777bba8d18) Fix lib/msun/test builds on platforms without 80-bit long doubles After d3338f3355a612cf385632291f46c5777bba8d18, the lib/msun test case 'hypotl_near_underflow' would fail to compile on platforms where long doubles weren't 80 bit, like on x86. Disable this particular test on such platforms for now. PR: 253313 (cherry picked from commit 25120662284466ecef976df8f86e97bafdedf991) contrib/netbsd-tests/lib/libm/t_hypot.c | 24 ++++++++++++++++++++++++ lib/msun/src/e_hypotl.c | 2 +- lib/msun/tests/Makefile | 1 + 3 files changed, 26 insertions(+), 1 deletion(-)