[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...
- To: commits@xxxxxxxxxx
- Subject: [Commits] r22511 - in /fsf/trunk/libc: ./ sysdeps/ieee754/dbl-64/ sysdeps/powerpc/powerpc32/power4/fpu/ sysdeps/powerpc/powerpc64/powe...
- From: eglibc@xxxxxxxxxx
- Date: Tue, 26 Feb 2013 00:01:55 -0000
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