[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[Commits] r22511 - in /fsf/trunk/libc: ./ sysdeps/ieee754/dbl-64/ sysdeps/powerpc/powerpc32/power4/fpu/ sysdeps/powerpc/powerpc64/powe...



Author: eglibc
Date: Tue Feb 26 00:01:54 2013
New Revision: 22511

Log:
Import glibc-mainline for 2013-02-26

Modified:
    fsf/trunk/libc/ChangeLog
    fsf/trunk/libc/sysdeps/ieee754/dbl-64/mpa.c
    fsf/trunk/libc/sysdeps/ieee754/dbl-64/sincos32.c
    fsf/trunk/libc/sysdeps/ieee754/dbl-64/slowexp.c
    fsf/trunk/libc/sysdeps/ieee754/dbl-64/slowpow.c
    fsf/trunk/libc/sysdeps/powerpc/powerpc32/power4/fpu/mpa.c
    fsf/trunk/libc/sysdeps/powerpc/powerpc64/power4/fpu/mpa.c

Modified: fsf/trunk/libc/ChangeLog
==============================================================================
--- fsf/trunk/libc/ChangeLog (original)
+++ fsf/trunk/libc/ChangeLog Tue Feb 26 00:01:54 2013
@@ -1,3 +1,34 @@
+2013-02-25  Siddhesh Poyarekar  <siddhesh@xxxxxxxxxx>
+
+	* sysdeps/ieee754/dbl-64/sincos32.c (ss32): Remove commented
+	code.
+	(cc32): Likewise.
+
+	* sysdeps/ieee754/dbl-64/mpa.c (mcr): Use long instead of int.
+	(__acr): Likewise.
+	(__cpy): Likewise.
+	(norm): Likewise.
+	(denorm): Likewise.
+	(__dbl_mp): Likewise.
+	(add_magnitudes): Likewise.
+	(sub_magnitudes): Likewise.
+	(__mul): Likewise.
+	(__inv): Likewise.
+
+	* sysdeps/ieee754/dbl-64/slowexp.c: Reformat in GNU coding
+	style.
+
+	* sysdeps/ieee754/dbl-64/slowpow.c: Reformat in GNU coding
+	style.
+
+	* sysdeps/ieee754/dbl-64/slowexp.c (__slowexp): Remove commented
+	code.
+
+	* sysdeps/powerpc/powerpc32/power4/fpu/mpa.c (__mp_dbl): Sync
+	up changes with default code.
+	* sysdeps/powerpc/powerpc64/power4/fpu/mpa.c (__mp_dbl):
+	Likewise.
+
 2013-02-24  Allan McRae  <allan@xxxxxxxxxxxxx>
 
 	* manual/socket.texi (The Internet Namespace): Order menu items

Modified: fsf/trunk/libc/sysdeps/ieee754/dbl-64/mpa.c
==============================================================================
--- fsf/trunk/libc/sysdeps/ieee754/dbl-64/mpa.c (original)
+++ fsf/trunk/libc/sysdeps/ieee754/dbl-64/mpa.c Tue Feb 26 00:01:54 2013
@@ -59,8 +59,9 @@
 static int
 mcr (const mp_no *x, const mp_no *y, int p)
 {
-  int i;
-  for (i = 1; i <= p; i++)
+  long i;
+  long p2 = p;
+  for (i = 1; i <= p2; i++)
     {
       if (X[i] == Y[i])
 	continue;
@@ -76,7 +77,7 @@
 int
 __acr (const mp_no *x, const mp_no *y, int p)
 {
-  int i;
+  long i;
 
   if (X[0] == ZERO)
     {
@@ -107,8 +108,10 @@
 void
 __cpy (const mp_no *x, mp_no *y, int p)
 {
+  long i;
+
   EY = EX;
-  for (int i = 0; i <= p; i++)
+  for (i = 0; i <= p; i++)
     Y[i] = X[i];
 }
 #endif
@@ -120,7 +123,7 @@
 norm (const mp_no *x, double *y, int p)
 {
 #define R RADIXI
-  int i;
+  long i;
   double a, c, u, v, z[5];
   if (p < 5)
     {
@@ -194,7 +197,8 @@
 static void
 denorm (const mp_no *x, double *y, int p)
 {
-  int i, k;
+  long i, k;
+  long p2 = p;
   double c, u, z[5];
 
 #define R RADIXI
@@ -204,7 +208,7 @@
       return;
     }
 
-  if (p == 1)
+  if (p2 == 1)
     {
       if (EX == -42)
 	{
@@ -228,7 +232,7 @@
 	  k = 1;
 	}
     }
-  else if (p == 2)
+  else if (p2 == 2)
     {
       if (EX == -42)
 	{
@@ -281,7 +285,7 @@
 
   if (u == z[3])
     {
-      for (i = k + 1; i <= p; i++)
+      for (i = k + 1; i <= p2; i++)
 	{
 	  if (X[i] == ZERO)
 	    continue;
@@ -323,7 +327,8 @@
 SECTION
 __dbl_mp (double x, mp_no *y, int p)
 {
-  int i, n;
+  long i, n;
+  long p2 = p;
   double u;
 
   /* Sign.  */
@@ -347,7 +352,7 @@
     x *= RADIX;
 
   /* Digits.  */
-  n = MIN (p, 4);
+  n = MIN (p2, 4);
   for (i = 1; i <= n; i++)
     {
       u = (x + TWO52) - TWO52;
@@ -357,7 +362,7 @@
       x -= u;
       x *= RADIX;
     }
-  for (; i <= p; i++)
+  for (; i <= p2; i++)
     Y[i] = ZERO;
 }
 
@@ -369,14 +374,15 @@
 SECTION
 add_magnitudes (const mp_no *x, const mp_no *y, mp_no *z, int p)
 {
-  int i, j, k;
+  long i, j, k;
+  long p2 = p;
   double zk;
 
   EZ = EX;
 
-  i = p;
-  j = p + EY - EX;
-  k = p + 1;
+  i = p2;
+  j = p2 + EY - EX;
+  k = p2 + 1;
 
   if (__glibc_unlikely (j < 1))
     {
@@ -418,7 +424,7 @@
 
   if (zk == ZERO)
     {
-      for (i = 1; i <= p; i++)
+      for (i = 1; i <= p2; i++)
 	Z[i] = Z[i + 1];
     }
   else
@@ -436,13 +442,14 @@
 SECTION
 sub_magnitudes (const mp_no *x, const mp_no *y, mp_no *z, int p)
 {
-  int i, j, k;
+  long i, j, k;
+  long p2 = p;
   double zk;
 
   EZ = EX;
-  i = p;
-  j = p + EY - EX;
-  k = p;
+  i = p2;
+  j = p2 + EY - EX;
+  k = p2;
 
   /* Y is too small compared to X, copy X over to the result.  */
   if (__glibc_unlikely (j < 1))
@@ -453,7 +460,7 @@
 
   /* The relevant least significant digit in Y is non-zero, so we factor it in
      to enhance accuracy.  */
-  if (j < p && Y[j + 1] > ZERO)
+  if (j < p2 && Y[j + 1] > ZERO)
     {
       Z[k + 1] = RADIX - Y[j + 1];
       zk = MONE;
@@ -496,9 +503,9 @@
   /* Normalize.  */
   for (i = 1; Z[i] == ZERO; i++);
   EZ = EZ - i + 1;
-  for (k = 1; i <= p + 1;)
+  for (k = 1; i <= p2 + 1;)
     Z[k++] = Z[i++];
-  for (; k <= p;)
+  for (; k <= p2;)
     Z[k++] = ZERO;
 }
 
@@ -610,7 +617,8 @@
 SECTION
 __mul (const mp_no *x, const mp_no *y, mp_no *z, int p)
 {
-  int i, j, k, ip, ip2;
+  long i, j, k, ip, ip2;
+  long p2 = p;
   double u, zk;
   const mp_no *a;
 
@@ -623,7 +631,7 @@
 
   /* We need not iterate through all X's and Y's since it's pointless to
      multiply zeroes.  Here, both are zero...  */
-  for (ip2 = p; ip2 > 0; ip2--)
+  for (ip2 = p2; ip2 > 0; ip2--)
     if (X[ip2] != ZERO || Y[ip2] != ZERO)
       break;
 
@@ -660,16 +668,16 @@
      'internal precision' of the input numbers, i.e. digits after ip and ip2
      are all ZERO.  */
 
-  k = (__glibc_unlikely (p < 3)) ? p + p : p + 3;
+  k = (__glibc_unlikely (p2 < 3)) ? p2 + p2 : p2 + 3;
 
   while (k > ip + ip2 + 1)
     Z[k--] = ZERO;
 
   zk = Z[k] = ZERO;
 
-  while (k > p)
-    {
-      for (i = k - p, j = p; i < p + 1; i++, j--)
+  while (k > p2)
+    {
+      for (i = k - p2, j = p2; i < p2 + 1; i++, j--)
 	zk += X[i] * Y[j];
 
       u = (zk + CUTTER) - CUTTER;
@@ -701,7 +709,7 @@
   /* Is there a carry beyond the most significant digit?  */
   if (__glibc_unlikely (Z[1] == ZERO))
     {
-      for (i = 1; i <= p; i++)
+      for (i = 1; i <= p2; i++)
 	Z[i] = Z[i + 1];
       e--;
     }
@@ -821,7 +829,7 @@
 SECTION
 __inv (const mp_no *x, mp_no *y, int p)
 {
-  int i;
+  long i;
   double t;
   mp_no z, w;
   static const int np1[] =

Modified: fsf/trunk/libc/sysdeps/ieee754/dbl-64/sincos32.c
==============================================================================
--- fsf/trunk/libc/sysdeps/ieee754/dbl-64/sincos32.c (original)
+++ fsf/trunk/libc/sysdeps/ieee754/dbl-64/sincos32.c Tue Feb 26 00:01:54 2013
@@ -57,14 +57,7 @@
 ss32(mp_no *x, mp_no *y, int p) {
   int i;
   double a;
-#if 0
-  double b;
-  static const mp_no mpone = {1,{1.0,1.0}};
-#endif
   mp_no mpt1,x2,gor,sum ,mpk={1,{1.0}};
-#if 0
-  mp_no mpt2;
-#endif
   for (i=1;i<=p;i++) mpk.d[i]=0;
 
   __sqr(x,&x2,p);
@@ -89,14 +82,7 @@
 cc32(mp_no *x, mp_no *y, int p) {
   int i;
   double a;
-#if 0
-  double b;
-  static const mp_no mpone = {1,{1.0,1.0}};
-#endif
   mp_no mpt1,x2,gor,sum ,mpk={1,{1.0}};
-#if 0
-  mp_no mpt2;
-#endif
   for (i=1;i<=p;i++) mpk.d[i]=0;
 
   __sqr(x,&x2,p);

Modified: fsf/trunk/libc/sysdeps/ieee754/dbl-64/slowexp.c
==============================================================================
--- fsf/trunk/libc/sysdeps/ieee754/dbl-64/slowexp.c (original)
+++ fsf/trunk/libc/sysdeps/ieee754/dbl-64/slowexp.c Tue Feb 26 00:01:54 2013
@@ -34,38 +34,36 @@
 # define SECTION
 #endif
 
-void __mpexp(mp_no *x, mp_no *y, int p);
+void __mpexp (mp_no *x, mp_no *y, int p);
 
 /*Converting from double precision to Multi-precision and calculating  e^x */
 double
 SECTION
-__slowexp(double x) {
-  double w,z,res,eps=3.0e-26;
-#if 0
-  double y;
-#endif
+__slowexp (double x)
+{
+  double w, z, res, eps = 3.0e-26;
   int p;
-#if 0
-  int orig,i;
-#endif
-  mp_no mpx, mpy, mpz,mpw,mpeps,mpcor;
+  mp_no mpx, mpy, mpz, mpw, mpeps, mpcor;
 
-  p=6;
-  __dbl_mp(x,&mpx,p); /* Convert a double precision number  x               */
-		    /* into a multiple precision number mpx with prec. p. */
-  __mpexp(&mpx, &mpy, p); /* Multi-Precision exponential function */
-  __dbl_mp(eps,&mpeps,p);
-  __mul(&mpeps,&mpy,&mpcor,p);
-  __add(&mpy,&mpcor,&mpw,p);
-  __sub(&mpy,&mpcor,&mpz,p);
-  __mp_dbl(&mpw, &w, p);
-  __mp_dbl(&mpz, &z, p);
-  if (w == z) return w;
-  else  {                   /* if calculating is not exactly   */
-    p = 32;
-    __dbl_mp(x,&mpx,p);
-    __mpexp(&mpx, &mpy, p);
-    __mp_dbl(&mpy, &res, p);
-    return res;
-  }
+  /* Use the multiple precision __MPEXP function to compute the exponential
+     First at 144 bits and if it is not accurate enough, at 768 bits.  */
+  p = 6;
+  __dbl_mp (x, &mpx, p);
+  __mpexp (&mpx, &mpy, p);
+  __dbl_mp (eps, &mpeps, p);
+  __mul (&mpeps, &mpy, &mpcor, p);
+  __add (&mpy, &mpcor, &mpw, p);
+  __sub (&mpy, &mpcor, &mpz, p);
+  __mp_dbl (&mpw, &w, p);
+  __mp_dbl (&mpz, &z, p);
+  if (w == z)
+    return w;
+  else
+    {
+      p = 32;
+      __dbl_mp (x, &mpx, p);
+      __mpexp (&mpx, &mpy, p);
+      __mp_dbl (&mpy, &res, p);
+      return res;
+    }
 }

Modified: fsf/trunk/libc/sysdeps/ieee754/dbl-64/slowpow.c
==============================================================================
--- fsf/trunk/libc/sysdeps/ieee754/dbl-64/slowpow.c (original)
+++ fsf/trunk/libc/sysdeps/ieee754/dbl-64/slowpow.c Tue Feb 26 00:01:54 2013
@@ -38,42 +38,59 @@
 # define SECTION
 #endif
 
-void __mpexp(mp_no *x, mp_no *y, int p);
-void __mplog(mp_no *x, mp_no *y, int p);
-double ulog(double);
-double __halfulp(double x,double y);
+void __mpexp (mp_no *x, mp_no *y, int p);
+void __mplog (mp_no *x, mp_no *y, int p);
+double ulog (double);
+double __halfulp (double x, double y);
 
 double
 SECTION
-__slowpow(double x, double y, double z) {
-  double res,res1;
-  mp_no mpx, mpy, mpz,mpw,mpp,mpr,mpr1;
-  static const mp_no eps = {-3,{1.0,4.0}};
+__slowpow (double x, double y, double z)
+{
+  double res, res1;
+  mp_no mpx, mpy, mpz, mpw, mpp, mpr, mpr1;
+  static const mp_no eps = {-3, {1.0, 4.0}};
   int p;
 
-  res = __halfulp(x,y);        /* halfulp() returns -10 or x^y             */
-  if (res >= 0) return res;  /* if result was really computed by halfulp */
-		  /*  else, if result was not really computed by halfulp */
-  p = 10;         /*  p=precision   */
-  __dbl_mp(x,&mpx,p);
-  __dbl_mp(y,&mpy,p);
-  __dbl_mp(z,&mpz,p);
-  __mplog(&mpx, &mpz, p);     /* log(x) = z   */
-  __mul(&mpy,&mpz,&mpw,p);    /*  y * z =w    */
-  __mpexp(&mpw, &mpp, p);     /*  e^w =pp     */
-  __add(&mpp,&eps,&mpr,p);    /*  pp+eps =r   */
-  __mp_dbl(&mpr, &res, p);
-  __sub(&mpp,&eps,&mpr1,p);   /*  pp -eps =r1 */
-  __mp_dbl(&mpr1, &res1, p);  /*  converting into double precision */
-  if (res == res1) return res;
+  /* __HALFULP returns -10 or X^Y.  */
+  res = __halfulp (x, y);
 
-  p = 32;     /* if we get here result wasn't calculated exactly, continue */
-  __dbl_mp(x,&mpx,p);                          /* for more exact calculation */
-  __dbl_mp(y,&mpy,p);
-  __dbl_mp(z,&mpz,p);
-  __mplog(&mpx, &mpz, p);   /* log(c)=z  */
-  __mul(&mpy,&mpz,&mpw,p);  /* y*z =w    */
-  __mpexp(&mpw, &mpp, p);   /* e^w=pp    */
-  __mp_dbl(&mpp, &res, p);  /* converting into double precision */
+  /* Return if the result was computed by __HALFULP.  */
+  if (res >= 0)
+    return res;
+
+  /* Or else, calculate using multiple precision.  P = 10 implies accuracy of
+     240 bits accuracy, since MP_NO has a radix of 2^24.  */
+  p = 10;
+  __dbl_mp (x, &mpx, p);
+  __dbl_mp (y, &mpy, p);
+  __dbl_mp (z, &mpz, p);
+
+  /* z = x ^ y
+     log (z) = y * log (x)
+     z = exp (y * log (x))  */
+  __mplog (&mpx, &mpz, p);
+  __mul (&mpy, &mpz, &mpw, p);
+  __mpexp (&mpw, &mpp, p);
+
+  /* Add and subtract EPS to ensure that the result remains unchanged, i.e. we
+     have last bit accuracy.  */
+  __add (&mpp, &eps, &mpr, p);
+  __mp_dbl (&mpr, &res, p);
+  __sub (&mpp, &eps, &mpr1, p);
+  __mp_dbl (&mpr1, &res1, p);
+  if (res == res1)
+    return res;
+
+  /* If we don't, then we repeat using a higher precision.  768 bits of
+     precision ought to be enough for anybody.  */
+  p = 32;
+  __dbl_mp (x, &mpx, p);
+  __dbl_mp (y, &mpy, p);
+  __dbl_mp (z, &mpz, p);
+  __mplog (&mpx, &mpz, p);
+  __mul (&mpy, &mpz, &mpw, p);
+  __mpexp (&mpw, &mpp, p);
+  __mp_dbl (&mpp, &res, p);
   return res;
 }

Modified: fsf/trunk/libc/sysdeps/powerpc/powerpc32/power4/fpu/mpa.c
==============================================================================
--- fsf/trunk/libc/sysdeps/powerpc/powerpc32/power4/fpu/mpa.c (original)
+++ fsf/trunk/libc/sysdeps/powerpc/powerpc32/power4/fpu/mpa.c Tue Feb 26 00:01:54 2013
@@ -304,9 +304,7 @@
       return;
     }
 
-  if (EX > -42)
-    norm (x, y, p);
-  else if (EX == -42 && X[1] >= TWO10)
+  if (__glibc_likely (EX > -42 || (EX == -42 && X[1] >= TWO10)))
     norm (x, y, p);
   else
     denorm (x, y, p);

Modified: fsf/trunk/libc/sysdeps/powerpc/powerpc64/power4/fpu/mpa.c
==============================================================================
--- fsf/trunk/libc/sysdeps/powerpc/powerpc64/power4/fpu/mpa.c (original)
+++ fsf/trunk/libc/sysdeps/powerpc/powerpc64/power4/fpu/mpa.c Tue Feb 26 00:01:54 2013
@@ -304,9 +304,7 @@
       return;
     }
 
-  if (EX > -42)
-    norm (x, y, p);
-  else if (EX == -42 && X[1] >= TWO10)
+  if (__glibc_likely (EX > -42 || (EX == -42 && X[1] >= TWO10)))
     norm (x, y, p);
   else
     denorm (x, y, p);

_______________________________________________
Commits mailing list
Commits@xxxxxxxxxx
http://eglibc.org/cgi-bin/mailman/listinfo/commits