[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Commits] r22724 - in /fsf/trunk/libc: ./ ports/ ports/sysdeps/arm/ sysdeps/ieee754/dbl-64/ sysdeps/powerpc/fpu/ sysdeps/powerpc/power...
- To: commits@xxxxxxxxxx
- Subject: [Commits] r22724 - in /fsf/trunk/libc: ./ ports/ ports/sysdeps/arm/ sysdeps/ieee754/dbl-64/ sysdeps/powerpc/fpu/ sysdeps/powerpc/power...
- From: eglibc@xxxxxxxxxx
- Date: Wed, 27 Mar 2013 00:01:58 -0000
Author: eglibc
Date: Wed Mar 27 00:01:57 2013
New Revision: 22724
Log:
Import glibc-mainline for 2013-03-27
Added:
fsf/trunk/libc/sysdeps/ieee754/dbl-64/mpa-arch.h
fsf/trunk/libc/sysdeps/powerpc/power4/fpu/mpa-arch.h
Modified:
fsf/trunk/libc/ChangeLog
fsf/trunk/libc/ports/ChangeLog.arm
fsf/trunk/libc/ports/sysdeps/arm/preconfigure
fsf/trunk/libc/ports/sysdeps/arm/preconfigure.in
fsf/trunk/libc/sysdeps/ieee754/dbl-64/branred.c
fsf/trunk/libc/sysdeps/ieee754/dbl-64/dosincos.c
fsf/trunk/libc/sysdeps/ieee754/dbl-64/e_asin.c
fsf/trunk/libc/sysdeps/ieee754/dbl-64/e_atan2.c
fsf/trunk/libc/sysdeps/ieee754/dbl-64/e_exp.c
fsf/trunk/libc/sysdeps/ieee754/dbl-64/e_log.c
fsf/trunk/libc/sysdeps/ieee754/dbl-64/e_pow.c
fsf/trunk/libc/sysdeps/ieee754/dbl-64/e_remainder.c
fsf/trunk/libc/sysdeps/ieee754/dbl-64/mpa.c
fsf/trunk/libc/sysdeps/ieee754/dbl-64/mpa.h
fsf/trunk/libc/sysdeps/ieee754/dbl-64/s_atan.c
fsf/trunk/libc/sysdeps/ieee754/dbl-64/s_sin.c
fsf/trunk/libc/sysdeps/ieee754/dbl-64/s_tan.c
fsf/trunk/libc/sysdeps/powerpc/fpu/s_llround.c
Modified: fsf/trunk/libc/ChangeLog
==============================================================================
--- fsf/trunk/libc/ChangeLog (original)
+++ fsf/trunk/libc/ChangeLog Wed Mar 27 00:01:57 2013
@@ -1,3 +1,62 @@
+2013-03-26 Siddhesh Poyarekar <siddhesh@xxxxxxxxxx>
+
+ * sysdeps/ieee754/dbl-64/mpa.c (__acr): Use integral
+ constants.
+ (norm): Likewise.
+ (denorm): Likewise.
+ (__dbl_mp): Likewise.
+ (add_magnitudes): Likewise.
+ (sub_magnitudes): Likewise.
+ (__add): Likewise.
+ (__sub): Likewise.
+ (__mul): Likewise.
+ (__sqr): Likewise.
+ (__inv): Likewise.
+ (__dvd): Likewise.
+
+ * sysdeps/ieee754/dbl-64/branred.c (__branred): Remove
+ commented code.
+ * sysdeps/ieee754/dbl-64/dosincos.c (__dubsin): Likewise.
+ (__dubcos): Likewise.
+ * sysdeps/ieee754/dbl-64/e_asin.c (__ieee754_asin): Likewise.
+ (__ieee754_acos): Likewise.
+ * sysdeps/ieee754/dbl-64/e_atan2.c (__ieee754_atan2): Likewise.
+ * sysdeps/ieee754/dbl-64/e_exp.c (__ieee754_exp): Likewise.
+ (__exp1): Likewise.
+ * sysdeps/ieee754/dbl-64/e_log.c (__ieee754_log): Likewise.
+ * sysdeps/ieee754/dbl-64/e_pow.c (__ieee754_pow): Likewise.
+ (log1): Likewise.
+ (my_log2): Likewise.
+ (checkint): Likewise.
+ * sysdeps/ieee754/dbl-64/e_remainder.c
+ (__ieee754_remainder): Likewise.
+ * sysdeps/ieee754/dbl-64/s_atan.c (atan): Likewise.
+ * sysdeps/ieee754/dbl-64/s_sin.c (__sin): Likewise.
+ (bsloww): Likewise.
+ * sysdeps/ieee754/dbl-64/s_tan.c (tan): Likewise.
+
+ * sysdeps/ieee754/dbl-64/mpa-arch.h: New file.
+ * sysdeps/ieee754/dbl-64/mpa.c (norm): Use MANTISSA_T and
+ MANTISSA_STORE_T to store computations on mantissa. Use
+ macros for rounding and division.
+ (denorm): Likewise.
+ (__dbl_mp): Likewise.
+ (add_magnitudes): Likewise.
+ (sub_magnitudes): Likewise.
+ (__mul): Likewise.
+ (__sqr): Likewise.
+ * sysdeps/ieee754/dbl-64/mpa.h: Include mpa-arch.h. Define
+ powers of two in terms of TWOPOW macro.
+ (mp_no): Make type of mantissa as MANTISSA_T.
+ [!RADIXI]: Define RADIXI.
+ [!TWO52]: Define TWO52.
+ * sysdeps/powerpc/power4/fpu/mpa-arch.h: New file.
+
+2013-03-25 Adhemerval Zanella <azanella@xxxxxxxxxxxxxxxxxx>
+
+ * sysdeps/powerpc/fpu/s_llround.c: Fix libm ABI issue with missing
+ llroundl symbol when building for PPC32.
+
2013-03-24 Mark H Weaver <mhw@xxxxxxxxxx>
* manual/arith.texi (Normalization Functions): Fix prototypes for
Modified: fsf/trunk/libc/ports/ChangeLog.arm
==============================================================================
--- fsf/trunk/libc/ports/ChangeLog.arm (original)
+++ fsf/trunk/libc/ports/ChangeLog.arm Wed Mar 27 00:01:57 2013
@@ -1,3 +1,8 @@
+2013-03-26 Mans Rullgard <mans@xxxxxxxxx>
+
+ * sysdeps/arm/preconfigure.in: Use "test" instead of [ ].
+ * sysdeps/arm/preconfigure: Regenerated.
+
2013-03-20 Joseph Myers <joseph@xxxxxxxxxxxxxxxx>
* sysdeps/arm/configure.in (default-abi): Set using
Modified: fsf/trunk/libc/ports/sysdeps/arm/preconfigure
==============================================================================
--- fsf/trunk/libc/ports/sysdeps/arm/preconfigure (original)
+++ fsf/trunk/libc/ports/sysdeps/arm/preconfigure Wed Mar 27 00:01:57 2013
@@ -10,7 +10,7 @@
# avoid this, add -fno-unwind-tables here and remove it in
# sysdeps/unix/sysv/linux/arm/configure.in after those tests have
# been run.
- if "${CFLAGS+set}" != "set" ; then
+ if test "${CFLAGS+set}" != "set"; then
CFLAGS="-g -O2"
fi
CFLAGS="$CFLAGS -fno-unwind-tables"
Modified: fsf/trunk/libc/ports/sysdeps/arm/preconfigure.in
==============================================================================
--- fsf/trunk/libc/ports/sysdeps/arm/preconfigure.in (original)
+++ fsf/trunk/libc/ports/sysdeps/arm/preconfigure.in Wed Mar 27 00:01:57 2013
@@ -10,7 +10,7 @@
# avoid this, add -fno-unwind-tables here and remove it in
# sysdeps/unix/sysv/linux/arm/configure.in after those tests have
# been run.
- if [ "${CFLAGS+set}" != "set" ]; then
+ if test "${CFLAGS+set}" != "set"; then
CFLAGS="-g -O2"
fi
CFLAGS="$CFLAGS -fno-unwind-tables"
Modified: fsf/trunk/libc/sysdeps/ieee754/dbl-64/branred.c
==============================================================================
--- fsf/trunk/libc/sysdeps/ieee754/dbl-64/branred.c (original)
+++ fsf/trunk/libc/sysdeps/ieee754/dbl-64/branred.c Wed Mar 27 00:01:57 2013
@@ -53,13 +53,7 @@
__branred(double x, double *a, double *aa)
{
int i,k;
-#if 0
- int n;
-#endif
mynumber u,gor;
-#if 0
- mynumber v;
-#endif
double r[6],s,t,sum,b,bb,sum1,sum2,b1,bb1,b2,bb2,x1,x2,t1,t2;
x*=tm600.x;
Modified: fsf/trunk/libc/sysdeps/ieee754/dbl-64/dosincos.c
==============================================================================
--- fsf/trunk/libc/sysdeps/ieee754/dbl-64/dosincos.c (original)
+++ fsf/trunk/libc/sysdeps/ieee754/dbl-64/dosincos.c Wed Mar 27 00:01:57 2013
@@ -63,9 +63,6 @@
#ifndef DLA_FMS
double p,hx,tx,hy,ty,q;
#endif
-#if 0
- double xx,y,yy,z,zz;
-#endif
mynumber u;
int4 k;
@@ -118,9 +115,6 @@
sn,ssn,cs,ccs,ds,dss,dc,dcc;
#ifndef DLA_FMS
double p,hx,tx,hy,ty,q;
-#endif
-#if 0
- double xx,y,yy,z,zz;
#endif
mynumber u;
int4 k;
Modified: fsf/trunk/libc/sysdeps/ieee754/dbl-64/e_asin.c
==============================================================================
--- fsf/trunk/libc/sysdeps/ieee754/dbl-64/e_asin.c (original)
+++ fsf/trunk/libc/sysdeps/ieee754/dbl-64/e_asin.c Wed Mar 27 00:01:57 2013
@@ -62,9 +62,6 @@
double x1,x2,xx,s1,s2,res1,p,t,res,r,cor,cc,y,c,z,w[2];
mynumber u,v;
int4 k,m,n;
-#if 0
- int4 nn;
-#endif
u.x = x;
m = u.i[HIGH_HALF];
@@ -344,14 +341,8 @@
__ieee754_acos(double x)
{
double x1,x2,xx,s1,s2,res1,p,t,res,r,cor,cc,y,c,z,w[2],eps;
-#if 0
- double fc;
-#endif
mynumber u,v;
int4 k,m,n;
-#if 0
- int4 nn;
-#endif
u.x = x;
m = u.i[HIGH_HALF];
k = 0x7fffffff&m;
Modified: fsf/trunk/libc/sysdeps/ieee754/dbl-64/e_atan2.c
==============================================================================
--- fsf/trunk/libc/sysdeps/ieee754/dbl-64/e_atan2.c (original)
+++ fsf/trunk/libc/sysdeps/ieee754/dbl-64/e_atan2.c Wed Mar 27 00:01:57 2013
@@ -67,22 +67,13 @@
__ieee754_atan2(double y,double x) {
int i,de,ux,dx,uy,dy;
-#if 0
- int p;
-#endif
static const int pr[MM]={6,8,10,20,32};
double ax,ay,u,du,u9,ua,v,vv,dv,t1,t2,t3,t7,t8,
z,zz,cor,s1,ss1,s2,ss2;
#ifndef DLA_FMS
double t4,t5,t6;
#endif
-#if 0
- double z1,z2;
-#endif
number num;
-#if 0
- mp_no mperr,mpt1,mpx,mpy,mpz,mpz1,mpz2;
-#endif
static const int ep= 59768832, /* 57*16**5 */
em=-59768832; /* -57*16**5 */
Modified: fsf/trunk/libc/sysdeps/ieee754/dbl-64/e_exp.c
==============================================================================
--- fsf/trunk/libc/sysdeps/ieee754/dbl-64/e_exp.c (original)
+++ fsf/trunk/libc/sysdeps/ieee754/dbl-64/e_exp.c Wed Mar 27 00:01:57 2013
@@ -55,9 +55,6 @@
__ieee754_exp(double x) {
double bexp, t, eps, del, base, y, al, bet, res, rem, cor;
mynumber junk1, junk2, binexp = {{0,0}};
-#if 0
- int4 k;
-#endif
int4 i,j,m,n,ex;
double retval;
@@ -174,9 +171,6 @@
__exp1(double x, double xx, double error) {
double bexp, t, eps, del, base, y, al, bet, res, rem, cor;
mynumber junk1, junk2, binexp = {{0,0}};
-#if 0
- int4 k;
-#endif
int4 i,j,m,n,ex;
junk1.x = x;
Modified: fsf/trunk/libc/sysdeps/ieee754/dbl-64/e_log.c
==============================================================================
--- fsf/trunk/libc/sysdeps/ieee754/dbl-64/e_log.c (original)
+++ fsf/trunk/libc/sysdeps/ieee754/dbl-64/e_log.c Wed Mar 27 00:01:57 2013
@@ -56,9 +56,6 @@
#define M 4
static const int pr[M]={8,10,18,32};
int i,j,n,ux,dx,p;
-#if 0
- int k;
-#endif
double dbl_n,u,p0,q,r0,w,nln2a,luai,lubi,lvaj,lvbj,
sij,ssij,ttij,A,B,B0,y,y1,y2,polI,polII,sa,sb,
t1,t2,t7,t8,t,ra,rb,ww,
Modified: fsf/trunk/libc/sysdeps/ieee754/dbl-64/e_pow.c
==============================================================================
--- fsf/trunk/libc/sysdeps/ieee754/dbl-64/e_pow.c (original)
+++ fsf/trunk/libc/sysdeps/ieee754/dbl-64/e_pow.c Wed Mar 27 00:01:57 2013
@@ -64,9 +64,6 @@
SECTION
__ieee754_pow(double x, double y) {
double z,a,aa,error, t,a1,a2,y1,y2;
-#if 0
- double gor=1.0;
-#endif
mynumber u,v;
int k;
int4 qx,qy;
@@ -206,13 +203,7 @@
SECTION
log1(double x, double *delta, double *error) {
int i,j,m;
-#if 0
- int n;
-#endif
double uu,vv,eps,nx,e,e1,e2,t,t1,t2,res,add=0;
-#if 0
- double cor;
-#endif
mynumber u,v;
#ifdef BIG_ENDI
mynumber
@@ -300,13 +291,7 @@
SECTION
my_log2(double x, double *delta, double *error) {
int i,j,m;
-#if 0
- int n;
-#endif
double uu,vv,eps,nx,e,e1,e2,t,t1,t2,res,add=0;
-#if 0
- double cor;
-#endif
double ou1,ou2,lu1,lu2,ov,lv1,lv2,a,a1,a2;
double y,yy,z,zz,j1,j2,j7,j8;
#ifndef DLA_FMS
@@ -397,9 +382,6 @@
checkint(double x) {
union {int4 i[2]; double x;} u;
int k,m,n;
-#if 0
- int l;
-#endif
u.x = x;
m = u.i[HIGH_HALF]&0x7fffffff; /* no sign */
if (m >= 0x7ff00000) return 0; /* x is +/-inf or NaN */
Modified: fsf/trunk/libc/sysdeps/ieee754/dbl-64/e_remainder.c
==============================================================================
--- fsf/trunk/libc/sysdeps/ieee754/dbl-64/e_remainder.c (original)
+++ fsf/trunk/libc/sysdeps/ieee754/dbl-64/e_remainder.c Wed Mar 27 00:01:57 2013
@@ -42,13 +42,7 @@
double __ieee754_remainder(double x, double y)
{
double z,d,xx;
-#if 0
- double yy;
-#endif
int4 kx,ky,n,nn,n1,m1,l;
-#if 0
- int4 m;
-#endif
mynumber u,t,w={{0,0}},v={{0,0}},ww={{0,0}},r;
u.x=x;
t.x=y;
Added: fsf/trunk/libc/sysdeps/ieee754/dbl-64/mpa-arch.h
==============================================================================
--- fsf/trunk/libc/sysdeps/ieee754/dbl-64/mpa-arch.h (added)
+++ fsf/trunk/libc/sysdeps/ieee754/dbl-64/mpa-arch.h Wed Mar 27 00:01:57 2013
@@ -1,0 +1,47 @@
+/* Overridable constants and operations.
+ Copyright (C) 2013 Free Software Foundation, Inc.
+
+ This program is free software; you can redistribute it and/or modify
+ it under the terms of the GNU Lesser General Public License as published by
+ the Free Software Foundation; either version 2.1 of the License, or
+ (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU Lesser General Public License for more details.
+
+ You should have received a copy of the GNU Lesser General Public License
+ along with this program; if not, see <http://www.gnu.org/licenses/>. */
+
+#include <stdint.h>
+
+typedef long mantissa_t;
+typedef int64_t mantissa_store_t;
+
+#define TWOPOW(i) (1L << i)
+
+#define RADIX_EXP 24
+#define RADIX TWOPOW (RADIX_EXP) /* 2^24 */
+
+/* Divide D by RADIX and put the remainder in R. D must be a non-negative
+ integral value. */
+#define DIV_RADIX(d, r) \
+ ({ \
+ r = d & (RADIX - 1); \
+ d >>= RADIX_EXP; \
+ })
+
+/* Put the integer component of a double X in R and retain the fraction in
+ X. This is used in extracting mantissa digits for MP_NO by using the
+ integer portion of the current value of the number as the current mantissa
+ digit and then scaling by RADIX to get the next mantissa digit in the same
+ manner. */
+#define INTEGER_OF(x, i) \
+ ({ \
+ i = (mantissa_t) x; \
+ x -= i; \
+ })
+
+/* Align IN down to F. The code assumes that F is a power of two. */
+#define ALIGN_DOWN_TO(in, f) ((in) & -(f))
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 Wed Mar 27 00:01:57 2013
@@ -80,14 +80,14 @@
{
long i;
- if (X[0] == ZERO)
- {
- if (Y[0] == ZERO)
+ if (X[0] == 0)
+ {
+ if (Y[0] == 0)
i = 0;
else
i = -1;
}
- else if (Y[0] == ZERO)
+ else if (Y[0] == 0)
i = 1;
else
{
@@ -125,7 +125,8 @@
{
#define R RADIXI
long i;
- double a, c, u, v, z[5];
+ double c;
+ mantissa_t a, u, v, z[5];
if (p < 5)
{
if (p == 1)
@@ -139,44 +140,41 @@
}
else
{
- for (a = ONE, z[1] = X[1]; z[1] < TWO23;)
- {
- a *= TWO;
- z[1] *= TWO;
+ for (a = 1, z[1] = X[1]; z[1] < TWO23;)
+ {
+ a *= 2;
+ z[1] *= 2;
}
for (i = 2; i < 5; i++)
{
- z[i] = X[i] * a;
- u = (z[i] + CUTTER) - CUTTER;
- if (u > z[i])
- u -= RADIX;
- z[i] -= u;
- z[i - 1] += u * RADIXI;
- }
-
- u = (z[3] + TWO71) - TWO71;
- if (u > z[3])
- u -= TWO19;
+ mantissa_store_t d, r;
+ d = X[i] * (mantissa_store_t) a;
+ DIV_RADIX (d, r);
+ z[i] = r;
+ z[i - 1] += d;
+ }
+
+ u = ALIGN_DOWN_TO (z[3], TWO19);
v = z[3] - u;
if (v == TWO18)
{
- if (z[4] == ZERO)
+ if (z[4] == 0)
{
for (i = 5; i <= p; i++)
{
- if (X[i] == ZERO)
+ if (X[i] == 0)
continue;
else
{
- z[3] += ONE;
+ z[3] += 1;
break;
}
}
}
else
- z[3] += ONE;
+ z[3] += 1;
}
c = (z[1] + R * (z[2] + R * z[3])) / a;
@@ -200,12 +198,13 @@
{
long i, k;
long p2 = p;
- double c, u, z[5];
+ double c;
+ mantissa_t u, z[5];
#define R RADIXI
if (EX < -44 || (EX == -44 && X[1] < TWO5))
{
- *y = ZERO;
+ *y = 0;
return;
}
@@ -214,21 +213,21 @@
if (EX == -42)
{
z[1] = X[1] + TWO10;
- z[2] = ZERO;
- z[3] = ZERO;
+ z[2] = 0;
+ z[3] = 0;
k = 3;
}
else if (EX == -43)
{
z[1] = TWO10;
z[2] = X[1];
- z[3] = ZERO;
+ z[3] = 0;
k = 2;
}
else
{
z[1] = TWO10;
- z[2] = ZERO;
+ z[2] = 0;
z[3] = X[1];
k = 1;
}
@@ -239,7 +238,7 @@
{
z[1] = X[1] + TWO10;
z[2] = X[2];
- z[3] = ZERO;
+ z[3] = 0;
k = 3;
}
else if (EX == -43)
@@ -252,7 +251,7 @@
else
{
z[1] = TWO10;
- z[2] = ZERO;
+ z[2] = 0;
z[3] = X[1];
k = 1;
}
@@ -274,25 +273,23 @@
else
{
z[1] = TWO10;
- z[2] = ZERO;
+ z[2] = 0;
k = 1;
}
z[3] = X[k];
}
- u = (z[3] + TWO57) - TWO57;
- if (u > z[3])
- u -= TWO5;
+ u = ALIGN_DOWN_TO (z[3], TWO5);
if (u == z[3])
{
for (i = k + 1; i <= p2; i++)
{
- if (X[i] == ZERO)
+ if (X[i] == 0)
continue;
else
{
- z[3] += ONE;
+ z[3] += 1;
break;
}
}
@@ -309,9 +306,9 @@
void
__mp_dbl (const mp_no *x, double *y, int p)
{
- if (X[0] == ZERO)
- {
- *y = ZERO;
+ if (X[0] == 0)
+ {
+ *y = 0;
return;
}
@@ -330,41 +327,36 @@
{
long i, n;
long p2 = p;
- double u;
/* Sign. */
- if (x == ZERO)
- {
- Y[0] = ZERO;
- return;
- }
- else if (x > ZERO)
- Y[0] = ONE;
+ if (x == 0)
+ {
+ Y[0] = 0;
+ return;
+ }
+ else if (x > 0)
+ Y[0] = 1;
else
{
- Y[0] = MONE;
+ Y[0] = -1;
x = -x;
}
/* Exponent. */
- for (EY = ONE; x >= RADIX; EY += ONE)
+ for (EY = 1; x >= RADIX; EY += 1)
x *= RADIXI;
- for (; x < ONE; EY -= ONE)
+ for (; x < 1; EY -= 1)
x *= RADIX;
/* Digits. */
n = MIN (p2, 4);
for (i = 1; i <= n; i++)
{
- u = (x + TWO52) - TWO52;
- if (u > x)
- u -= ONE;
- Y[i] = u;
- x -= u;
+ INTEGER_OF (x, Y[i]);
x *= RADIX;
}
for (; i <= p2; i++)
- Y[i] = ZERO;
+ Y[i] = 0;
}
/* Add magnitudes of *X and *Y assuming that abs (*X) >= abs (*Y) > 0. The
@@ -377,7 +369,7 @@
{
long i, j, k;
long p2 = p;
- double zk;
+ mantissa_t zk;
EZ = EX;
@@ -391,7 +383,7 @@
return;
}
- zk = ZERO;
+ zk = 0;
for (; j > 0; i--, j--)
{
@@ -399,12 +391,12 @@
if (zk >= RADIX)
{
Z[k--] = zk - RADIX;
- zk = ONE;
+ zk = 1;
}
else
{
Z[k--] = zk;
- zk = ZERO;
+ zk = 0;
}
}
@@ -414,16 +406,16 @@
if (zk >= RADIX)
{
Z[k--] = zk - RADIX;
- zk = ONE;
+ zk = 1;
}
else
{
Z[k--] = zk;
- zk = ZERO;
- }
- }
-
- if (zk == ZERO)
+ zk = 0;
+ }
+ }
+
+ if (zk == 0)
{
for (i = 1; i <= p2; i++)
Z[i] = Z[i + 1];
@@ -431,7 +423,7 @@
else
{
Z[1] = zk;
- EZ += ONE;
+ EZ += 1;
}
}
@@ -445,7 +437,7 @@
{
long i, j, k;
long p2 = p;
- double zk;
+ mantissa_t zk;
EZ = EX;
i = p2;
@@ -461,27 +453,27 @@
/* The relevant least significant digit in Y is non-zero, so we factor it in
to enhance accuracy. */
- if (j < p2 && Y[j + 1] > ZERO)
+ if (j < p2 && Y[j + 1] > 0)
{
Z[k + 1] = RADIX - Y[j + 1];
- zk = MONE;
+ zk = -1;
}
else
- zk = Z[k + 1] = ZERO;
+ zk = Z[k + 1] = 0;
/* Subtract and borrow. */
for (; j > 0; i--, j--)
{
zk += (X[i] - Y[j]);
- if (zk < ZERO)
+ if (zk < 0)
{
Z[k--] = zk + RADIX;
- zk = MONE;
+ zk = -1;
}
else
{
Z[k--] = zk;
- zk = ZERO;
+ zk = 0;
}
}
@@ -489,25 +481,25 @@
for (; i > 0; i--)
{
zk += X[i];
- if (zk < ZERO)
+ if (zk < 0)
{
Z[k--] = zk + RADIX;
- zk = MONE;
+ zk = -1;
}
else
{
Z[k--] = zk;
- zk = ZERO;
+ zk = 0;
}
}
/* Normalize. */
- for (i = 1; Z[i] == ZERO; i++);
+ for (i = 1; Z[i] == 0; i++);
EZ = EZ - i + 1;
for (k = 1; i <= p2 + 1;)
Z[k++] = Z[i++];
for (; k <= p2;)
- Z[k++] = ZERO;
+ Z[k++] = 0;
}
/* Add *X and *Y and store the result in *Z. X and Y may overlap, but not X
@@ -519,12 +511,12 @@
{
int n;
- if (X[0] == ZERO)
+ if (X[0] == 0)
{
__cpy (y, z, p);
return;
}
- else if (Y[0] == ZERO)
+ else if (Y[0] == 0)
{
__cpy (x, z, p);
return;
@@ -556,7 +548,7 @@
Z[0] = Y[0];
}
else
- Z[0] = ZERO;
+ Z[0] = 0;
}
}
@@ -569,13 +561,13 @@
{
int n;
- if (X[0] == ZERO)
+ if (X[0] == 0)
{
__cpy (y, z, p);
Z[0] = -Z[0];
return;
}
- else if (Y[0] == ZERO)
+ else if (Y[0] == 0)
{
__cpy (x, z, p);
return;
@@ -607,7 +599,7 @@
Z[0] = -Y[0];
}
else
- Z[0] = ZERO;
+ Z[0] = 0;
}
}
@@ -621,28 +613,28 @@
{
long i, j, k, ip, ip2;
long p2 = p;
- double u, zk;
+ mantissa_store_t zk;
const mp_no *a;
- double *diag;
+ mantissa_store_t *diag;
/* Is z=0? */
- if (__glibc_unlikely (X[0] * Y[0] == ZERO))
- {
- Z[0] = ZERO;
+ if (__glibc_unlikely (X[0] * Y[0] == 0))
+ {
+ Z[0] = 0;
return;
}
/* We need not iterate through all X's and Y's since it's pointless to
multiply zeroes. Here, both are zero... */
for (ip2 = p2; ip2 > 0; ip2--)
- if (X[ip2] != ZERO || Y[ip2] != ZERO)
+ if (X[ip2] != 0 || Y[ip2] != 0)
break;
- a = X[ip2] != ZERO ? y : x;
+ a = X[ip2] != 0 ? y : x;
/* ... and here, at least one of them is still zero. */
for (ip = ip2; ip > 0; ip--)
- if (a->d[ip] != ZERO)
+ if (a->d[ip] != 0)
break;
/* The product looks like this for p = 3 (as an example):
@@ -669,22 +661,22 @@
Another thing that becomes evident is that only the most significant
ip+ip2 digits of the result are non-zero, where ip and ip2 are the
'internal precision' of the input numbers, i.e. digits after ip and ip2
- are all ZERO. */
+ are all 0. */
k = (__glibc_unlikely (p2 < 3)) ? p2 + p2 : p2 + 3;
while (k > ip + ip2 + 1)
- Z[k--] = ZERO;
-
- zk = ZERO;
+ Z[k--] = 0;
+
+ zk = 0;
/* Precompute sums of diagonal elements so that we can directly use them
later. See the next comment to know we why need them. */
- diag = alloca (k * sizeof (double));
- double d = ZERO;
+ diag = alloca (k * sizeof (mantissa_store_t));
+ mantissa_store_t d = 0;
for (i = 1; i <= ip; i++)
{
- d += X[i] * Y[i];
+ d += X[i] * (mantissa_store_t) Y[i];
diag[i] = d;
}
while (i < k)
@@ -697,18 +689,15 @@
if (k % 2 == 0)
/* We want to add this only once, but since we subtract it in the sum
of products above, we add twice. */
- zk += 2 * X[lim] * Y[lim];
+ zk += 2 * X[lim] * (mantissa_store_t) Y[lim];
for (i = k - p2, j = p2; i < j; i++, j--)
- zk += (X[i] + X[j]) * (Y[i] + Y[j]);
+ zk += (X[i] + X[j]) * (mantissa_store_t) (Y[i] + Y[j]);
zk -= diag[k - 1];
- u = (zk + CUTTER) - CUTTER;
- if (u > zk)
- u -= RADIX;
- Z[k--] = zk - u;
- zk = u * RADIXI;
+ DIV_RADIX (zk, Z[k]);
+ k--;
}
/* The real deal. Mantissa digit Z[k] is the sum of all X[i] * Y[j] where i
@@ -731,18 +720,15 @@
if (k % 2 == 0)
/* We want to add this only once, but since we subtract it in the sum
of products above, we add twice. */
- zk += 2 * X[lim] * Y[lim];
+ zk += 2 * X[lim] * (mantissa_store_t) Y[lim];
for (i = 1, j = k - 1; i < j; i++, j--)
- zk += (X[i] + X[j]) * (Y[i] + Y[j]);
+ zk += (X[i] + X[j]) * (mantissa_store_t) (Y[i] + Y[j]);
zk -= diag[k - 1];
- u = (zk + CUTTER) - CUTTER;
- if (u > zk)
- u -= RADIX;
- Z[k--] = zk - u;
- zk = u * RADIXI;
+ DIV_RADIX (zk, Z[k]);
+ k--;
}
Z[k] = zk;
@@ -752,7 +738,7 @@
int e = EX + EY;
/* Is there a carry beyond the most significant digit? */
- if (__glibc_unlikely (Z[1] == ZERO))
+ if (__glibc_unlikely (Z[1] == 0))
{
for (i = 1; i <= p2; i++)
Z[i] = Z[i + 1];
@@ -774,35 +760,35 @@
__sqr (const mp_no *x, mp_no *y, int p)
{
long i, j, k, ip;
- double u, yk;
+ mantissa_store_t yk;
/* Is z=0? */
- if (__glibc_unlikely (X[0] == ZERO))
- {
- Y[0] = ZERO;
+ if (__glibc_unlikely (X[0] == 0))
+ {
+ Y[0] = 0;
return;
}
/* We need not iterate through all X's since it's pointless to
multiply zeroes. */
for (ip = p; ip > 0; ip--)
- if (X[ip] != ZERO)
+ if (X[ip] != 0)
break;
k = (__glibc_unlikely (p < 3)) ? p + p : p + 3;
while (k > 2 * ip + 1)
- Y[k--] = ZERO;
-
- yk = ZERO;
+ Y[k--] = 0;
+
+ yk = 0;
while (k > p)
{
- double yk2 = 0.0;
+ mantissa_store_t yk2 = 0;
long lim = k / 2;
if (k % 2 == 0)
- yk += X[lim] * X[lim];
+ yk += X[lim] * (mantissa_store_t) X[lim];
/* In __mul, this loop (and the one within the next while loop) run
between a range to calculate the mantissa as follows:
@@ -814,41 +800,35 @@
result. For cases where the range size is even, the mid-point needs
to be added separately (above). */
for (i = k - p, j = p; i < j; i++, j--)
- yk2 += X[i] * X[j];
-
- yk += 2.0 * yk2;
-
- u = (yk + CUTTER) - CUTTER;
- if (u > yk)
- u -= RADIX;
- Y[k--] = yk - u;
- yk = u * RADIXI;
+ yk2 += X[i] * (mantissa_store_t) X[j];
+
+ yk += 2 * yk2;
+
+ DIV_RADIX (yk, Y[k]);
+ k--;
}
while (k > 1)
{
- double yk2 = 0.0;
+ mantissa_store_t yk2 = 0;
long lim = k / 2;
if (k % 2 == 0)
- yk += X[lim] * X[lim];
+ yk += X[lim] * (mantissa_store_t) X[lim];
/* Likewise for this loop. */
for (i = 1, j = k - 1; i < j; i++, j--)
- yk2 += X[i] * X[j];
-
- yk += 2.0 * yk2;
-
- u = (yk + CUTTER) - CUTTER;
- if (u > yk)
- u -= RADIX;
- Y[k--] = yk - u;
- yk = u * RADIXI;
+ yk2 += X[i] * (mantissa_store_t) X[j];
+
+ yk += 2 * yk2;
+
+ DIV_RADIX (yk, Y[k]);
+ k--;
}
Y[k] = yk;
/* Squares are always positive. */
- Y[0] = 1.0;
+ Y[0] = 1;
/* Get the exponent sum into an intermediate variable. This is a subtle
optimization, where given enough registers, all operations on the exponent
@@ -856,7 +836,7 @@
int e = EX * 2;
/* Is there a carry beyond the most significant digit? */
- if (__glibc_unlikely (Y[1] == ZERO))
+ if (__glibc_unlikely (Y[1] == 0))
{
for (i = 1; i <= p; i++)
Y[i] = Y[i + 1];
@@ -888,7 +868,7 @@
__cpy (x, &z, p);
z.e = 0;
__mp_dbl (&z, &t, p);
- t = ONE / t;
+ t = 1 / t;
__dbl_mp (t, y, p);
EY -= EX;
@@ -914,8 +894,8 @@
{
mp_no w;
- if (X[0] == ZERO)
- Z[0] = ZERO;
+ if (X[0] == 0)
+ Z[0] = 0;
else
{
__inv (y, &w, p);
Modified: fsf/trunk/libc/sysdeps/ieee754/dbl-64/mpa.h
==============================================================================
--- fsf/trunk/libc/sysdeps/ieee754/dbl-64/mpa.h (original)
+++ fsf/trunk/libc/sysdeps/ieee754/dbl-64/mpa.h Wed Mar 27 00:01:57 2013
@@ -35,6 +35,7 @@
/* Common types and definition */
/************************************************************************/
+#include <mpa-arch.h>
/* The mp_no structure holds the details of a multi-precision floating point
number.
@@ -61,7 +62,7 @@
typedef struct
{
int e;
- double d[40];
+ mantissa_t d[40];
} mp_no;
typedef union
@@ -82,9 +83,13 @@
#define ABS(x) ((x) < 0 ? -(x) : (x))
-#define RADIX 0x1.0p24 /* 2^24 */
-#define RADIXI 0x1.0p-24 /* 2^-24 */
-#define CUTTER 0x1.0p76 /* 2^76 */
+#ifndef RADIXI
+# define RADIXI 0x1.0p-24 /* 2^-24 */
+#endif
+
+#ifndef TWO52
+# define TWO52 0x1.0p52 /* 2^52 */
+#endif
#define ZERO 0.0 /* 0 */
#define MZERO -0.0 /* 0 with the sign bit set */
@@ -92,13 +97,13 @@
#define MONE -1.0 /* -1 */
#define TWO 2.0 /* 2 */
-#define TWO5 0x1.0p5 /* 2^5 */
-#define TWO8 0x1.0p8 /* 2^52 */
-#define TWO10 0x1.0p10 /* 2^10 */
-#define TWO18 0x1.0p18 /* 2^18 */
-#define TWO19 0x1.0p19 /* 2^19 */
-#define TWO23 0x1.0p23 /* 2^23 */
-#define TWO52 0x1.0p52 /* 2^52 */
+#define TWO5 TWOPOW (5) /* 2^5 */
+#define TWO8 TWOPOW (8) /* 2^52 */
+#define TWO10 TWOPOW (10) /* 2^10 */
+#define TWO18 TWOPOW (18) /* 2^18 */
+#define TWO19 TWOPOW (19) /* 2^19 */
+#define TWO23 TWOPOW (23) /* 2^23 */
+
#define TWO57 0x1.0p57 /* 2^57 */
#define TWO71 0x1.0p71 /* 2^71 */
#define TWOM1032 0x1.0p-1032 /* 2^-1032 */
Modified: fsf/trunk/libc/sysdeps/ieee754/dbl-64/s_atan.c
==============================================================================
--- fsf/trunk/libc/sysdeps/ieee754/dbl-64/s_atan.c (original)
+++ fsf/trunk/libc/sysdeps/ieee754/dbl-64/s_atan.c Wed Mar 27 00:01:57 2013
@@ -62,18 +62,9 @@
#ifndef DLA_FMS
double t4,t5,t6;
#endif
-#if 0
- double y1,y2;
-#endif
int i,ux,dx;
-#if 0
- int p;
-#endif
static const int pr[M]={6,8,10,32};
number num;
-#if 0
- mp_no mpt1,mpx,mpy,mpy1,mpy2,mperr;
-#endif
num.d = x; ux = num.i[HIGH_HALF]; dx = num.i[LOW_HALF];
Modified: fsf/trunk/libc/sysdeps/ieee754/dbl-64/s_sin.c
==============================================================================
--- fsf/trunk/libc/sysdeps/ieee754/dbl-64/s_sin.c (original)
+++ fsf/trunk/libc/sysdeps/ieee754/dbl-64/s_sin.c Wed Mar 27 00:01:57 2013
@@ -100,14 +100,8 @@
SECTION
__sin(double x){
double xx,res,t,cor,y,s,c,sn,ssn,cs,ccs,xn,a,da,db,eps,xn1,xn2;
-#if 0
- double w[2];
-#endif
mynumber u,v;
int4 k,m,n;
-#if 0
- int4 nn;
-#endif
double retval = 0;
SET_RESTORE_ROUND_53BIT (FE_TONEAREST);
@@ -902,10 +896,6 @@
bsloww(double x,double dx, double orig,int n) {
static const double th2_36 = 206158430208.0; /* 1.5*2**37 */
double y,x1,x2,xx,r,t,res,cor,w[2];
-#if 0
- double a,da,xn;
- union {int4 i[2]; double x;} v;
-#endif
x1=(x+th2_36)-th2_36;
y = aa.x*x1*x1*x1;
r=x+y;
Modified: fsf/trunk/libc/sysdeps/ieee754/dbl-64/s_tan.c
==============================================================================
--- fsf/trunk/libc/sysdeps/ieee754/dbl-64/s_tan.c (original)
+++ fsf/trunk/libc/sysdeps/ieee754/dbl-64/s_tan.c Wed Mar 27 00:01:57 2013
@@ -64,9 +64,6 @@
int p;
number num,v;
mp_no mpa,mpt1,mpt2;
-#if 0
- mp_no mpy;
-#endif
double retval;
Modified: fsf/trunk/libc/sysdeps/powerpc/fpu/s_llround.c
==============================================================================
--- fsf/trunk/libc/sysdeps/powerpc/fpu/s_llround.c (original)
+++ fsf/trunk/libc/sysdeps/powerpc/fpu/s_llround.c Wed Mar 27 00:01:57 2013
@@ -17,6 +17,7 @@
<http://www.gnu.org/licenses/>. */
#include <math.h>
+#include <math_ldbl_opt.h>
/* I think that what this routine is supposed to do is round a value
to the nearest integer, with values exactly on the boundary rounded
@@ -47,3 +48,6 @@
strong_alias (__llround, __llroundl)
weak_alias (__llround, llroundl)
#endif
+#if LONG_DOUBLE_COMPAT (libm, GLIBC_2_1)
+compat_symbol (libm, __llround, llroundl, GLIBC_2_1);
+#endif
Added: fsf/trunk/libc/sysdeps/powerpc/power4/fpu/mpa-arch.h
==============================================================================
--- fsf/trunk/libc/sysdeps/powerpc/power4/fpu/mpa-arch.h (added)
+++ fsf/trunk/libc/sysdeps/powerpc/power4/fpu/mpa-arch.h Wed Mar 27 00:01:57 2013
@@ -1,0 +1,56 @@
+/* Overridable constants and operations.
+ Copyright (C) 2013 Free Software Foundation, Inc.
+
+ This program is free software; you can redistribute it and/or modify
+ it under the terms of the GNU Lesser General Public License as published by
+ the Free Software Foundation; either version 2.1 of the License, or
+ (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU Lesser General Public License for more details.
+
+ You should have received a copy of the GNU Lesser General Public License
+ along with this program; if not, see <http://www.gnu.org/licenses/>. */
+
+typedef double mantissa_t;
+typedef double mantissa_store_t;
+
+#define TWOPOW(i) (0x1.0p##i)
+
+#define RADIX TWOPOW (24) /* 2^24 */
+#define CUTTER TWOPOW (76) /* 2^76 */
+#define RADIXI 0x1.0p-24 /* 2^-24 */
+#define TWO52 TWOPOW (52) /* 2^52 */
+
+/* Divide D by RADIX and put the remainder in R. */
+#define DIV_RADIX(d,r) \
+ ({ \
+ double u = ((d) + CUTTER) - CUTTER; \
+ if (u > (d)) \
+ u -= RADIX; \
+ r = (d) - u; \
+ (d) = u * RADIXI; \
+ })
+
+/* Put the integer component of a double X in R and retain the fraction in
+ X. */
+#define INTEGER_OF(x, r) \
+ ({ \
+ double u = ((x) + TWO52) - TWO52; \
+ if (u > (x)) \
+ u -= ONE; \
+ (r) = u; \
+ (x) -= u; \
+ })
+
+/* Align IN down to a multiple of F, where F is a power of two. */
+#define ALIGN_DOWN_TO(in, f) \
+ ({ \
+ double factor = f * TWO52; \
+ double u = (in + factor) - factor; \
+ if (u > in) \
+ u -= f; \
+ u; \
+ })
_______________________________________________
Commits mailing list
Commits@xxxxxxxxxx
http://eglibc.org/cgi-bin/mailman/listinfo/commits