This is the mail archive of the
libc-alpha@sourceware.org
mailing list for the glibc project.
[PATCH] Convert mantissa storage in mp_no to int
- From: Siddhesh Poyarekar <siddhesh at redhat dot com>
- To: libc-alpha at sourceware dot org
- Date: Fri, 21 Dec 2012 00:43:20 +0530
- Subject: [PATCH] Convert mantissa storage in mp_no to int
Hi,
The mantissa in mp_no is currently stored as an array of doubles.
This is unnecessary and slow since:
1) The largest number stored is 2^24
2) Only integers are stored.
The side-effect of the storage being double is that the generated code
is all fp arithmetic, which is a lost opportunity for optimization.
Attached patch converts this double array to int and adjusts mp code
accordingly. This patch is not exhaustive in terms of converting all
possible fp operations to integer operations, but it ensures that
correctness of operations is maintained. I'll be following up with
patches for individual functions to tighten them up and eliminate fp
operations wherever possible.
I have verified that the testsuite runs clean after this change. Like
my last patch, I tested the patch for performance improvement with the
help of the old trusty pow function test in the following blog post:
http://entropymine.com/imageworsener/slowpow/
and the command line:
time ./powtest 1000000 1.0000000000000020 1.5000000050000000
The patch takes a whole second off the execution time for the above
command.
OK for 2.18?
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..5ea13b8 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];
+ long tmp = Z[k];
+ for (i=i1,j=i2-1; i<i2; i++,j--) tmp += (long) 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 */