This is the mail archive of the
libc-alpha@sourceware.org
mailing list for the glibc project.
[PATCH v2] Convert mantissa storage in mp_no to int
Hi,
Here's v2 of the patch with the problem of truncation of an
intermediate value fixed as Joseph pointed out. The code is not used
by i386, so I couldn't test it to see it break and then fixed.
Siddhesh
ChangeLog:
2012-12-21 Siddhesh Poyarekar <siddhesh@redhat.com>
* sysdeps/ieee754/dbl-64/mpa.c (I_RADIX): New macro.
(TWOPOW): Likewise.
(__acr): Juggle logic a bit to make it faster and use int
instead of doubles.
(__mp_dbl): Use int values instead of double.
(__dbl_mp): Likewise.
(add_magnitudes): Likewise.
(sub_magnitudes): Likewise.
(__add): Likewise.
(__sub): Likewise.
(__mul): Likewise.
* sysdeps/ieee754/dbl-64/mpa.h (mp_no): Change type of D to int
array.
* sysdeps/ieee754/dbl-64/mpa.c (__mpatan): Use int values
instead of double.
* sysdeps/ieee754/dbl-64/mpatan2.c (__mpatan2): Likewise.
* sysdeps/ieee754/dbl-64/mpexp.c (__mpexp): Likewise.
* sysdeps/ieee754/dbl-64/mpsqrt.c (__mpsqrt): Likewise.
* sysdeps/ieee754/dbl-64/mptan.c (__mptan): Likewise.
* sysdeps/ieee754/dbl-64/sincos32.c (__c32): Likewise.
(__mpranred): Likewise.
* sysdeps/ieee754/dbl-64/sincos32.h (oofac27): Likewise.
(pi): Likewise.
(hp): Likewise.
* sysdeps/ieee754/dbl-64/slowpow.c (__slowpow): Likewise.
diff --git a/sysdeps/ieee754/dbl-64/mpa.c b/sysdeps/ieee754/dbl-64/mpa.c
index b5d25ed..127ee62 100644
--- a/sysdeps/ieee754/dbl-64/mpa.c
+++ b/sysdeps/ieee754/dbl-64/mpa.c
@@ -1,7 +1,7 @@
/*
* IBM Accurate Mathematical Library
* written by International Business Machines Corp.
- * Copyright (C) 2001, 2011 Free Software Foundation
+ * Copyright (C) 2001-2012 Free Software Foundation
*
* 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
@@ -47,6 +47,9 @@
#include "mpa2.h"
#include <sys/param.h> /* For MIN() */
+#define I_RADIX (1 << 24)
+#define TWOPOW(n) (1 << n)
+
#ifndef SECTION
# define SECTION
#endif
@@ -76,16 +79,12 @@ int
__acr(const mp_no *x, const mp_no *y, int p) {
int i;
- if (X[0] == ZERO) {
- if (Y[0] == ZERO) i= 0;
- else i=-1;
- }
- else if (Y[0] == ZERO) i= 1;
- else {
- if (EX > EY) i= 1;
- else if (EX < EY) i=-1;
- else i= mcr(x,y,p);
- }
+ if (X[0] == 0 || Y[0] == 0)
+ return X[0] - Y[0];
+
+ if (EX > EY) i= 1;
+ else if (EX < EY) i=-1;
+ else i= mcr(x,y,p);
return i;
}
@@ -245,10 +244,10 @@ void __mp_dbl(const mp_no *x, double *y, int p) {
double a,c,u,v,z[5];
#endif
- if (X[0] == ZERO) {*y = ZERO; return; }
+ if (X[0] == 0) {*y = ZERO; return; }
if (EX> -42) norm(x,y,p);
- else if (EX==-42 && X[1]>=TWO10) norm(x,y,p);
+ else if (EX==-42 && X[1]>=TWOPOW(10)) norm(x,y,p);
else denorm(x,y,p);
}
#endif
@@ -263,24 +262,21 @@ SECTION
__dbl_mp(double x, mp_no *y, int p) {
int i,n;
- double u;
/* Sign */
- if (x == ZERO) {Y[0] = ZERO; return; }
- else if (x > ZERO) Y[0] = ONE;
- else {Y[0] = MONE; x=-x; }
+ if (x == ZERO) {Y[0] = 0; return; }
+ else if (x > ZERO) Y[0] = 1;
+ else {Y[0] = -1; x=-x; }
/* Exponent */
- for (EY=ONE; x >= RADIX; EY += ONE) x *= RADIXI;
- for ( ; x < ONE; EY -= ONE) x *= RADIX;
+ for (EY=1; x >= RADIX; EY++) x *= RADIXI;
+ for ( ; x < ONE; EY--) x *= RADIX;
/* Digits */
n=MIN(p,4);
for (i=1; i<=n; i++) {
- u = (x + TWO52) - TWO52;
- if (u>x) u -= ONE;
- Y[i] = u; x -= u; x *= RADIX; }
- for ( ; i<=p; i++) Y[i] = ZERO;
+ Y[i] = x; x -= Y[i]; x *= RADIX; }
+ for ( ; i<=p; i++) Y[i] = 0;
}
@@ -302,29 +298,29 @@ add_magnitudes(const mp_no *x, const mp_no *y, mp_no *z, int p) {
if (j<1)
{__cpy(x,z,p); return; }
- else Z[k] = ZERO;
+ else Z[k] = 0;
for (; j>0; i--,j--) {
Z[k] += X[i] + Y[j];
- if (Z[k] >= RADIX) {
- Z[k] -= RADIX;
- Z[--k] = ONE; }
+ if (Z[k] >= I_RADIX) {
+ Z[k] -= I_RADIX;
+ Z[--k] = 1; }
else
- Z[--k] = ZERO;
+ Z[--k] = 0;
}
for (; i>0; i--) {
Z[k] += X[i];
- if (Z[k] >= RADIX) {
- Z[k] -= RADIX;
- Z[--k] = ONE; }
+ if (Z[k] >= I_RADIX) {
+ Z[k] -= I_RADIX;
+ Z[--k] = 1; }
else
- Z[--k] = ZERO;
+ Z[--k] = 0;
}
- if (Z[1] == ZERO) {
+ if (Z[1] == 0) {
for (i=1; i<=p; i++) Z[i] = Z[i+1]; }
- else EZ += ONE;
+ else EZ++;
}
@@ -344,45 +340,45 @@ sub_magnitudes(const mp_no *x, const mp_no *y, mp_no *z, int p) {
if (EX == EY) {
i=j=k=p;
- Z[k] = Z[k+1] = ZERO; }
+ Z[k] = Z[k+1] = 0; }
else {
j= EX - EY;
if (j > p) {__cpy(x,z,p); return; }
else {
i=p; j=p+1-j; k=p;
- if (Y[j] > ZERO) {
- Z[k+1] = RADIX - Y[j--];
- Z[k] = MONE; }
+ if (Y[j] > 0) {
+ Z[k+1] = I_RADIX - Y[j--];
+ Z[k] = -1; }
else {
- Z[k+1] = ZERO;
- Z[k] = ZERO; j--;}
+ Z[k+1] = 0;
+ Z[k] = 0; j--;}
}
}
for (; j>0; i--,j--) {
Z[k] += (X[i] - Y[j]);
- if (Z[k] < ZERO) {
- Z[k] += RADIX;
- Z[--k] = MONE; }
+ if (Z[k] < 0) {
+ Z[k] += I_RADIX;
+ Z[--k] = -1; }
else
- Z[--k] = ZERO;
+ Z[--k] = 0;
}
for (; i>0; i--) {
Z[k] += X[i];
- if (Z[k] < ZERO) {
- Z[k] += RADIX;
- Z[--k] = MONE; }
+ if (Z[k] < 0) {
+ Z[k] += I_RADIX;
+ Z[--k] = -1; }
else
- Z[--k] = ZERO;
+ Z[--k] = 0;
}
- for (i=1; Z[i] == ZERO; i++) ;
+ for (i=1; Z[i] == 0; i++) ;
EZ = EZ - i + 1;
for (k=1; i <= p+1; )
Z[k++] = Z[i++];
for (; k <= p; )
- Z[k++] = ZERO;
+ Z[k++] = 0;
}
@@ -396,8 +392,8 @@ __add(const mp_no *x, const mp_no *y, mp_no *z, int p) {
int n;
- if (X[0] == ZERO) {__cpy(y,z,p); return; }
- else if (Y[0] == ZERO) {__cpy(x,z,p); return; }
+ if (X[0] == 0) {__cpy(y,z,p); return; }
+ else if (Y[0] == 0) {__cpy(x,z,p); return; }
if (X[0] == Y[0]) {
if (__acr(x,y,p) > 0) {add_magnitudes(x,y,z,p); Z[0] = X[0]; }
@@ -406,7 +402,7 @@ __add(const mp_no *x, const mp_no *y, mp_no *z, int p) {
else {
if ((n=__acr(x,y,p)) == 1) {sub_magnitudes(x,y,z,p); Z[0] = X[0]; }
else if (n == -1) {sub_magnitudes(y,x,z,p); Z[0] = Y[0]; }
- else Z[0] = ZERO;
+ else Z[0] = 0;
}
}
@@ -421,8 +417,8 @@ __sub(const mp_no *x, const mp_no *y, mp_no *z, int p) {
int n;
- if (X[0] == ZERO) {__cpy(y,z,p); Z[0] = -Z[0]; return; }
- else if (Y[0] == ZERO) {__cpy(x,z,p); return; }
+ if (X[0] == 0) {__cpy(y,z,p); Z[0] = -Z[0]; return; }
+ else if (Y[0] == 0) {__cpy(x,z,p); return; }
if (X[0] != Y[0]) {
if (__acr(x,y,p) > 0) {add_magnitudes(x,y,z,p); Z[0] = X[0]; }
@@ -431,7 +427,7 @@ __sub(const mp_no *x, const mp_no *y, mp_no *z, int p) {
else {
if ((n=__acr(x,y,p)) == 1) {sub_magnitudes(x,y,z,p); Z[0] = X[0]; }
else if (n == -1) {sub_magnitudes(y,x,z,p); Z[0] = -Y[0]; }
- else Z[0] = ZERO;
+ else Z[0] = 0;
}
}
@@ -446,28 +442,26 @@ SECTION
__mul(const mp_no *x, const mp_no *y, mp_no *z, int p) {
int i, i1, i2, j, k, k2;
- double u;
/* Is z=0? */
- if (X[0]*Y[0]==ZERO)
- { Z[0]=ZERO; return; }
+ if (X[0]*Y[0] == 0)
+ { Z[0]=0; return; }
/* Multiply, add and carry */
k2 = (p<3) ? p+p : p+3;
- Z[k2]=ZERO;
+ Z[k2]=0;
for (k=k2; k>1; ) {
if (k > p) {i1=k-p; i2=p+1; }
else {i1=1; i2=k; }
- for (i=i1,j=i2-1; i<i2; i++,j--) Z[k] += X[i]*Y[j];
+ int64_t tmp = Z[k];
+ for (i=i1,j=i2-1; i<i2; i++,j--) tmp += (int64_t) X[i]*Y[j];
- u = (Z[k] + CUTTER)-CUTTER;
- if (u > Z[k]) u -= RADIX;
- Z[k] -= u;
- Z[--k] = u*RADIXI;
+ Z[k] = (int) (tmp % (1 << 24));
+ Z[--k] = (int) (tmp / (1 << 24));
}
/* Is there a carry beyond the most significant digit? */
- if (Z[1] == ZERO) {
+ if (Z[1] == 0) {
for (i=1; i<=p; i++) Z[i]=Z[i+1];
EZ = EX + EY - 1; }
else
diff --git a/sysdeps/ieee754/dbl-64/mpa.h b/sysdeps/ieee754/dbl-64/mpa.h
index 4fdecb6..c09689a 100644
--- a/sysdeps/ieee754/dbl-64/mpa.h
+++ b/sysdeps/ieee754/dbl-64/mpa.h
@@ -1,7 +1,7 @@
/*
* IBM Accurate Mathematical Library
* Written by International Business Machines Corp.
- * Copyright (C) 2001, 2011 Free Software Foundation, Inc.
+ * Copyright (C) 2001-2012 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
@@ -41,7 +41,7 @@
typedef struct {/* This structure holds the details of a multi-precision */
int e; /* floating point number, x: d[0] holds its sign (-1,0 or 1) */
- double d[40]; /* e holds its exponent (...,-2,-1,0,1,2,...) and */
+ int d[40]; /* e holds its exponent (...,-2,-1,0,1,2,...) and */
} mp_no; /* d[1]...d[p] hold its mantissa digits. The value of x is, */
/* x = d[1]*r**(e-1) + d[2]*r**(e-2) + ... + d[p]*r**(e-p). */
/* Here r = 2**24, 0 <= d[i] < r and 1 <= p <= 32. */
diff --git a/sysdeps/ieee754/dbl-64/mpatan.c b/sysdeps/ieee754/dbl-64/mpatan.c
index d897bbb..9119e58 100644
--- a/sysdeps/ieee754/dbl-64/mpatan.c
+++ b/sysdeps/ieee754/dbl-64/mpatan.c
@@ -1,7 +1,7 @@
/*
* IBM Accurate Mathematical Library
* written by International Business Machines Corp.
- * Copyright (C) 2001, 2011 Free Software Foundation
+ * Copyright (C) 2001-2012 Free Software Foundation
*
* 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
@@ -48,12 +48,12 @@ __mpatan(mp_no *x, mp_no *y, int p) {
int i,m,n;
double dx;
mp_no
- mptwo = {0,{0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,
- 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,
- 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0}},
- mptwoim1 = {0,{0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,
- 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,
- 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0}};
+ mptwo = {0,{0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+ 0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+ 0,0,0,0,0,0,0,0,0,0,0,0}},
+ mptwoim1 = {0,{0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+ 0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+ 0,0,0,0,0,0,0,0,0,0,0,0}};
mp_no mps,mpsm,mpt,mpt1,mpt2,mpt3;
@@ -66,8 +66,8 @@ __mpatan(mp_no *x, mp_no *y, int p) {
{if (dx>__atan_xm[m].d) break;}
}
mptwo.e = mptwoim1.e = 1;
- mptwo.d[0] = mptwoim1.d[0] = ONE;
- mptwo.d[1] = TWO;
+ mptwo.d[0] = mptwoim1.d[0] = 1;
+ mptwo.d[1] = 2;
/* Reduce x m times */
__mul(x,x,&mpsm,p);
@@ -89,7 +89,7 @@ __mpatan(mp_no *x, mp_no *y, int p) {
n=__atan_np[p]; mptwoim1.d[1] = __atan_twonm1[p].d;
__dvd(&mpsm,&mptwoim1,&mpt,p);
for (i=n-1; i>1; i--) {
- mptwoim1.d[1] -= TWO;
+ mptwoim1.d[1] -= 2;
__dvd(&mpsm,&mptwoim1,&mpt1,p);
__mul(&mpsm,&mpt,&mpt2,p);
__sub(&mpt1,&mpt2,&mpt,p);
diff --git a/sysdeps/ieee754/dbl-64/mpatan2.c b/sysdeps/ieee754/dbl-64/mpatan2.c
index 42cb6a7..3e3252b 100644
--- a/sysdeps/ieee754/dbl-64/mpatan2.c
+++ b/sysdeps/ieee754/dbl-64/mpatan2.c
@@ -1,7 +1,7 @@
/*
* IBM Accurate Mathematical Library
* written by International Business Machines Corp.
- * Copyright (C) 2001, 2011 Free Software Foundation
+ * Copyright (C) 2001-2012 Free Software Foundation
*
* 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
@@ -49,14 +49,12 @@ void
SECTION
__mpatan2(mp_no *y, mp_no *x, mp_no *z, int p) {
- static const double ZERO = 0.0, ONE = 1.0;
-
mp_no mpt1,mpt2,mpt3;
- if (X[0] <= ZERO) {
+ if (X[0] <= 0) {
__dvd(x,y,&mpt1,p); __mul(&mpt1,&mpt1,&mpt2,p);
- if (mpt1.d[0] != ZERO) mpt1.d[0] = ONE;
+ if (mpt1.d[0] != 0) mpt1.d[0] = 1;
__add(&mpt2,&mpone,&mpt3,p); __mpsqrt(&mpt3,&mpt2,p);
__add(&mpt1,&mpt2,&mpt3,p); mpt3.d[0]=Y[0];
__mpatan(&mpt3,&mpt1,p); __add(&mpt1,&mpt1,z,p);
diff --git a/sysdeps/ieee754/dbl-64/mpexp.c b/sysdeps/ieee754/dbl-64/mpexp.c
index 2762ba9..87391dc 100644
--- a/sysdeps/ieee754/dbl-64/mpexp.c
+++ b/sysdeps/ieee754/dbl-64/mpexp.c
@@ -1,7 +1,7 @@
/*
* IBM Accurate Mathematical Library
* written by International Business Machines Corp.
- * Copyright (C) 2001, 2011 Free Software Foundation
+ * Copyright (C) 2001-2012 Free Software Foundation
*
* 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
@@ -56,9 +56,9 @@ __mpexp(mp_no *x, mp_no *y, int p) {
{ 0, 0, 0, 0, 0, 0,23,28,33,38,42,47,52,57,62,66, 0, 0},
{ 0, 0, 0, 0, 0, 0, 0, 0,27, 0, 0,39,43,47,51,55,59,63},
{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,43,47,50,54}};
- mp_no mpk = {0,{0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,
- 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,
- 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0}};
+ mp_no mpk = {0,{0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+ 0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+ 0,0,0,0,0,0,0,0,0,0,0,0}};
mp_no mps,mpak,mpt1,mpt2;
/* Choose m,n and compute a=2**(-m) */
@@ -81,7 +81,7 @@ __mpexp(mp_no *x, mp_no *y, int p) {
__mul(x,&mpt1,&mps,p);
/* Evaluate the polynomial. Put result in mpt2 */
- mpk.e = 1; mpk.d[0] = ONE; mpk.d[1]=__mpexp_nn[n].d;
+ mpk.e = 1; mpk.d[0] = 1; mpk.d[1]=__mpexp_nn[n].d;
__dvd(&mps,&mpk,&mpt1,p);
__add(&mpone,&mpt1,&mpak,p);
for (k=n-1; k>1; k--) {
diff --git a/sysdeps/ieee754/dbl-64/mpsqrt.c b/sysdeps/ieee754/dbl-64/mpsqrt.c
index 92bf5ef..e29105b 100644
--- a/sysdeps/ieee754/dbl-64/mpsqrt.c
+++ b/sysdeps/ieee754/dbl-64/mpsqrt.c
@@ -1,7 +1,7 @@
/*
* IBM Accurate Mathematical Library
* written by International Business Machines Corp.
- * Copyright (C) 2001, 2011 Free Software Foundation
+ * Copyright (C) 2001-2012 Free Software Foundation
*
* 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
@@ -53,17 +53,17 @@ __mpsqrt(mp_no *x, mp_no *y, int p) {
int i,m,ey;
double dx,dy;
mp_no
- mphalf = {0,{0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,
- 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,
- 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0}},
- mp3halfs = {0,{0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,
- 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,
- 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0}};
+ mphalf = {0,{0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+ 0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+ 0,0,0,0,0,0,0,0,0,0,0,0}},
+ mp3halfs = {0,{0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+ 0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+ 0,0,0,0,0,0,0,0,0,0,0,0}};
mp_no mpxn,mpz,mpu,mpt1,mpt2;
/* Prepare multi-precision 1/2 and 3/2 */
- mphalf.e =0; mphalf.d[0] =ONE; mphalf.d[1] =HALFRAD;
- mp3halfs.e=1; mp3halfs.d[0]=ONE; mp3halfs.d[1]=ONE; mp3halfs.d[2]=HALFRAD;
+ mphalf.e =0; mphalf.d[0] =1; mphalf.d[1] =(1 << 23);
+ mp3halfs.e=1; mp3halfs.d[0]=1; mp3halfs.d[1]=1; mp3halfs.d[2]=(1 << 23);
ey=EX/2; __cpy(x,&mpxn,p); mpxn.e -= (ey+ey);
__mp_dbl(&mpxn,&dx,p); dy=fastiroot(dx); __dbl_mp(dy,&mpu,p);
diff --git a/sysdeps/ieee754/dbl-64/mptan.c b/sysdeps/ieee754/dbl-64/mptan.c
index e75b2da..356bc04 100644
--- a/sysdeps/ieee754/dbl-64/mptan.c
+++ b/sysdeps/ieee754/dbl-64/mptan.c
@@ -1,7 +1,7 @@
/*
* IBM Accurate Mathematical Library
* written by International Business Machines Corp.
- * Copyright (C) 2001, 2011 Free Software Foundation
+ * Copyright (C) 2001-2012 Free Software Foundation
*
* 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
@@ -47,8 +47,6 @@ void
SECTION
__mptan(double x, mp_no *mpy, int p) {
- static const double MONE = -1.0;
-
int n;
mp_no mpw, mpc, mps;
@@ -56,7 +54,7 @@ __mptan(double x, mp_no *mpy, int p) {
__c32(&mpw, &mpc, &mps, p); /* computing sin(x) and cos(x) */
if (n) /* second or fourth quarter of unit circle */
{ __dvd(&mpc,&mps,mpy,p);
- mpy->d[0] *= MONE;
+ mpy->d[0] *= -1;
} /* tan is negative in this area */
else __dvd(&mps,&mpc,mpy,p);
diff --git a/sysdeps/ieee754/dbl-64/sincos32.c b/sysdeps/ieee754/dbl-64/sincos32.c
index 038d1bd..3890815 100644
--- a/sysdeps/ieee754/dbl-64/sincos32.c
+++ b/sysdeps/ieee754/dbl-64/sincos32.c
@@ -1,7 +1,7 @@
/*
* IBM Accurate Mathematical Library
* written by International Business Machines Corp.
- * Copyright (C) 2001, 2011 Free Software Foundation
+ * Copyright (C) 2001-2012 Free Software Foundation
*
* 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
@@ -119,7 +119,7 @@ cc32(mp_no *x, mp_no *y, int p) {
void
SECTION
__c32(mp_no *x, mp_no *y, mp_no *z, int p) {
- static const mp_no mpt={1,{1.0,2.0}}, one={1,{1.0,1.0}};
+ static const mp_no mpt={1,{1,2}}, one={1,{1,1}};
mp_no u,t,t1,t2,c,s;
int i;
__cpy(x,&u,p);
@@ -251,7 +251,7 @@ __mpranred(double x, mp_no *y, int p)
number v;
double t,xn;
int i,k,n;
- static const mp_no one = {1,{1.0,1.0}};
+ static const mp_no one = {1,{1,1}};
mp_no a,b,c;
if (ABS(x) < 2.8e14) {
@@ -267,18 +267,18 @@ __mpranred(double x, mp_no *y, int p)
}
else { /* if x is very big more precision required */
__dbl_mp(x,&a,p);
- a.d[0]=1.0;
+ a.d[0]=1;
k = a.e-5;
if (k < 0) k=0;
b.e = -k;
- b.d[0] = 1.0;
+ b.d[0] = 1;
for (i=0;i<p;i++) b.d[i+1] = toverp[i+k];
__mul(&a,&b,&c,p);
t = c.d[c.e];
for (i=1;i<=p-c.e;i++) c.d[i]=c.d[i+c.e];
for (i=p+1-c.e;i<=p;i++) c.d[i]=0;
c.e=0;
- if (c.d[1] >= 8388608.0)
+ if (c.d[1] >= 8388608)
{ t +=1.0;
__sub(&c,&one,&b,p);
__mul(&b,&hp,y,p);
diff --git a/sysdeps/ieee754/dbl-64/sincos32.h b/sysdeps/ieee754/dbl-64/sincos32.h
index 6efe4d4..14562fb 100644
--- a/sysdeps/ieee754/dbl-64/sincos32.h
+++ b/sysdeps/ieee754/dbl-64/sincos32.h
@@ -1,7 +1,7 @@
/*
* IBM Accurate Mathematical Library
* Written by International Business Machines Corp.
- * Copyright (C) 2001 Free Software Foundation, Inc.
+ * Copyright (C) 2001-2012 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
@@ -42,24 +42,24 @@ static const number
#endif
static const mp_no
- oofac27 = {-3,{1.0,7.0,4631664.0,12006312.0,13118056.0,6538613.0,646354.0,
- 8508025.0,9131256.0,7548776.0,2529842.0,8864927.0,660489.0,15595125.0,12777885.0,
- 11618489.0,13348664.0,5486686.0,514518.0,11275535.0,4727621.0,3575562.0,
- 13579710.0,5829745.0,7531862.0,9507898.0,6915060.0,4079264.0,1907586.0,
- 6078398.0,13789314.0,5504104.0,14136.0}},
- pi = {1,{1.0,3.0,
- 2375530.0,8947107.0,578323.0,1673774.0,225395.0,4498441.0,3678761.0,
- 10432976.0,536314.0,10021966.0,7113029.0,2630118.0,3723283.0,7847508.0,
- 6737716.0,15273068.0,12626985.0,12044668.0,5299519.0,8705461.0,11880201.0,
- 1544726.0,14014857.0,7994139.0,13709579.0,10918111.0,11906095.0,16610011.0,
- 13638367.0,12040417.0,11529578.0,2522774.0}},
- hp = {1,{1.0, 1.0,
- 9576373.0,4473553.0,8677769.0,9225495.0,112697.0,10637828.0,
- 10227988.0,13605096.0,268157.0,5010983.0,3556514.0,9703667.0,
- 1861641.0,12312362.0,3368858.0,7636534.0,6313492.0,14410942.0,
- 2649759.0,12741338.0,14328708.0,9160971.0,7007428.0,12385677.0,
- 15243397.0,13847663.0,14341655.0,16693613.0,15207791.0,14408816.0,
- 14153397.0,1261387.0,6110792.0,2291862.0,4181138.0,5295267.0}};
+ oofac27 = {-3,{1,7,4631664,12006312,13118056,6538613,646354,
+ 8508025,9131256,7548776,2529842,8864927,660489,15595125,12777885,
+ 11618489,13348664,5486686,514518,11275535,4727621,3575562,
+ 13579710,5829745,7531862,9507898,6915060,4079264,1907586,
+ 6078398,13789314,5504104,14136}},
+ pi = {1,{1,3,
+ 2375530,8947107,578323,1673774,225395,4498441,3678761,
+ 10432976,536314,10021966,7113029,2630118,3723283,7847508,
+ 6737716,15273068,12626985,12044668,5299519,8705461,11880201,
+ 1544726,14014857,7994139,13709579,10918111,11906095,16610011,
+ 13638367,12040417,11529578,2522774}},
+ hp = {1,{1, 1,
+ 9576373,4473553,8677769,9225495,112697,10637828,
+ 10227988,13605096,268157,5010983,3556514,9703667,
+ 1861641,12312362,3368858,7636534,6313492,14410942,
+ 2649759,12741338,14328708,9160971,7007428,12385677,
+ 15243397,13847663,14341655,16693613,15207791,14408816,
+ 14153397,1261387,6110792,2291862,4181138,5295267}};
static const double toverp[75] = {
10680707.0, 7228996.0, 1387004.0, 2578385.0, 16069853.0,
diff --git a/sysdeps/ieee754/dbl-64/slowpow.c b/sysdeps/ieee754/dbl-64/slowpow.c
index 7c829c7..be2808f 100644
--- a/sysdeps/ieee754/dbl-64/slowpow.c
+++ b/sysdeps/ieee754/dbl-64/slowpow.c
@@ -1,7 +1,7 @@
/*
* IBM Accurate Mathematical Library
* written by International Business Machines Corp.
- * Copyright (C) 2001, 2011 Free Software Foundation
+ * Copyright (C) 2001-2012 Free Software Foundation
*
* 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
@@ -48,7 +48,7 @@ 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}};
+ static const mp_no eps = {-3,{1,4}};
int p;
res = __halfulp(x,y); /* halfulp() returns -10 or x^y */