summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--ChangeLog24
-rw-r--r--sysdeps/ieee754/dbl-64/branred.c2
-rw-r--r--sysdeps/ieee754/dbl-64/doasin.c26
-rw-r--r--sysdeps/ieee754/dbl-64/dosincos.c14
-rw-r--r--sysdeps/ieee754/dbl-64/e_asin.c4
-rw-r--r--sysdeps/ieee754/dbl-64/e_atan2.c12
-rw-r--r--sysdeps/ieee754/dbl-64/e_lgamma_r.c2
-rw-r--r--sysdeps/ieee754/dbl-64/e_log.c6
-rw-r--r--sysdeps/ieee754/dbl-64/e_pow.c44
-rw-r--r--sysdeps/ieee754/dbl-64/e_remainder.c6
-rw-r--r--sysdeps/ieee754/dbl-64/halfulp.c52
-rw-r--r--sysdeps/ieee754/dbl-64/mpa.c58
-rw-r--r--sysdeps/ieee754/dbl-64/mpa.h20
-rw-r--r--sysdeps/ieee754/dbl-64/mpatan.c38
-rw-r--r--sysdeps/ieee754/dbl-64/mpatan2.c18
-rw-r--r--sysdeps/ieee754/dbl-64/mpexp.c28
-rw-r--r--sysdeps/ieee754/dbl-64/mplog.c16
-rw-r--r--sysdeps/ieee754/dbl-64/mpsqrt.c20
-rw-r--r--sysdeps/ieee754/dbl-64/mptan.c10
-rw-r--r--sysdeps/ieee754/dbl-64/s_atan.c6
-rw-r--r--sysdeps/ieee754/dbl-64/s_sin.c57
-rw-r--r--sysdeps/ieee754/dbl-64/s_tan.c12
-rw-r--r--sysdeps/ieee754/dbl-64/sincos32.c120
-rw-r--r--sysdeps/ieee754/dbl-64/slowexp.c20
-rw-r--r--sysdeps/ieee754/dbl-64/slowpow.c38
25 files changed, 350 insertions, 303 deletions
diff --git a/ChangeLog b/ChangeLog
index 964320481d..0741b0f34b 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,10 +1,26 @@
+2001-03-12 Ulrich Drepper <drepper@redhat.com>
+
+ * sysdeps/ieee754/dbl-64/e_remainder.c: Fix handling of boundary
+ conditions.
+
+ * sysdeps/ieee754/dbl-64/e_pow.c: Fix handling of boundary
+ conditions.
+
+ * sysdeps/ieee754/dbl-64/s_sin.c (__sin): Handle Inf and NaN
+ correctly.
+ (__cos): Likewise.
+
+ * sysdeps/ieee754/dbl-64/e_asin.c (__ieee754_asin): Handle NaN
+ correctly.
+ (__ieee754_acos): Likewise.
+
2001-03-12 Andreas Jaeger <aj@suse.de>
* sysdeps/unix/sysv/linux/s390/sysdep.h (_LINUX_S390_SYSDEP_H):
Fix typo. Patch by Martin Schwidefsky <schwidefsky@de.ibm.com>.
* sysdeps/s390/bits/string.h: Protect __STRING_INLINE against
- redefinition.
+ redefinition.
2001-03-11 Roland McGrath <roland@frob.com>
@@ -12,6 +28,12 @@
2001-03-11 Ulrich Drepper <drepper@redhat.com>
+ * sysdeps/ieee754/dbl-64/endian.h: Define also one of BIG_ENDI and
+ LITTLE_ENDI.
+
+ * sysdeps/ieee754/dbl-64/MathLib.h (Init_Lib): Use void as
+ parameter list.
+
Last-bit accurate math library implementation by IBM Haifa.
Contributed by Abraham Ziv <ziv@il.ibm.com>, Moshe Olshansky
<olshansk@il.ibm.com>, Ealan Henis <ealan@il.ibm.com>, and
diff --git a/sysdeps/ieee754/dbl-64/branred.c b/sysdeps/ieee754/dbl-64/branred.c
index 2bceb947c0..3253d34aba 100644
--- a/sysdeps/ieee754/dbl-64/branred.c
+++ b/sysdeps/ieee754/dbl-64/branred.c
@@ -45,7 +45,7 @@
/* x=n*pi/2+(a+aa), abs(a+aa)<pi/4, n=0,+-1,+-2,.... */
/* Routine return integer (n mod 4) */
/*******************************************************************/
-int branred(double x, double *a, double *aa)
+int __branred(double x, double *a, double *aa)
{
int i,k;
#if 0
diff --git a/sysdeps/ieee754/dbl-64/doasin.c b/sysdeps/ieee754/dbl-64/doasin.c
index 4ea108b9d7..d3cc88bb19 100644
--- a/sysdeps/ieee754/dbl-64/doasin.c
+++ b/sysdeps/ieee754/dbl-64/doasin.c
@@ -5,9 +5,9 @@
*
* 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 of the License, or
+ * the Free Software Foundation; either version 2 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
@@ -15,12 +15,12 @@
*
* You should have received a copy of the GNU Lesser General Public License
* along with this program; if not, write to the Free Software
- * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
+ * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
*/
/**********************************************************************/
/* MODULE_NAME: doasin.c */
/* */
-/* FUNCTION: doasin */
+/* FUNCTION: doasin */
/* */
/* FILES NEEDED:endian.h mydefs.h dla.h doasin.h */
/* mpa.c */
@@ -31,17 +31,17 @@
#include "endian.h"
#include "mydefs.h"
-#include "dla.h"
+#include "dla.h"
/********************************************************************/
/* Compute arcsin(x,dx,v) of double-length number (x+dx) the result */
/* stored in v where v= v[0]+v[1] =arcsin(x+dx) */
/********************************************************************/
-void doasin(double x, double dx, double v[]) {
-
+void __doasin(double x, double dx, double v[]) {
+
#include "doasin.h"
- static const double
+ static const double
d5 = 0.22372159090911789889975459505194491E-01,
d6 = 0.17352764422456822913014975683014622E-01,
d7 = 0.13964843843786693521653681033981614E-01,
@@ -49,10 +49,10 @@ void doasin(double x, double dx, double v[]) {
d9 = 0.97622386568166960207425666787248914E-02,
d10 = 0.83638737193775788576092749009744976E-02,
d11 = 0.79470250400727425881446981833568758E-02;
-
+
double xx,p,pp,u,uu,r,s;
double hx,tx,hy,ty,tp,tq,tc,tcc;
-
+
/* Taylor series for arcsin for Double-Length numbers */
xx = x*x+2.0*x*dx;
@@ -73,9 +73,3 @@ void doasin(double x, double dx, double v[]) {
v[0]=p;
v[1]=pp; /* arcsin(x+dx)=v[0]+v[1] */
}
-
-
-
-
-
-
diff --git a/sysdeps/ieee754/dbl-64/dosincos.c b/sysdeps/ieee754/dbl-64/dosincos.c
index 692fc85889..498cc60929 100644
--- a/sysdeps/ieee754/dbl-64/dosincos.c
+++ b/sysdeps/ieee754/dbl-64/dosincos.c
@@ -45,7 +45,7 @@
/*(x+dx) between 0 and PI/4 */
/***********************************************************************/
-void dubsin(double x, double dx, double v[]) {
+void __dubsin(double x, double dx, double v[]) {
double r,s,p,hx,tx,hy,ty,q,c,cc,d,dd,d2,dd2,e,ee,
sn,ssn,cs,ccs,ds,dss,dc,dcc;
#if 0
@@ -96,7 +96,7 @@ void dubsin(double x, double dx, double v[]) {
/*(x+dx) between 0 and PI/4 */
/**********************************************************************/
-void dubcos(double x, double dx, double v[]) {
+void __dubcos(double x, double dx, double v[]) {
double r,s,p,hx,tx,hy,ty,q,c,cc,d,dd,d2,dd2,e,ee,
sn,ssn,cs,ccs,ds,dss,dc,dcc;
#if 0
@@ -159,20 +159,20 @@ void dubcos(double x, double dx, double v[]) {
/* Routine receive Double-Length number (x+dx) and computes cos(x+dx) */
/* as Double-Length number and store it in array v */
/**********************************************************************/
-void docos(double x, double dx, double v[]) {
+void __docos(double x, double dx, double v[]) {
double y,yy,p,w[2];
if (x>0) {y=x; yy=dx;}
else {y=-x; yy=-dx;}
if (y<0.5*hp0.x) /* y< PI/4 */
- {dubcos(y,yy,w); v[0]=w[0]; v[1]=w[1];}
+ {__dubcos(y,yy,w); v[0]=w[0]; v[1]=w[1];}
else if (y<1.5*hp0.x) { /* y< 3/4 * PI */
p=hp0.x-y; /* p = PI/2 - y */
yy=hp1.x-yy;
y=p+yy;
yy=(p-y)+yy;
- if (y>0) {dubsin(y,yy,w); v[0]=w[0]; v[1]=w[1];}
+ if (y>0) {__dubsin(y,yy,w); v[0]=w[0]; v[1]=w[1];}
/* cos(x) = sin ( 90 - x ) */
- else {dubsin(-y,-yy,w); v[0]=-w[0]; v[1]=-w[1];
+ else {__dubsin(-y,-yy,w); v[0]=-w[0]; v[1]=-w[1];
}
}
else { /* y>= 3/4 * PI */
@@ -180,7 +180,7 @@ void docos(double x, double dx, double v[]) {
yy=2.0*hp1.x-yy;
y=p+yy;
yy=(p-y)+yy;
- dubcos(y,yy,w);
+ __dubcos(y,yy,w);
v[0]=-w[0];
v[1]=-w[1];
}
diff --git a/sysdeps/ieee754/dbl-64/e_asin.c b/sysdeps/ieee754/dbl-64/e_asin.c
index 2efa63166e..a2b309e4d5 100644
--- a/sysdeps/ieee754/dbl-64/e_asin.c
+++ b/sysdeps/ieee754/dbl-64/e_asin.c
@@ -312,6 +312,8 @@ double __ieee754_asin(double x){
} /* else if (k < 0x3ff00000) */
/*---------------------------- |x|>=1 -------------------------------*/
else if (k==0x3ff00000 && u.i[LOW_HALF]==0) return (m>0)?hp0.x:-hp0.x;
+ else
+ if (k>0x7ff00000 || (k == 0x7ff00000 && u.i[LOW_HALF] != 0)) return x;
else {
u.i[HIGH_HALF]=0x7ff00000;
v.i[HIGH_HALF]=0x7ff00000;
@@ -621,6 +623,8 @@ double __ieee754_acos(double x)
/*---------------------------- |x|>=1 -----------------------*/
else
if (k==0x3ff00000 && u.i[LOW_HALF]==0) return (m>0)?0:2.0*hp0.x;
+ else
+ if (k>0x7ff00000 || (k == 0x7ff00000 && u.i[LOW_HALF] != 0)) return x;
else {
u.i[HIGH_HALF]=0x7ff00000;
v.i[HIGH_HALF]=0x7ff00000;
diff --git a/sysdeps/ieee754/dbl-64/e_atan2.c b/sysdeps/ieee754/dbl-64/e_atan2.c
index b77505db36..adf7a0d11b 100644
--- a/sysdeps/ieee754/dbl-64/e_atan2.c
+++ b/sysdeps/ieee754/dbl-64/e_atan2.c
@@ -374,9 +374,9 @@ static double normalized(double ax,double ay,double y, double z)
{ int p;
mp_no mpx,mpy,mpz,mperr,mpz2,mpt1;
p=6;
- dbl_mp(ax,&mpx,p); dbl_mp(ay,&mpy,p); dvd(&mpy,&mpx,&mpz,p);
- dbl_mp(ue.d,&mpt1,p); mul(&mpz,&mpt1,&mperr,p);
- sub(&mpz,&mperr,&mpz2,p); __mp_dbl(&mpz2,&z,p);
+ __dbl_mp(ax,&mpx,p); __dbl_mp(ay,&mpy,p); __dvd(&mpy,&mpx,&mpz,p);
+ __dbl_mp(ue.d,&mpt1,p); __mul(&mpz,&mpt1,&mperr,p);
+ __sub(&mpz,&mperr,&mpz2,p); __mp_dbl(&mpz2,&z,p);
return signArctan2(y,z);
}
/* Fix the sign and return after stage 1 or stage 2 */
@@ -392,10 +392,10 @@ static double atan2Mp(double x,double y,const int pr[])
mp_no mpx,mpy,mpz,mpz1,mpz2,mperr,mpt1;
for (i=0; i<MM; i++) {
p = pr[i];
- dbl_mp(x,&mpx,p); dbl_mp(y,&mpy,p);
+ __dbl_mp(x,&mpx,p); __dbl_mp(y,&mpy,p);
__mpatan2(&mpy,&mpx,&mpz,p);
- dbl_mp(ud[i].d,&mpt1,p); mul(&mpz,&mpt1,&mperr,p);
- add(&mpz,&mperr,&mpz1,p); sub(&mpz,&mperr,&mpz2,p);
+ __dbl_mp(ud[i].d,&mpt1,p); __mul(&mpz,&mpt1,&mperr,p);
+ __add(&mpz,&mperr,&mpz1,p); __sub(&mpz,&mperr,&mpz2,p);
__mp_dbl(&mpz1,&z1,p); __mp_dbl(&mpz2,&z2,p);
if (z1==z2) return z1;
}
diff --git a/sysdeps/ieee754/dbl-64/e_lgamma_r.c b/sysdeps/ieee754/dbl-64/e_lgamma_r.c
index e5b402a023..863eaa4362 100644
--- a/sysdeps/ieee754/dbl-64/e_lgamma_r.c
+++ b/sysdeps/ieee754/dbl-64/e_lgamma_r.c
@@ -175,7 +175,7 @@ static double zero= 0.00000000000000000000e+00;
GET_HIGH_WORD(ix,x);
ix &= 0x7fffffff;
- if(ix<0x3fd00000) return __kernel_sin(pi*x,zero,0);
+ if(ix<0x3fd00000) return sin(pi*x);
y = -x; /* x is assume negative */
/*
diff --git a/sysdeps/ieee754/dbl-64/e_log.c b/sysdeps/ieee754/dbl-64/e_log.c
index 814ac13720..1a099726d5 100644
--- a/sysdeps/ieee754/dbl-64/e_log.c
+++ b/sysdeps/ieee754/dbl-64/e_log.c
@@ -189,10 +189,10 @@ double __ieee754_log(double x) {
for (i=0; i<M; i++) {
p = pr[i];
- dbl_mp(x,&mpx,p); dbl_mp(y,&mpy,p);
+ __dbl_mp(x,&mpx,p); __dbl_mp(y,&mpy,p);
__mplog(&mpx,&mpy,p);
- dbl_mp(e[i].d,&mperr,p);
- add(&mpy,&mperr,&mpy1,p); sub(&mpy,&mperr,&mpy2,p);
+ __dbl_mp(e[i].d,&mperr,p);
+ __add(&mpy,&mperr,&mpy1,p); __sub(&mpy,&mperr,&mpy2,p);
__mp_dbl(&mpy1,&y1,p); __mp_dbl(&mpy2,&y2,p);
if (y1==y2) return y1;
}
diff --git a/sysdeps/ieee754/dbl-64/e_pow.c b/sysdeps/ieee754/dbl-64/e_pow.c
index cd6f27f33e..dc92922d18 100644
--- a/sysdeps/ieee754/dbl-64/e_pow.c
+++ b/sysdeps/ieee754/dbl-64/e_pow.c
@@ -45,7 +45,7 @@
double __exp1(double x, double xx, double error);
static double log1(double x, double *delta, double *error);
static double log2(double x, double *delta, double *error);
-double slowpow(double x, double y,double z);
+double __slowpow(double x, double y,double z);
static double power1(double x, double y);
static int checkint(double x);
@@ -69,8 +69,8 @@ double __ieee754_pow(double x, double y) {
if (((qx==0x7ff00000)&&(u.i[LOW_HALF]!=0))||(qx>0x7ff00000)) return NaNQ.x;
if (y == 1.0) return x;
if (y == 2.0) return x*x;
- if (y == -1.0) return (x!=0)?1.0/x:NaNQ.x;
- if (y == 0) return ((x>0)&&(qx<0x7ff00000))?1.0:NaNQ.x;
+ if (y == -1.0) return 1.0/x;
+ if (y == 0) return 1.0;
}
/* else */
if(((u.i[HIGH_HALF]>0 && u.i[HIGH_HALF]<0x7ff00000)|| /* x>0 and not x->0 */
@@ -94,16 +94,37 @@ double __ieee754_pow(double x, double y) {
}
if (x == 0) {
- if (ABS(y) > 1.0e20) return (y>0)?0:NaNQ.x;
+ if (((v.i[HIGH_HALF] & 0x7fffffff) == 0x7ff00000 && v.i[LOW_HALF] != 0)
+ || (v.i[HIGH_HALF] & 0x7fffffff) > 0x7ff00000)
+ return y;
+ if (ABS(y) > 1.0e20) return (y>0)?0:INF.x;
k = checkint(y);
- if (k == 0 || y < 0) return NaNQ.x;
- else return (k==1)?0:x; /* return 0 */
+ if (k == -1)
+ return y < 0 ? 1.0/x : x;
+ else
+ return y < 0 ? 1.0/ABS(x) : 0.0; /* return 0 */
}
/* if x<0 */
if (u.i[HIGH_HALF] < 0) {
k = checkint(y);
- if (k==0) return NaNQ.x; /* y not integer and x<0 */
- return (k==1)?upow(-x,y):-upow(-x,y); /* if y even or odd */
+ if (k==0) {
+ if ((v.i[HIGH_HALF] & 0x7fffffff) == 0x7ff00000 && v.i[LOW_HALF] == 0) {
+ if (x == -1.0) return 1.0;
+ else if (x > -1.0) return v.i[HIGH_HALF] < 0 ? INF.x : 0.0;
+ else return v.i[HIGH_HALF] < 0 ? 0.0 : INF.x;
+ }
+ else if (u.i[HIGH_HALF] == 0xfff00000 && u.i[LOW_HALF] == 0)
+ return y < 0 ? 0.0 : INF.x;
+ return NaNQ.x; /* y not integer and x<0 */
+ }
+ else if (u.i[HIGH_HALF] == 0xfff00000 && u.i[LOW_HALF] == 0)
+ {
+ if (k < 0)
+ return y < 0 ? nZERO.x : nINF.x;
+ else
+ return y < 0 ? 0.0 : INF.x;
+ }
+ return (k==1)?__ieee754_pow(-x,y):-__ieee754_pow(-x,y); /* if y even or odd */
}
/* x>0 */
qx = u.i[HIGH_HALF]&0x7fffffff; /* no sign */
@@ -111,7 +132,8 @@ double __ieee754_pow(double x, double y) {
if (qx > 0x7ff00000 || (qx == 0x7ff00000 && u.i[LOW_HALF] != 0)) return NaNQ.x;
/* if 0<x<2^-0x7fe */
- if (qy > 0x7ff00000 || (qy == 0x7ff00000 && v.i[LOW_HALF] != 0)) return NaNQ.x;
+ if (qy > 0x7ff00000 || (qy == 0x7ff00000 && v.i[LOW_HALF] != 0))
+ return x == 1.0 ? 1.0 : NaNQ.x;
/* if y<2^-0x7fe */
if (qx == 0x7ff00000) /* x= 2^-0x3ff */
@@ -124,7 +146,7 @@ double __ieee754_pow(double x, double y) {
if (y<0) return (x<1.0)?INF.x:0;
}
- if (x == 1.0) return NaNQ.x;
+ if (x == 1.0) return 1.0;
if (y>0) return (x>1.0)?INF.x:0;
if (y<0) return (x<1.0)?INF.x:0;
return 0; /* unreachable, to make the compiler happy */
@@ -148,7 +170,7 @@ static double power1(double x, double y) {
a2 = (a-a1)+aa;
error = error*ABS(y);
t = __exp1(a1,a2,1.9e16*error);
- return (t >= 0)?t:slowpow(x,y,z);
+ return (t >= 0)?t:__slowpow(x,y,z);
}
/****************************************************************************/
diff --git a/sysdeps/ieee754/dbl-64/e_remainder.c b/sysdeps/ieee754/dbl-64/e_remainder.c
index f025fdf1b5..c59e5895d8 100644
--- a/sysdeps/ieee754/dbl-64/e_remainder.c
+++ b/sysdeps/ieee754/dbl-64/e_remainder.c
@@ -103,11 +103,13 @@ double __ieee754_remainder(double x, double y)
else {
if (kx<0x7ff00000&&ky<0x7ff00000&&(ky>0||t.i[LOW_HALF]!=0)) {
y=ABS(y)*t128.x;
- z=uremainder(x,y)*t128.x;
- z=uremainder(z,y)*tm128.x;
+ z=__ieee754_remainder(x,y)*t128.x;
+ z=__ieee754_remainder(z,y)*tm128.x;
return z;
}
else { /* if x is too big */
+ if (kx == 0x7ff00000 && u.i[LOW_HALF] == 0 && y == 1.0)
+ return x / x;
if (kx>=0x7ff00000||(ky==0&&t.i[LOW_HALF]==0)||ky>0x7ff00000||
(ky==0x7ff00000&&t.i[LOW_HALF]!=0))
return (u.i[HIGH_HALF]&0x80000000)?nNAN.x:NAN.x;
diff --git a/sysdeps/ieee754/dbl-64/halfulp.c b/sysdeps/ieee754/dbl-64/halfulp.c
index 7901ec4e36..929ca91c50 100644
--- a/sysdeps/ieee754/dbl-64/halfulp.c
+++ b/sysdeps/ieee754/dbl-64/halfulp.c
@@ -5,9 +5,9 @@
*
* 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 of the License, or
+ * the Free Software Foundation; either version 2 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
@@ -15,12 +15,12 @@
*
* You should have received a copy of the GNU Lesser General Public License
* along with this program; if not, write to the Free Software
- * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
+ * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
*/
/************************************************************************/
/* */
-/* MODULE_NAME:halfulp.c */
-/* */
+/* MODULE_NAME:halfulp.c */
+/* */
/* FUNCTIONS:halfulp */
/* FILES NEEDED: mydefs.h dla.h endian.h */
/* uroot.c */
@@ -35,11 +35,11 @@
/*3. if x can be represented by x=2**n for some integer n. */
/************************************************************************/
-#include "endian.h"
+#include "endian.h"
#include "mydefs.h"
#include "dla.h"
-double usqrt(double x);
+double __ieee754_sqrt(double x);
int4 tab54[32] = {
262143, 11585, 1782, 511, 210, 107, 63, 42,
@@ -48,17 +48,17 @@ int4 tab54[32] = {
3, 3, 3, 3, 3, 3, 3, 3 };
-double halfulp(double x, double y)
+double __halfulp(double x, double y)
{
mynumber v;
double z,u,uu,j1,j2,j3,j4,j5;
int4 k,l,m,n;
if (y <= 0) { /*if power is negative or zero */
v.x = y;
- if (v.i[LOW_HALF] != 0) return -10.0;
+ if (v.i[LOW_HALF] != 0) return -10.0;
v.x = x;
- if (v.i[LOW_HALF] != 0) return -10.0;
- if ((v.i[HIGH_HALF]&0x000fffff) != 0) return -10; /* if x =2 ^ n */
+ if (v.i[LOW_HALF] != 0) return -10.0;
+ if ((v.i[HIGH_HALF]&0x000fffff) != 0) return -10; /* if x =2 ^ n */
k = ((v.i[HIGH_HALF]&0x7fffffff)>>20)-1023; /* find this n */
z = (double) k;
return (z*y == -1075.0)?0: -10.0;
@@ -66,19 +66,19 @@ double halfulp(double x, double y)
/* if y > 0 */
v.x = y;
if (v.i[LOW_HALF] != 0) return -10.0;
-
+
v.x=x;
- /* case where x = 2**n for some integer n */
+ /* case where x = 2**n for some integer n */
if (((v.i[HIGH_HALF]&0x000fffff)|v.i[LOW_HALF]) == 0) {
k=(v.i[HIGH_HALF]>>20)-1023;
return (((double) k)*y == -1075.0)?0:-10.0;
- }
-
+ }
+
v.x = y;
k = v.i[HIGH_HALF];
m = k<<12;
l = 0;
- while (m)
+ while (m)
{m = m<<1; l++; }
n = (k&0x000fffff)|0x00100000;
n = n>>(20-l); /* n is the odd integer of y */
@@ -88,19 +88,19 @@ double halfulp(double x, double y)
if (n > 34) return -10.0;
k = -k;
if (k>5) return -10.0;
-
+
/* now treat x */
while (k>0) {
- z = usqrt(x);
+ z = __ieee754_sqrt(x);
EMULV(z,z,u,uu,j1,j2,j3,j4,j5);
if (((u-x)+uu) != 0) break;
x = z;
k--;
}
- if (k) return -10.0;
-
+ if (k) return -10.0;
+
/* it is impossible that n == 2, so the mantissa of x must be short */
-
+
v.x = x;
if (v.i[LOW_HALF]) return -10.0;
k = v.i[HIGH_HALF];
@@ -109,16 +109,14 @@ double halfulp(double x, double y)
while (m) {m = m<<1; l++; }
m = (k&0x000fffff)|0x00100000;
m = m>>(20-l); /* m is the odd integer of x */
-
+
/* now check whether the length of m**n is at most 54 bits */
-
+
if (m > tab54[n-3]) return -10.0;
-
+
/* yes, it is - now compute x**n by simple multiplications */
-
+
u = x;
for (k=1;k<n;k++) u = u*x;
return u;
}
-
-
diff --git a/sysdeps/ieee754/dbl-64/mpa.c b/sysdeps/ieee754/dbl-64/mpa.c
index 1586d91bc2..e9840b037f 100644
--- a/sysdeps/ieee754/dbl-64/mpa.c
+++ b/sysdeps/ieee754/dbl-64/mpa.c
@@ -62,7 +62,7 @@ static int mcr(const mp_no *x, const mp_no *y, int p) {
/* acr() compares the absolute values of two multiple precision numbers */
-int acr(const mp_no *x, const mp_no *y, int p) {
+int __acr(const mp_no *x, const mp_no *y, int p) {
int i;
if (X[0] == ZERO) {
@@ -81,20 +81,20 @@ int acr(const mp_no *x, const mp_no *y, int p) {
/* cr90 compares the values of two multiple precision numbers */
-int cr(const mp_no *x, const mp_no *y, int p) {
+int __cr(const mp_no *x, const mp_no *y, int p) {
int i;
if (X[0] > Y[0]) i= 1;
else if (X[0] < Y[0]) i=-1;
- else if (X[0] < ZERO ) i= acr(y,x,p);
- else i= acr(x,y,p);
+ else if (X[0] < ZERO ) i= __acr(y,x,p);
+ else i= __acr(x,y,p);
return i;
}
/* Copy a multiple precision number. Set *y=*x. x=y is permissible. */
-void cpy(const mp_no *x, mp_no *y, int p) {
+void __cpy(const mp_no *x, mp_no *y, int p) {
int i;
EY = EX;
@@ -110,7 +110,7 @@ void cpy(const mp_no *x, mp_no *y, int p) {
/* n<m, the digits of x beyond the n'th are ignored. */
/* x=y is permissible. */
-void cpymn(const mp_no *x, int m, mp_no *y, int n) {
+void __cpymn(const mp_no *x, int m, mp_no *y, int n) {
int i,k;
@@ -246,7 +246,7 @@ void __mp_dbl(const mp_no *x, double *y, int p) {
/* number *y. If the precision p is too small the result is truncated. x is */
/* left unchanged. */
-void dbl_mp(double x, mp_no *y, int p) {
+void __dbl_mp(double x, mp_no *y, int p) {
int i,n;
double u;
@@ -286,7 +286,7 @@ static void add_magnitudes(const mp_no *x, const mp_no *y, mp_no *z, int p) {
i=p; j=p+ EY - EX; k=p+1;
if (j<1)
- {cpy(x,z,p); return; }
+ {__cpy(x,z,p); return; }
else Z[k] = ZERO;
for (; j>0; i--,j--) {
@@ -330,7 +330,7 @@ static void sub_magnitudes(const mp_no *x, const mp_no *y, mp_no *z, int p) {
Z[k] = Z[k+1] = ZERO; }
else {
j= EX - EY;
- if (j > p) {cpy(x,z,p); return; }
+ if (j > p) {__cpy(x,z,p); return; }
else {
i=p; j=p+1-j; k=p;
if (Y[j] > ZERO) {
@@ -375,19 +375,19 @@ static void sub_magnitudes(const mp_no *x, const mp_no *y, mp_no *z, int p) {
/* but not x&z or y&z. One guard digit is used. The error is less than */
/* one ulp. *x & *y are left unchanged. */
-void add(const mp_no *x, const mp_no *y, mp_no *z, int p) {
+void __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] == ZERO) {__cpy(y,z,p); return; }
+ else if (Y[0] == ZERO) {__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]; }
+ if (__acr(x,y,p) > 0) {add_magnitudes(x,y,z,p); Z[0] = X[0]; }
else {add_magnitudes(y,x,z,p); Z[0] = Y[0]; }
}
else {
- if ((n=acr(x,y,p)) == 1) {sub_magnitudes(x,y,z,p); Z[0] = X[0]; }
+ 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;
}
@@ -399,19 +399,19 @@ void add(const mp_no *x, const mp_no *y, mp_no *z, int p) {
/* overlap but not x&z or y&z. One guard digit is used. The error is */
/* less than one ulp. *x & *y are left unchanged. */
-void sub(const mp_no *x, const mp_no *y, mp_no *z, int p) {
+void __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] == ZERO) {__cpy(y,z,p); Z[0] = -Z[0]; return; }
+ else if (Y[0] == ZERO) {__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]; }
+ if (__acr(x,y,p) > 0) {add_magnitudes(x,y,z,p); Z[0] = X[0]; }
else {add_magnitudes(y,x,z,p); Z[0] = -Y[0]; }
}
else {
- if ((n=acr(x,y,p)) == 1) {sub_magnitudes(x,y,z,p); Z[0] = X[0]; }
+ 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;
}
@@ -424,7 +424,7 @@ void sub(const mp_no *x, const mp_no *y, mp_no *z, int p) {
/* truncated to p digits. In case p>3 the error is bounded by 1.001 ulp. */
/* *x & *y are left unchanged. */
-void mul(const mp_no *x, const mp_no *y, mp_no *z, int p) {
+void __mul(const mp_no *x, const mp_no *y, mp_no *z, int p) {
int i, i1, i2, j, k, k2;
double u;
@@ -464,7 +464,7 @@ void mul(const mp_no *x, const mp_no *y, mp_no *z, int p) {
/* 2.001*r**(1-p) for p>3. */
/* *x=0 is not permissible. *x is left unchanged. */
-void inv(const mp_no *x, mp_no *y, int p) {
+void __inv(const mp_no *x, mp_no *y, int p) {
int i;
#if 0
int l;
@@ -478,14 +478,14 @@ void inv(const mp_no *x, mp_no *y, int p) {
0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,
0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0}};
- cpy(x,&z,p); z.e=0; __mp_dbl(&z,&t,p);
- t=ONE/t; dbl_mp(t,y,p); EY -= EX;
+ __cpy(x,&z,p); z.e=0; __mp_dbl(&z,&t,p);
+ t=ONE/t; __dbl_mp(t,y,p); EY -= EX;
for (i=0; i<np1[p]; i++) {
- cpy(y,&w,p);
- mul(x,&w,y,p);
- sub(&mptwo,y,&z,p);
- mul(&w,&z,y,p);
+ __cpy(y,&w,p);
+ __mul(x,&w,y,p);
+ __sub(&mptwo,y,&z,p);
+ __mul(&w,&z,y,p);
}
return;
}
@@ -496,11 +496,11 @@ void inv(const mp_no *x, mp_no *y, int p) {
/* Relative error bound = 2.001*r**(1-p) for p=2, 2.063*r**(1-p) for p=3 */
/* and 3.001*r**(1-p) for p>3. *y=0 is not permissible. */
-void dvd(const mp_no *x, const mp_no *y, mp_no *z, int p) {
+void __dvd(const mp_no *x, const mp_no *y, mp_no *z, int p) {
mp_no w;
if (X[0] == ZERO) Z[0] = ZERO;
- else {inv(y,&w,p); mul(x,&w,z,p);}
+ else {__inv(y,&w,p); __mul(x,&w,z,p);}
return;
}
diff --git a/sysdeps/ieee754/dbl-64/mpa.h b/sysdeps/ieee754/dbl-64/mpa.h
index bb63e69002..d136283993 100644
--- a/sysdeps/ieee754/dbl-64/mpa.h
+++ b/sysdeps/ieee754/dbl-64/mpa.h
@@ -65,14 +65,14 @@ typedef union { int i[2]; double d; } number;
#define MIN(x,y) ((x) < (y) ? (x) : (y))
#define ABS(x) ((x) < 0 ? -(x) : (x))
-int acr(const mp_no *, const mp_no *, int);
-int cr(const mp_no *, const mp_no *, int);
-void cpy(const mp_no *, mp_no *, int);
-void cpymn(const mp_no *, int, mp_no *, int);
+int __acr(const mp_no *, const mp_no *, int);
+int __cr(const mp_no *, const mp_no *, int);
+void __cpy(const mp_no *, mp_no *, int);
+void __cpymn(const mp_no *, int, mp_no *, int);
void __mp_dbl(const mp_no *, double *, int);
-void dbl_mp(double, mp_no *, int);
-void add(const mp_no *, const mp_no *, mp_no *, int);
-void sub(const mp_no *, const mp_no *, mp_no *, int);
-void mul(const mp_no *, const mp_no *, mp_no *, int);
-void inv(const mp_no *, mp_no *, int);
-void dvd(const mp_no *, const mp_no *, mp_no *, int);
+void __dbl_mp(double, mp_no *, int);
+void __add(const mp_no *, const mp_no *, mp_no *, int);
+void __sub(const mp_no *, const mp_no *, mp_no *, int);
+void __mul(const mp_no *, const mp_no *, mp_no *, int);
+void __inv(const mp_no *, mp_no *, int);
+void __dvd(const mp_no *, const mp_no *, mp_no *, int);
diff --git a/sysdeps/ieee754/dbl-64/mpatan.c b/sysdeps/ieee754/dbl-64/mpatan.c
index 416014bab7..0c8556db60 100644
--- a/sysdeps/ieee754/dbl-64/mpatan.c
+++ b/sysdeps/ieee754/dbl-64/mpatan.c
@@ -33,9 +33,9 @@
#include "endian.h"
#include "mpa.h"
-void mpsqrt(mp_no *, mp_no *, int);
+void __mpsqrt(mp_no *, mp_no *, int);
-void mpatan(mp_no *x, mp_no *y, int p) {
+void __mpatan(mp_no *x, mp_no *y, int p) {
#include "mpatan.h"
int i,m,n;
@@ -66,36 +66,36 @@ void mpatan(mp_no *x, mp_no *y, int p) {
mptwo.d[1] = TWO;
/* Reduce x m times */
- mul(x,x,&mpsm,p);
- if (m==0) cpy(x,&mps,p);
+ __mul(x,x,&mpsm,p);
+ if (m==0) __cpy(x,&mps,p);
else {
for (i=0; i<m; i++) {
- add(&mpone,&mpsm,&mpt1,p);
- mpsqrt(&mpt1,&mpt2,p);
- add(&mpt2,&mpt2,&mpt1,p);
- add(&mptwo,&mpsm,&mpt2,p);
- add(&mpt1,&mpt2,&mpt3,p);
- dvd(&mpsm,&mpt3,&mpt1,p);
- cpy(&mpt1,&mpsm,p);
+ __add(&mpone,&mpsm,&mpt1,p);
+ __mpsqrt(&mpt1,&mpt2,p);
+ __add(&mpt2,&mpt2,&mpt1,p);
+ __add(&mptwo,&mpsm,&mpt2,p);
+ __add(&mpt1,&mpt2,&mpt3,p);
+ __dvd(&mpsm,&mpt3,&mpt1,p);
+ __cpy(&mpt1,&mpsm,p);
}
- mpsqrt(&mpsm,&mps,p); mps.d[0] = X[0];
+ __mpsqrt(&mpsm,&mps,p); mps.d[0] = X[0];
}
/* Evaluate a truncated power series for Atan(s) */
n=np[p]; mptwoim1.d[1] = twonm1[p].d;
- dvd(&mpsm,&mptwoim1,&mpt,p);
+ __dvd(&mpsm,&mptwoim1,&mpt,p);
for (i=n-1; i>1; i--) {
mptwoim1.d[1] -= TWO;
- dvd(&mpsm,&mptwoim1,&mpt1,p);
- mul(&mpsm,&mpt,&mpt2,p);
- sub(&mpt1,&mpt2,&mpt,p);
+ __dvd(&mpsm,&mptwoim1,&mpt1,p);
+ __mul(&mpsm,&mpt,&mpt2,p);
+ __sub(&mpt1,&mpt2,&mpt,p);
}
- mul(&mps,&mpt,&mpt1,p);
- sub(&mps,&mpt1,&mpt,p);
+ __mul(&mps,&mpt,&mpt1,p);
+ __sub(&mps,&mpt1,&mpt,p);
/* Compute Atan(x) */
mptwoim1.d[1] = twom[m].d;
- mul(&mptwoim1,&mpt,y,p);
+ __mul(&mptwoim1,&mpt,y,p);
return;
}
diff --git a/sysdeps/ieee754/dbl-64/mpatan2.c b/sysdeps/ieee754/dbl-64/mpatan2.c
index e41350e72a..2d1625b823 100644
--- a/sysdeps/ieee754/dbl-64/mpatan2.c
+++ b/sysdeps/ieee754/dbl-64/mpatan2.c
@@ -37,12 +37,12 @@
#include "mpa.h"
-void mpsqrt(mp_no *, mp_no *, int);
-void mpatan(mp_no *, mp_no *, int);
+void __mpsqrt(mp_no *, mp_no *, int);
+void __mpatan(mp_no *, mp_no *, int);
/* Multi-Precision Atan2(y,x) function subroutine, for p >= 4. */
/* y=0 is not permitted if x<=0. No error messages are given. */
-void mpatan2(mp_no *y, mp_no *x, mp_no *z, int p) {
+void __mpatan2(mp_no *y, mp_no *x, mp_no *z, int p) {
static const double ZERO = 0.0, ONE = 1.0;
@@ -54,15 +54,15 @@ void mpatan2(mp_no *y, mp_no *x, mp_no *z, int p) {
if (X[0] <= ZERO) {
mpone.e = 1; mpone.d[0] = mpone.d[1] = ONE;
- dvd(x,y,&mpt1,p); mul(&mpt1,&mpt1,&mpt2,p);
+ __dvd(x,y,&mpt1,p); __mul(&mpt1,&mpt1,&mpt2,p);
if (mpt1.d[0] != ZERO) mpt1.d[0] = ONE;
- 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);
+ __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);
}
else
- { dvd(y,x,&mpt1,p);
- mpatan(&mpt1,z,p);
+ { __dvd(y,x,&mpt1,p);
+ __mpatan(&mpt1,z,p);
}
return;
diff --git a/sysdeps/ieee754/dbl-64/mpexp.c b/sysdeps/ieee754/dbl-64/mpexp.c
index 8b7c1dc4bb..2f0b826605 100644
--- a/sysdeps/ieee754/dbl-64/mpexp.c
+++ b/sysdeps/ieee754/dbl-64/mpexp.c
@@ -35,7 +35,7 @@
/* Multi-Precision exponential function subroutine (for p >= 4, */
/* 2**(-55) <= abs(x) <= 1024). */
-void mpexp(mp_no *x, mp_no *y, int p) {
+void __mpexp(mp_no *x, mp_no *y, int p) {
int i,j,k,m,m1,m2,n;
double a,b;
@@ -75,30 +75,30 @@ void mpexp(mp_no *x, mp_no *y, int p) {
}
/* Compute s=x*2**(-m). Put result in mps */
- dbl_mp(a,&mpt1,p);
- mul(x,&mpt1,&mps,p);
+ __dbl_mp(a,&mpt1,p);
+ __mul(x,&mpt1,&mps,p);
/* Evaluate the polynomial. Put result in mpt2 */
mpone.e=1; mpone.d[0]=ONE; mpone.d[1]=ONE;
mpk.e = 1; mpk.d[0] = ONE; mpk.d[1]=nn[n].d;
- dvd(&mps,&mpk,&mpt1,p);
- add(&mpone,&mpt1,&mpak,p);
+ __dvd(&mps,&mpk,&mpt1,p);
+ __add(&mpone,&mpt1,&mpak,p);
for (k=n-1; k>1; k--) {
- mul(&mps,&mpak,&mpt1,p);
+ __mul(&mps,&mpak,&mpt1,p);
mpk.d[1]=nn[k].d;
- dvd(&mpt1,&mpk,&mpt2,p);
- add(&mpone,&mpt2,&mpak,p);
+ __dvd(&mpt1,&mpk,&mpt2,p);
+ __add(&mpone,&mpt2,&mpak,p);
}
- mul(&mps,&mpak,&mpt1,p);
- add(&mpone,&mpt1,&mpt2,p);
+ __mul(&mps,&mpak,&mpt1,p);
+ __add(&mpone,&mpt1,&mpt2,p);
/* Raise polynomial value to the power of 2**m. Put result in y */
for (k=0,j=0; k<m; ) {
- mul(&mpt2,&mpt2,&mpt1,p); k++;
+ __mul(&mpt2,&mpt2,&mpt1,p); k++;
if (k==m) { j=1; break; }
- mul(&mpt1,&mpt1,&mpt2,p); k++;
+ __mul(&mpt1,&mpt1,&mpt2,p); k++;
}
- if (j) cpy(&mpt1,y,p);
- else cpy(&mpt2,y,p);
+ if (j) __cpy(&mpt1,y,p);
+ else __cpy(&mpt2,y,p);
return;
}
diff --git a/sysdeps/ieee754/dbl-64/mplog.c b/sysdeps/ieee754/dbl-64/mplog.c
index 5957c43cc2..cca0c7e4ed 100644
--- a/sysdeps/ieee754/dbl-64/mplog.c
+++ b/sysdeps/ieee754/dbl-64/mplog.c
@@ -37,9 +37,9 @@
#include "endian.h"
#include "mpa.h"
-void mpexp(mp_no *, mp_no *, int);
+void __mpexp(mp_no *, mp_no *, int);
-void mplog(mp_no *x, mp_no *y, int p) {
+void __mplog(mp_no *x, mp_no *y, int p) {
#include "mplog.h"
int i,m;
#if 0
@@ -58,14 +58,14 @@ void mplog(mp_no *x, mp_no *y, int p) {
/* Perform m newton iterations to solve for y: exp(y)-x=0. */
/* The iterations formula is: y(n+1)=y(n)+(x*exp(-y(n))-1). */
- cpy(y,&mpt1,p);
+ __cpy(y,&mpt1,p);
for (i=0; i<m; i++) {
mpt1.d[0]=-mpt1.d[0];
- mpexp(&mpt1,&mpt2,p);
- mul(x,&mpt2,&mpt1,p);
- sub(&mpt1,&mpone,&mpt2,p);
- add(y,&mpt2,&mpt1,p);
- cpy(&mpt1,y,p);
+ __mpexp(&mpt1,&mpt2,p);
+ __mul(x,&mpt2,&mpt1,p);
+ __sub(&mpt1,&mpone,&mpt2,p);
+ __add(y,&mpt2,&mpt1,p);
+ __cpy(&mpt1,y,p);
}
return;
}
diff --git a/sysdeps/ieee754/dbl-64/mpsqrt.c b/sysdeps/ieee754/dbl-64/mpsqrt.c
index 3a4632bdc9..824c03f746 100644
--- a/sysdeps/ieee754/dbl-64/mpsqrt.c
+++ b/sysdeps/ieee754/dbl-64/mpsqrt.c
@@ -42,7 +42,7 @@
double fastiroot(double);
-void mpsqrt(mp_no *x, mp_no *y, int p) {
+void __mpsqrt(mp_no *x, mp_no *y, int p) {
#include "mpsqrt.h"
int i,m,ex,ey;
@@ -60,19 +60,19 @@ void mpsqrt(mp_no *x, mp_no *y, int p) {
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;
- ex=EX; ey=EX/2; cpy(x,&mpxn,p); mpxn.e -= (ey+ey);
- __mp_dbl(&mpxn,&dx,p); dy=fastiroot(dx); dbl_mp(dy,&mpu,p);
- mul(&mpxn,&mphalf,&mpz,p);
+ ex=EX; ey=EX/2; __cpy(x,&mpxn,p); mpxn.e -= (ey+ey);
+ __mp_dbl(&mpxn,&dx,p); dy=fastiroot(dx); __dbl_mp(dy,&mpu,p);
+ __mul(&mpxn,&mphalf,&mpz,p);
m=mp[p];
for (i=0; i<m; i++) {
- mul(&mpu,&mpu,&mpt1,p);
- mul(&mpt1,&mpz,&mpt2,p);
- sub(&mp3halfs,&mpt2,&mpt1,p);
- mul(&mpu,&mpt1,&mpt2,p);
- cpy(&mpt2,&mpu,p);
+ __mul(&mpu,&mpu,&mpt1,p);
+ __mul(&mpt1,&mpz,&mpt2,p);
+ __sub(&mp3halfs,&mpt2,&mpt1,p);
+ __mul(&mpu,&mpt1,&mpt2,p);
+ __cpy(&mpt2,&mpu,p);
}
- mul(&mpxn,&mpu,y,p); EY += ey;
+ __mul(&mpxn,&mpu,y,p); EY += ey;
return;
}
diff --git a/sysdeps/ieee754/dbl-64/mptan.c b/sysdeps/ieee754/dbl-64/mptan.c
index cbfbd0260b..b51a9c0901 100644
--- a/sysdeps/ieee754/dbl-64/mptan.c
+++ b/sysdeps/ieee754/dbl-64/mptan.c
@@ -37,23 +37,23 @@
#include "endian.h"
#include "mpa.h"
-int mpranred(double, mp_no *, int);
+int __mpranred(double, mp_no *, int);
void __c32(mp_no *, mp_no *, mp_no *, int);
-void mptan(double x, mp_no *mpy, int p) {
+void __mptan(double x, mp_no *mpy, int p) {
static const double MONE = -1.0;
int n;
mp_no mpw, mpc, mps;
- n = mpranred(x, &mpw, p) & 0x00000001; /* negative or positive result */
+ n = __mpranred(x, &mpw, p) & 0x00000001; /* negative or positive result */
__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);
+ { __dvd(&mpc,&mps,mpy,p);
mpy->d[0] *= MONE;
} /* tan is negative in this area */
- else dvd(&mps,&mpc,mpy,p);
+ else __dvd(&mps,&mpc,mpy,p);
return;
}
diff --git a/sysdeps/ieee754/dbl-64/s_atan.c b/sysdeps/ieee754/dbl-64/s_atan.c
index 0b71c0d5c8..2872c75db7 100644
--- a/sysdeps/ieee754/dbl-64/s_atan.c
+++ b/sysdeps/ieee754/dbl-64/s_atan.c
@@ -214,9 +214,9 @@ static double atanMp(double x,const int pr[]){
for (i=0; i<M; i++) {
p = pr[i];
- dbl_mp(x,&mpx,p); __mpatan(&mpx,&mpy,p);
- dbl_mp(u9[i].d,&mpt1,p); mul(&mpy,&mpt1,&mperr,p);
- add(&mpy,&mperr,&mpy1,p); sub(&mpy,&mperr,&mpy2,p);
+ __dbl_mp(x,&mpx,p); __mpatan(&mpx,&mpy,p);
+ __dbl_mp(u9[i].d,&mpt1,p); __mul(&mpy,&mpt1,&mperr,p);
+ __add(&mpy,&mperr,&mpy1,p); __sub(&mpy,&mperr,&mpy2,p);
__mp_dbl(&mpy1,&y1,p); __mp_dbl(&mpy2,&y2,p);
if (y1==y2) return y1;
}
diff --git a/sysdeps/ieee754/dbl-64/s_sin.c b/sysdeps/ieee754/dbl-64/s_sin.c
index da9ea27d77..ff6cf01bb8 100644
--- a/sysdeps/ieee754/dbl-64/s_sin.c
+++ b/sysdeps/ieee754/dbl-64/s_sin.c
@@ -60,8 +60,8 @@ static const double
cs4 = -4.16666666666664434524222570944589E-02,
cs6 = 1.38888874007937613028114285595617E-03;
-void dubsin(double x, double dx, double w[]);
-void docos(double x, double dx, double w[]);
+void __dubsin(double x, double dx, double w[]);
+void __docos(double x, double dx, double w[]);
double __mpsin(double x, double dx);
double __mpcos(double x, double dx);
double __mpsin1(double x);
@@ -75,7 +75,7 @@ static double sloww2(double x, double dx, double orig, int n);
static double bsloww(double x, double dx, double orig, int n);
static double bsloww1(double x, double dx, double orig, int n);
static double bsloww2(double x, double dx, double orig, int n);
-int branred(double x, double *a, double *aa);
+int __branred(double x, double *a, double *aa);
static double cslow2(double x);
static double csloww(double x, double dx, double orig);
static double csloww1(double x, double dx, double orig);
@@ -84,7 +84,7 @@ static double csloww2(double x, double dx, double orig, int n);
/* An ultimate sin routine. Given an IEEE double machine number x */
/* it computes the correctly rounded (to nearest) value of sin(x) */
/*******************************************************************/
-double sin(double x){
+double __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];
@@ -307,7 +307,7 @@ double sin(double x){
/* -----------------281474976710656 <|x| <2^1024----------------------------*/
else if (k < 0x7ff00000) {
- n = branred(x,&a,&da);
+ n = __branred(x,&a,&da);
switch (n) {
case 0:
if (a*a < 0.01588) return bsloww(a,da,x,n);
@@ -327,7 +327,7 @@ double sin(double x){
} /* else if (k < 0x7ff00000 ) */
/*--------------------- |x| > 2^1024 ----------------------------------*/
- else return NAN.x;
+ else return x / x;
return 0; /* unreachable */
}
@@ -337,7 +337,7 @@ double sin(double x){
/* it computes the correctly rounded (to nearest) value of cos(x) */
/*******************************************************************/
-double cos(double x)
+double __cos(double x)
{
double y,xx,res,t,cor,s,c,sn,ssn,cs,ccs,xn,a,da,db,eps,xn1,xn2;
mynumber u,v;
@@ -548,7 +548,7 @@ double cos(double x)
else if (k < 0x7ff00000) {/* 281474976710656 <|x| <2^1024 */
- n = branred(x,&a,&da);
+ n = __branred(x,&a,&da);
switch (n) {
case 1:
if (a*a < 0.01588) return bsloww(-a,-da,x,n);
@@ -570,7 +570,7 @@ double cos(double x)
- else return NAN.x; /* |x| > 2^1024 */
+ else return x / x; /* |x| > 2^1024 */
return 0;
}
@@ -594,7 +594,7 @@ static const double th2_36 = 206158430208.0; /* 1.5*2**37 */
cor = (r-res)+t;
if (res == res + 1.0007*cor) return res;
else {
- dubsin(ABS(x),0,w);
+ __dubsin(ABS(x),0,w);
if (w[0] == w[0]+1.000000001*w[1]) return (x>0)?w[0]:-w[0];
else return (x>0)?__mpsin(x,0):-__mpsin(-x,0);
}
@@ -631,7 +631,7 @@ static double slow1(double x) {
cor=(y-res)+cor;
if (res == res+1.0005*cor) return (x>0)?res:-res;
else {
- dubsin(ABS(x),0,w);
+ __dubsin(ABS(x),0,w);
if (w[0] == w[0]+1.000000005*w[1]) return (x>0)?w[0]:-w[0];
else return (x>0)?__mpsin(x,0):-__mpsin(-x,0);
}
@@ -679,7 +679,7 @@ static double slow2(double x) {
y=ABS(x)-hp0.x;
y1=y-hp1.x;
y2=(y-y1)-hp1.x;
- docos(y1,y2,w);
+ __docos(y1,y2,w);
if (w[0] == w[0]+1.000000005*w[1]) return (x>0)?w[0]:-w[0];
else return (x>0)?__mpsin(x,0):-__mpsin(-x,0);
}
@@ -709,7 +709,7 @@ static double sloww(double x,double dx, double orig) {
cor = (cor>0)? 1.0005*cor+ABS(orig)*3.1e-30 : 1.0005*cor-ABS(orig)*3.1e-30;
if (res == res + cor) return res;
else {
- (x>0)? dubsin(x,dx,w) : dubsin(-x,-dx,w);
+ (x>0)? __dubsin(x,dx,w) : __dubsin(-x,-dx,w);
cor = (w[1]>0)? 1.000000001*w[1] + ABS(orig)*1.1e-30 : 1.000000001*w[1] - ABS(orig)*1.1e-30;
if (w[0] == w[0]+cor) return (x>0)?w[0]:-w[0];
else {
@@ -725,7 +725,7 @@ static double sloww(double x,double dx, double orig) {
a = t - y;
da = ((t-a)-y)+da;
if (n&2) {a=-a; da=-da;}
- (a>0)? dubsin(a,da,w) : dubsin(-a,-da,w);
+ (a>0)? __dubsin(a,da,w) : __dubsin(-a,-da,w);
cor = (w[1]>0)? 1.000000001*w[1] + ABS(orig)*1.1e-40 : 1.000000001*w[1] - ABS(orig)*1.1e-40;
if (w[0] == w[0]+cor) return (a>0)?w[0]:-w[0];
else return __mpsin1(orig);
@@ -768,7 +768,7 @@ static double sloww1(double x, double dx, double orig) {
cor = (cor>0)? 1.0005*cor+3.1e-30*ABS(orig) : 1.0005*cor-3.1e-30*ABS(orig);
if (res == res + cor) return (x>0)?res:-res;
else {
- dubsin(ABS(x),dx,w);
+ __dubsin(ABS(x),dx,w);
cor = (w[1]>0)? 1.000000005*w[1]+1.1e-30*ABS(orig) : 1.000000005*w[1]-1.1e-30*ABS(orig);
if (w[0] == w[0]+cor) return (x>0)?w[0]:-w[0];
else return __mpsin1(orig);
@@ -811,7 +811,7 @@ static double sloww2(double x, double dx, double orig, int n) {
cor = (cor>0)? 1.0005*cor+3.1e-30*ABS(orig) : 1.0005*cor-3.1e-30*ABS(orig);
if (res == res + cor) return (n&2)?-res:res;
else {
- docos(ABS(x),dx,w);
+ __docos(ABS(x),dx,w);
cor = (w[1]>0)? 1.000000005*w[1]+1.1e-30*ABS(orig) : 1.000000005*w[1]-1.1e-30*ABS(orig);
if (w[0] == w[0]+cor) return (n&2)?-w[0]:w[0];
else return __mpsin1(orig);
@@ -844,7 +844,7 @@ static double bsloww(double x,double dx, double orig,int n) {
cor = (cor>0)? 1.0005*cor+1.1e-24 : 1.0005*cor-1.1e-24;
if (res == res + cor) return res;
else {
- (x>0)? dubsin(x,dx,w) : dubsin(-x,-dx,w);
+ (x>0)? __dubsin(x,dx,w) : __dubsin(-x,-dx,w);
cor = (w[1]>0)? 1.000000001*w[1] + 1.1e-24 : 1.000000001*w[1] - 1.1e-24;
if (w[0] == w[0]+cor) return (x>0)?w[0]:-w[0];
else return (n&1)?__mpcos1(orig):__mpsin1(orig);
@@ -887,7 +887,7 @@ mynumber u;
cor = (cor>0)? 1.0005*cor+1.1e-24 : 1.0005*cor-1.1e-24;
if (res == res + cor) return (x>0)?res:-res;
else {
- dubsin(ABS(x),dx,w);
+ __dubsin(ABS(x),dx,w);
cor = (w[1]>0)? 1.000000005*w[1]+1.1e-24: 1.000000005*w[1]-1.1e-24;
if (w[0] == w[0]+cor) return (x>0)?w[0]:-w[0];
else return (n&1)?__mpcos1(orig):__mpsin1(orig);
@@ -931,7 +931,7 @@ mynumber u;
cor = (cor>0)? 1.0005*cor+1.1e-24 : 1.0005*cor-1.1e-24;
if (res == res + cor) return (n&2)?-res:res;
else {
- docos(ABS(x),dx,w);
+ __docos(ABS(x),dx,w);
cor = (w[1]>0)? 1.000000005*w[1]+1.1e-24 : 1.000000005*w[1]-1.1e-24;
if (w[0] == w[0]+cor) return (n&2)?-w[0]:w[0];
else return (n&1)?__mpsin1(orig):__mpcos1(orig);
@@ -972,7 +972,7 @@ static double cslow2(double x) {
return res;
else {
y=ABS(x);
- docos(y,0,w);
+ __docos(y,0,w);
if (w[0] == w[0]+1.000000005*w[1]) return w[0];
else return __mpcos(x,0);
}
@@ -1004,7 +1004,7 @@ static double csloww(double x,double dx, double orig) {
cor = (cor>0)? 1.0005*cor+ABS(orig)*3.1e-30 : 1.0005*cor-ABS(orig)*3.1e-30;
if (res == res + cor) return res;
else {
- (x>0)? dubsin(x,dx,w) : dubsin(-x,-dx,w);
+ (x>0)? __dubsin(x,dx,w) : __dubsin(-x,-dx,w);
cor = (w[1]>0)? 1.000000001*w[1] + ABS(orig)*1.1e-30 : 1.000000001*w[1] - ABS(orig)*1.1e-30;
if (w[0] == w[0]+cor) return (x>0)?w[0]:-w[0];
else {
@@ -1020,7 +1020,7 @@ static double csloww(double x,double dx, double orig) {
a = t - y;
da = ((t-a)-y)+da;
if (n==1) {a=-a; da=-da;}
- (a>0)? dubsin(a,da,w) : dubsin(-a,-da,w);
+ (a>0)? __dubsin(a,da,w) : __dubsin(-a,-da,w);
cor = (w[1]>0)? 1.000000001*w[1] + ABS(orig)*1.1e-40 : 1.000000001*w[1] - ABS(orig)*1.1e-40;
if (w[0] == w[0]+cor) return (a>0)?w[0]:-w[0];
else return __mpcos1(orig);
@@ -1064,7 +1064,7 @@ static double csloww1(double x, double dx, double orig) {
cor = (cor>0)? 1.0005*cor+3.1e-30*ABS(orig) : 1.0005*cor-3.1e-30*ABS(orig);
if (res == res + cor) return (x>0)?res:-res;
else {
- dubsin(ABS(x),dx,w);
+ __dubsin(ABS(x),dx,w);
cor = (w[1]>0)? 1.000000005*w[1]+1.1e-30*ABS(orig) : 1.000000005*w[1]-1.1e-30*ABS(orig);
if (w[0] == w[0]+cor) return (x>0)?w[0]:-w[0];
else return __mpcos1(orig);
@@ -1109,14 +1109,19 @@ static double csloww2(double x, double dx, double orig, int n) {
cor = (cor>0)? 1.0005*cor+3.1e-30*ABS(orig) : 1.0005*cor-3.1e-30*ABS(orig);
if (res == res + cor) return (n)?-res:res;
else {
- docos(ABS(x),dx,w);
+ __docos(ABS(x),dx,w);
cor = (w[1]>0)? 1.000000005*w[1]+1.1e-30*ABS(orig) : 1.000000005*w[1]-1.1e-30*ABS(orig);
if (w[0] == w[0]+cor) return (n)?-w[0]:w[0];
else return __mpcos1(orig);
}
}
+weak_alias (__cos, cos)
+weak_alias (__sin, sin)
+
#ifdef NO_LONG_DOUBLE
-weak_alias (sin, sinl)
-weak_alias (cos, cosl)
+strong_alias (__sin, __sinl)
+weak_alias (__sin, sinl)
+strong_alias (__cos, __cosl)
+weak_alias (__cos, cosl)
#endif
diff --git a/sysdeps/ieee754/dbl-64/s_tan.c b/sysdeps/ieee754/dbl-64/s_tan.c
index 26d53b1757..7b5dc4f3da 100644
--- a/sysdeps/ieee754/dbl-64/s_tan.c
+++ b/sysdeps/ieee754/dbl-64/s_tan.c
@@ -53,8 +53,8 @@ double tan(double x) {
mp_no mpy;
#endif
- int branred(double, double *, double *);
- int mpranred(double, mp_no *, int);
+ int __branred(double, double *, double *);
+ int __mpranred(double, mp_no *, int);
/* x=+-INF, x=NaN */
num.d = x; ux = num.i[HIGH_HALF];
@@ -361,7 +361,7 @@ double tan(double x) {
/* (---) The case 1e8 < abs(x) < 2**1024 */
/* Range reduction by algorithm iii */
- n = (branred(x,&a,&da)) & 0x00000001;
+ n = (__branred(x,&a,&da)) & 0x00000001;
EADD(a,da,t1,t2) a=t1; da=t2;
if (a<ZERO) {ya=-a; yya=-da; sy=MONE;}
else {ya= a; yya= da; sy= ONE;}
@@ -384,9 +384,9 @@ double tan(double x) {
/* Second stage */
/* Reduction by algorithm iv */
- p=10; n = (mpranred(x,&mpa,p)) & 0x00000001;
- __mp_dbl(&mpa,&a,p); dbl_mp(a,&mpt1,p);
- sub(&mpa,&mpt1,&mpt2,p); __mp_dbl(&mpt2,&da,p);
+ p=10; n = (__mpranred(x,&mpa,p)) & 0x00000001;
+ __mp_dbl(&mpa,&a,p); __dbl_mp(a,&mpt1,p);
+ __sub(&mpa,&mpt1,&mpt2,p); __mp_dbl(&mpt2,&da,p);
MUL2(a,da,a,da,x2,xx2,t1,t2,t3,t4,t5,t6,t7,t8)
c1 = x2*(a15.d+x2*(a17.d+x2*(a19.d+x2*(a21.d+x2*(a23.d+x2*(a25.d+
diff --git a/sysdeps/ieee754/dbl-64/sincos32.c b/sysdeps/ieee754/dbl-64/sincos32.c
index 6e8c9d5530..0fee643a54 100644
--- a/sysdeps/ieee754/dbl-64/sincos32.c
+++ b/sysdeps/ieee754/dbl-64/sincos32.c
@@ -61,17 +61,17 @@ static void ss32(mp_no *x, mp_no *y, int p) {
#endif
for (i=1;i<=p;i++) mpk.d[i]=0;
- mul(x,x,&x2,p);
- cpy(&oofac27,&gor,p);
- cpy(&gor,&sum,p);
+ __mul(x,x,&x2,p);
+ __cpy(&oofac27,&gor,p);
+ __cpy(&gor,&sum,p);
for (a=27.0;a>1.0;a-=2.0) {
mpk.d[1]=a*(a-1.0);
- mul(&gor,&mpk,&mpt1,p);
- cpy(&mpt1,&gor,p);
- mul(&x2,&sum,&mpt1,p);
- sub(&gor,&mpt1,&sum,p);
+ __mul(&gor,&mpk,&mpt1,p);
+ __cpy(&mpt1,&gor,p);
+ __mul(&x2,&sum,&mpt1,p);
+ __sub(&gor,&mpt1,&sum,p);
}
- mul(x,&sum,y,p);
+ __mul(x,&sum,y,p);
}
/**********************************************************************/
@@ -91,18 +91,18 @@ static void cc32(mp_no *x, mp_no *y, int p) {
#endif
for (i=1;i<=p;i++) mpk.d[i]=0;
- mul(x,x,&x2,p);
+ __mul(x,x,&x2,p);
mpk.d[1]=27.0;
- mul(&oofac27,&mpk,&gor,p);
- cpy(&gor,&sum,p);
+ __mul(&oofac27,&mpk,&gor,p);
+ __cpy(&gor,&sum,p);
for (a=26.0;a>2.0;a-=2.0) {
mpk.d[1]=a*(a-1.0);
- mul(&gor,&mpk,&mpt1,p);
- cpy(&mpt1,&gor,p);
- mul(&x2,&sum,&mpt1,p);
- sub(&gor,&mpt1,&sum,p);
+ __mul(&gor,&mpk,&mpt1,p);
+ __cpy(&mpt1,&gor,p);
+ __mul(&x2,&sum,&mpt1,p);
+ __sub(&gor,&mpt1,&sum,p);
}
- mul(&x2,&sum,y,p);
+ __mul(&x2,&sum,y,p);
}
/***************************************************************************/
@@ -112,20 +112,20 @@ void __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}};
mp_no u,t,t1,t2,c,s;
int i;
- cpy(x,&u,p);
+ __cpy(x,&u,p);
u.e=u.e-1;
cc32(&u,&c,p);
ss32(&u,&s,p);
for (i=0;i<24;i++) {
- mul(&c,&s,&t,p);
- sub(&s,&t,&t1,p);
- add(&t1,&t1,&s,p);
- sub(&mpt,&c,&t1,p);
- mul(&t1,&c,&t2,p);
- add(&t2,&t2,&c,p);
+ __mul(&c,&s,&t,p);
+ __sub(&s,&t,&t1,p);
+ __add(&t1,&t1,&s,p);
+ __sub(&mpt,&c,&t1,p);
+ __mul(&t1,&c,&t2,p);
+ __add(&t2,&t2,&c,p);
}
- sub(&one,&c,y,p);
- cpy(&s,z,p);
+ __sub(&one,&c,y,p);
+ __cpy(&s,z,p);
}
/************************************************************************/
@@ -137,16 +137,16 @@ double __sin32(double x, double res, double res1) {
int p;
mp_no a,b,c;
p=32;
- dbl_mp(res,&a,p);
- dbl_mp(0.5*(res1-res),&b,p);
- add(&a,&b,&c,p);
+ __dbl_mp(res,&a,p);
+ __dbl_mp(0.5*(res1-res),&b,p);
+ __add(&a,&b,&c,p);
if (x>0.8)
- { sub(&hp,&c,&a,p);
+ { __sub(&hp,&c,&a,p);
__c32(&a,&b,&c,p);
}
else __c32(&c,&a,&b,p); /* b=sin(0.5*(res+res1)) */
- dbl_mp(x,&c,p); /* c = x */
- sub(&b,&c,&a,p);
+ __dbl_mp(x,&c,p); /* c = x */
+ __sub(&b,&c,&a,p);
/* if a>0 return min(res,res1), otherwise return max(res,res1) */
if (a.d[0]>0) return (res<res1)?res:res1;
else return (res>res1)?res:res1;
@@ -161,21 +161,21 @@ double __cos32(double x, double res, double res1) {
int p;
mp_no a,b,c;
p=32;
- dbl_mp(res,&a,p);
- dbl_mp(0.5*(res1-res),&b,p);
- add(&a,&b,&c,p);
+ __dbl_mp(res,&a,p);
+ __dbl_mp(0.5*(res1-res),&b,p);
+ __add(&a,&b,&c,p);
if (x>2.4)
- { sub(&pi,&c,&a,p);
+ { __sub(&pi,&c,&a,p);
__c32(&a,&b,&c,p);
b.d[0]=-b.d[0];
}
else if (x>0.8)
- { sub(&hp,&c,&a,p);
+ { __sub(&hp,&c,&a,p);
__c32(&a,&c,&b,p);
}
else __c32(&c,&b,&a,p); /* b=cos(0.5*(res+res1)) */
- dbl_mp(x,&c,p); /* c = x */
- sub(&b,&c,&a,p);
+ __dbl_mp(x,&c,p); /* c = x */
+ __sub(&b,&c,&a,p);
/* if a>0 return max(res,res1), otherwise return min(res,res1) */
if (a.d[0]>0) return (res>res1)?res:res1;
else return (res<res1)?res:res1;
@@ -190,10 +190,10 @@ double __mpsin(double x, double dx) {
double y;
mp_no a,b,c;
p=32;
- dbl_mp(x,&a,p);
- dbl_mp(dx,&b,p);
- add(&a,&b,&c,p);
- if (x>0.8) { sub(&hp,&c,&a,p); __c32(&a,&b,&c,p); }
+ __dbl_mp(x,&a,p);
+ __dbl_mp(dx,&b,p);
+ __add(&a,&b,&c,p);
+ if (x>0.8) { __sub(&hp,&c,&a,p); __c32(&a,&b,&c,p); }
else __c32(&c,&a,&b,p); /* b = sin(x+dx) */
__mp_dbl(&b,&y,p);
return y;
@@ -208,11 +208,11 @@ double __mpcos(double x, double dx) {
double y;
mp_no a,b,c;
p=32;
- dbl_mp(x,&a,p);
- dbl_mp(dx,&b,p);
- add(&a,&b,&c,p);
+ __dbl_mp(x,&a,p);
+ __dbl_mp(dx,&b,p);
+ __add(&a,&b,&c,p);
if (x>0.8)
- { sub(&hp,&c,&b,p);
+ { __sub(&hp,&c,&b,p);
__c32(&b,&a,&c,p);
}
else __c32(&c,&a,&b,p); /* a = cos(x+dx) */
@@ -226,7 +226,7 @@ double __mpcos(double x, double dx) {
/* n=0,+-1,+-2,.... */
/* Return int which indicates in which quarter of circle x is */
/******************************************************************/
-int mpranred(double x, mp_no *y, int p)
+int __mpranred(double x, mp_no *y, int p)
{
number v;
double t,xn;
@@ -239,31 +239,31 @@ int mpranred(double x, mp_no *y, int p)
xn = t - toint.d;
v.d = t;
n =v.i[LOW_HALF]&3;
- dbl_mp(xn,&a,p);
- mul(&a,&hp,&b,p);
- dbl_mp(x,&c,p);
- sub(&c,&b,y,p);
+ __dbl_mp(xn,&a,p);
+ __mul(&a,&hp,&b,p);
+ __dbl_mp(x,&c,p);
+ __sub(&c,&b,y,p);
return n;
}
else { /* if x is very big more precision required */
- dbl_mp(x,&a,p);
+ __dbl_mp(x,&a,p);
a.d[0]=1.0;
k = a.e-5;
if (k < 0) k=0;
b.e = -k;
b.d[0] = 1.0;
for (i=0;i<p;i++) b.d[i+1] = toverp[i+k];
- mul(&a,&b,&c,p);
+ __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)
{ t +=1.0;
- sub(&c,&one,&b,p);
- mul(&b,&hp,y,p);
+ __sub(&c,&one,&b,p);
+ __mul(&b,&hp,y,p);
}
- else mul(&c,&hp,y,p);
+ else __mul(&c,&hp,y,p);
n = (int) t;
if (x < 0) { y->d[0] = - y->d[0]; n = -n; }
return (n&3);
@@ -274,14 +274,14 @@ int mpranred(double x, mp_no *y, int p)
/* Multi-Precision sin() function subroutine, for p=32. It is */
/* based on the routines mpranred() and c32(). */
/*******************************************************************/
-double mpsin1(double x)
+double __mpsin1(double x)
{
int p;
int n;
mp_no u,s,c;
double y;
p=32;
- n=mpranred(x,&u,p); /* n is 0, 1, 2 or 3 */
+ n=__mpranred(x,&u,p); /* n is 0, 1, 2 or 3 */
__c32(&u,&c,&s,p);
switch (n) { /* in which quarter of unit circle y is*/
case 0:
@@ -313,7 +313,7 @@ double mpsin1(double x)
/* based on the routines mpranred() and c32(). */
/*****************************************************************/
-double mpcos1(double x)
+double __mpcos1(double x)
{
int p;
int n;
@@ -321,7 +321,7 @@ double mpcos1(double x)
double y;
p=32;
- n=mpranred(x,&u,p); /* n is 0, 1, 2 or 3 */
+ n=__mpranred(x,&u,p); /* n is 0, 1, 2 or 3 */
__c32(&u,&c,&s,p);
switch (n) { /* in what quarter of unit circle y is*/
diff --git a/sysdeps/ieee754/dbl-64/slowexp.c b/sysdeps/ieee754/dbl-64/slowexp.c
index 1c2779b0e9..3a3758b305 100644
--- a/sysdeps/ieee754/dbl-64/slowexp.c
+++ b/sysdeps/ieee754/dbl-64/slowexp.c
@@ -30,10 +30,10 @@
/**************************************************************************/
#include "mpa.h"
-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 slowexp(double x) {
+double __slowexp(double x) {
double w,z,res,eps=3.0e-26;
#if 0
double y;
@@ -45,20 +45,20 @@ double slowexp(double x) {
mp_no mpx, mpy, mpz,mpw,mpeps,mpcor;
p=6;
- dbl_mp(x,&mpx,p); /* Convert a double precision number x */
+ __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);
+ __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);
+ __dbl_mp(x,&mpx,p);
+ __mpexp(&mpx, &mpy, p);
__mp_dbl(&mpy, &res, p);
return res;
}
diff --git a/sysdeps/ieee754/dbl-64/slowpow.c b/sysdeps/ieee754/dbl-64/slowpow.c
index 11da7e43d9..3412197aef 100644
--- a/sysdeps/ieee754/dbl-64/slowpow.c
+++ b/sysdeps/ieee754/dbl-64/slowpow.c
@@ -34,40 +34,40 @@
#include "mpa.h"
-void mpexp(mp_no *x, mp_no *y, int p);
-void mplog(mp_no *x, mp_no *y, int p);
+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 __halfulp(double x,double y);
-double slowpow(double x, double y, double z) {
+double __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 */
+ 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 */
+ __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 */
+ __sub(&mpp,&eps,&mpr1,p); /* pp -eps =r1 */
__mp_dbl(&mpr1, &res1, p); /* converting into double precision */
if (res == res1) return res;
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 */
+ __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 res;
}