scalbnl patch
Stephen L Moshier
steve@moshier.net
Thu Jun 13 06:32:00 GMT 2002
The long double scalbnl function seems to have several bugs. I don't know
if it actually gets used in glibc any more but here are a test program and
a patch.
--------------
extern long double scalbnl (long double, int);
main()
{
long double x, y, p2, z;
int i, j;
x = 4.0;
p2 = 8.0;
for (i = 0; i < (16384+64); i++)
{
p2 = 0.5L * p2;
y = scalbnl (x, -i);
if (y != p2)
printf ("2**-%d: p2 = %.4Le, y = %.4Le\n", i, p2, y);
z = scalbnl (p2, i);
if (z != 4.0)
printf ("? scalb(%.4Le, %d) = %.4Le\n", p2, i, z);
}
printf ("last: p2 = %.4Le, y = %.4Le\n", p2, y);
printf ("last: scalb(%.4Le, %d) = %.4Le\n", p2, i, z);
}
--------------
* sysdeps/ieee754/ldbl-96/s_scalbnl.c: Fix cases in which
argument or result is subnormal.
*** s_scalbnl.c 1999/07/14 00:15:06 1.1
--- s_scalbnl.c 2002/06/13 13:00:07
***************
*** 33,40 ****
#else
static long double
#endif
! two63 = 4.50359962737049600000e+15,
! twom63 = 1.08420217248550443400e-19,
huge = 1.0e+4900L,
tiny = 1.0e-4900L;
--- 33,40 ----
#else
static long double
#endif
! two64 = 1.8446744073709551616e19L,
! twom64 = 5.421010862427522170037e-20L,
huge = 1.0e+4900L,
tiny = 1.0e-4900L;
***************
*** 50,58 ****
k = es&0x7fff; /* extract exponent */
if (k==0) { /* 0 or subnormal x */
if ((lx|(hx&0x7fffffff))==0) return x; /* +-0 */
! x *= two63;
! GET_LDOUBLE_EXP(es,x);
! k = (hx&0x7fff) - 63;
}
if (k==0x7fff) return x+x; /* NaN or Inf */
k = k+n;
--- 50,58 ----
k = es&0x7fff; /* extract exponent */
if (k==0) { /* 0 or subnormal x */
if ((lx|(hx&0x7fffffff))==0) return x; /* +-0 */
! x *= two64;
! GET_LDOUBLE_EXP(hx,x);
! k = (hx&0x7fff) - 64;
}
if (k==0x7fff) return x+x; /* NaN or Inf */
k = k+n;
***************
*** 62,71 ****
return tiny*__copysignl(tiny,x);
if (k > 0) /* normal result */
{SET_LDOUBLE_EXP(x,(es&0x8000)|k); return x;}
! if (k <= -63)
return tiny*__copysignl(tiny,x); /*underflow*/
! k += 63; /* subnormal result */
SET_LDOUBLE_EXP(x,(es&0x8000)|k);
! return x*twom63;
}
weak_alias (__scalbnl, scalbnl)
--- 62,71 ----
return tiny*__copysignl(tiny,x);
if (k > 0) /* normal result */
{SET_LDOUBLE_EXP(x,(es&0x8000)|k); return x;}
! if (k <= -64)
return tiny*__copysignl(tiny,x); /*underflow*/
! k += 64; /* subnormal result */
SET_LDOUBLE_EXP(x,(es&0x8000)|k);
! return x*twom64;
}
weak_alias (__scalbnl, scalbnl)
More information about the Libc-alpha
mailing list