Skip site navigation (1)Skip section navigation (2) 
Date:      Sat, 06 Feb 2021 21:29:02 +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@https.bugs.freebsd.org/bugzilla/>
next in thread | raw e-mail | index | archive | help 
https://bugs.freebsd.org/bugzilla/show_bug.cgi?id=3D253313

            Bug ID: 253313
           Summary: lib/msun: hypotl(3) mishandles subnormal numbers
           Product: Base System
           Version: CURRENT
          Hardware: amd64
                OS: Any
            Status: New
          Severity: Affects Some People
          Priority: ---
         Component: bin
          Assignee: bugs@FreeBSD.org
          Reporter: dim@FreeBSD.org

>From https://github.com/JuliaMath/openlibm/issues/224, via Steve Kargl:

Test program:
=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=
=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=
=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D
#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 =3D 0x3.6526795f4176ep-16384L;
  y =3D 0x3.ffffffffffffep-16352L;
  z =3D hypotl(x, y);

  if (x > y) {
     a =3D y / x;
     a =3D x * sqrtl(1 + a);
  } else {
     a =3D x / y;
     a =3D y * sqrtl(a + 1);
  }

  printf("Via scaling and sqrtl\n");
  printf("x=3D%Le y=3D%Le a=3D%Le\n", x, y, a);
  printf("x=3D%La y=3D%La a=3D%La\n", x, y, a);

  printf("\nVia hypotl\n");
  printf("x=3D%Le y=3D%Le z=3D%Le\n", x, y, z);
  printf("x=3D%La y=3D%La z=3D%La\n", x, y, z);

  return 0;
}
=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=
=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=
=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D

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=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=3Dnan
x=3D0x1.b2933cafa0bb7p-16383 y=3D0x1.fffffffffffffp-16351 z=3Dnan

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=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=3Dnan
x=3D0x1.b2933cafa0bb7p-16383 y=3D0x1.fffffffffffffp-16351 z=3Dnan

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=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

clang 11.0 apparently doesn't have a built-in hypotl, so the output isn't
different.

--=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>