Update.

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.

	redefinition.
	* 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.

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;
 }