hypothl(x) mishandles subnormal numbers.

From: Steve Kargl <sgk_at_troutmask.apl.washington.edu>
Date: Sat, 6 Feb 2021 12:39:29 -0800
I've long forgotten by freebsd bugzilla password.  
So, if someone would like to submit a bug report,
here's a test program.

% cc -o z -O2 -fno-builtin z.c -lm
% ./z
Via scaling and sqrtl
x=2.853684e-4932 y=1.444012e-4922 a=1.444012e-4922
x=0x1.b2933cafa0bb8p-16383 y=0x1.fffffffffffffp-16351 a=0x1.000000006ca4c8p-16350

Via hypotl
x=2.853684e-4932 y=1.444012e-4922 z=nan
x=0x1.b2933cafa0bb8p-16383 y=0x1.fffffffffffffp-16351 z=nan

% cat z.c
#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;
}

-- 
Steve
Received on Sat Feb 06 2021 - 19:39:31 UTC

This archive was generated by hypermail 2.4.0 : Wed May 19 2021 - 11:41:27 UTC