Bug 253313 - lib/msun: hypotl(3) mishandles subnormal numbers
Summary: lib/msun: hypotl(3) mishandles subnormal numbers
Status: Closed FIXED
Alias: None
Product: Base System
Classification: Unclassified
Component: bin (show other bugs)
Version: CURRENT
Hardware: amd64 Any
: --- Affects Some People
Assignee: freebsd-bugs (Nobody)
URL:
Keywords:
Depends on:
Blocks:
 
Reported: 2021-02-06 21:29 UTC by Dimitry Andric
Modified: 2021-02-18 18:15 UTC (History)
2 users (show)

See Also:


Attachments

Note You need to log in before you can comment on or make changes to this bug.
Description Dimitry Andric freebsd_committer freebsd_triage 2021-02-06 21:29:02 UTC
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.
Comment 1 Dimitry Andric freebsd_committer freebsd_triage 2021-02-06 21:49:37 UTC
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.
Comment 2 Dimitry Andric freebsd_committer freebsd_triage 2021-02-07 16:17:17 UTC
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.)
Comment 3 commit-hook freebsd_committer freebsd_triage 2021-02-10 22:31:04 UTC
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(-)
Comment 4 Mark Millard 2021-02-11 05:00:07 UTC
(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.
Comment 5 commit-hook freebsd_committer freebsd_triage 2021-02-11 11:02:03 UTC
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(+)
Comment 6 commit-hook freebsd_committer freebsd_triage 2021-02-18 18:13:50 UTC
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(-)
Comment 7 commit-hook freebsd_committer freebsd_triage 2021-02-18 18:13:51 UTC
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(-)
Comment 8 commit-hook freebsd_committer freebsd_triage 2021-02-18 18:13:52 UTC
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(-)
Comment 9 commit-hook freebsd_committer freebsd_triage 2021-02-18 18:13:53 UTC
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(-)
Comment 10 commit-hook freebsd_committer freebsd_triage 2021-02-18 18:14:56 UTC
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(-)
Comment 11 commit-hook freebsd_committer freebsd_triage 2021-02-18 18:14:57 UTC
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(-)