This is the mail archive of the libc-alpha@sources.redhat.com mailing list for the glibc project.


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]
Other format: [Raw text]

scalbnl patch



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)


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]