[PATCH v2] Convert mantissa storage in mp_no to int
Siddhesh Poyarekar
siddhesh@redhat.com
Fri Dec 21 02:48:00 GMT 2012
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 */
More information about the Libc-alpha
mailing list