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

[bugzilla-noreply at freebsd.org]

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

            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 at FreeBSD.org
          Reporter: dim at FreeBSD.org

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

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