Merge changes between r15223 and r15532 from /fsf/trunk.

git-svn-id: svn://svn.eglibc.org/trunk@15545 7b3dc134-2b1b-0410-93df-9e9f96275f8d

48 files changed, 1278 insertions, 1266 deletions

diff --git a/libc/sysdeps/ieee754/dbl-64/dla.h b/libc/sysdeps/ieee754/dbl-64/dla.h
index bf73fa902..cb12dbc8f 100644
--- a/libc/sysdeps/ieee754/dbl-64/dla.h
+++ b/libc/sysdeps/ieee754/dbl-64/dla.h
@@ -1,7 +1,7 @@
 /*
  * IBM Accurate Mathematical Library
  * Written by International Business Machines Corp.
- * Copyright (C) 2001 Free Software Foundation, Inc.
+ * Copyright (C) 2001, 2011 Free Software Foundation, Inc.
  *
  * This program is free software; you can redistribute it and/or modify
  * it under the terms of the GNU Lesser General Public License as published by
@@ -44,7 +44,7 @@
 /* z+zz = x+y exactly.                                                 */
 
 #define  EADD(x,y,z,zz)  \
-           z=(x)+(y);  zz=(ABS(x)>ABS(y)) ? (((x)-(z))+(y)) : (((y)-(z))+(x));
+	   z=(x)+(y);  zz=(ABS(x)>ABS(y)) ? (((x)-(z))+(y)) : (((y)-(z))+(x));
 
 
 /* Exact subtraction of two single-length floating point numbers, Dekker. */
@@ -52,7 +52,7 @@
 /* z+zz = x-y exactly.                                                    */
 
 #define  ESUB(x,y,z,zz)  \
-           z=(x)-(y);  zz=(ABS(x)>ABS(y)) ? (((x)-(z))-(y)) : ((x)-((y)+(z)));
+	   z=(x)-(y);  zz=(ABS(x)>ABS(y)) ? (((x)-(z))-(y)) : ((x)-((y)+(z)));
 
 
 /* Exact multiplication of two single-length floating point numbers,   */
@@ -60,10 +60,15 @@
 /* satisfies z+zz = x*y exactly. p,hx,tx,hy,ty are temporary           */
 /* storage variables of type double.                                   */
 
-#define  EMULV(x,y,z,zz,p,hx,tx,hy,ty)          \
-           p=CN*(x);  hx=((x)-p)+p;  tx=(x)-hx; \
-           p=CN*(y);  hy=((y)-p)+p;  ty=(y)-hy; \
-           z=(x)*(y); zz=(((hx*hy-z)+hx*ty)+tx*hy)+tx*ty;
+#ifdef DLA_FMS
+# define  EMULV(x,y,z,zz,p,hx,tx,hy,ty)          \
+	   z=x*y; zz=DLA_FMS(x,y,z);
+#else
+# define  EMULV(x,y,z,zz,p,hx,tx,hy,ty)          \
+	   p=CN*(x);  hx=((x)-p)+p;  tx=(x)-hx; \
+	   p=CN*(y);  hy=((y)-p)+p;  ty=(y)-hy; \
+	   z=(x)*(y); zz=(((hx*hy-z)+hx*ty)+tx*hy)+tx*ty;
+#endif
 
 
 /* Exact multiplication of two single-length floating point numbers, Dekker. */
@@ -71,10 +76,15 @@
 /* that satisfies z+zz = x*y exactly. p,hx,tx,hy,ty,q are temporary          */
 /* storage variables of type double.                                         */
 
-#define  MUL12(x,y,z,zz,p,hx,tx,hy,ty,q)        \
-           p=CN*(x);  hx=((x)-p)+p;  tx=(x)-hx; \
-           p=CN*(y);  hy=((y)-p)+p;  ty=(y)-hy; \
-           p=hx*hy;  q=hx*ty+tx*hy; z=p+q;  zz=((p-z)+q)+tx*ty;
+#ifdef DLA_FMS
+# define  MUL12(x,y,z,zz,p,hx,tx,hy,ty,q)        \
+	   EMULV(x,y,z,zz,p,hx,tx,hy,ty)
+#else
+# define  MUL12(x,y,z,zz,p,hx,tx,hy,ty,q)        \
+	   p=CN*(x);  hx=((x)-p)+p;  tx=(x)-hx; \
+	   p=CN*(y);  hy=((y)-p)+p;  ty=(y)-hy; \
+	   p=hx*hy;  q=hx*ty+tx*hy; z=p+q;  zz=((p-z)+q)+tx*ty;
+#endif
 
 
 /* Double-length addition, Dekker. The macro produces a double-length   */
@@ -84,10 +94,10 @@
 /* storage variables of type double.                                    */
 
 #define  ADD2(x,xx,y,yy,z,zz,r,s)                    \
-           r=(x)+(y);  s=(ABS(x)>ABS(y)) ?           \
-                       (((((x)-r)+(y))+(yy))+(xx)) : \
-                       (((((y)-r)+(x))+(xx))+(yy));  \
-           z=r+s;  zz=(r-z)+s;
+	   r=(x)+(y);  s=(ABS(x)>ABS(y)) ?           \
+		       (((((x)-r)+(y))+(yy))+(xx)) : \
+		       (((((y)-r)+(x))+(xx))+(yy));  \
+	   z=r+s;  zz=(r-z)+s;
 
 
 /* Double-length subtraction, Dekker. The macro produces a double-length  */
@@ -97,10 +107,10 @@
 /* storage variables of type double.                                      */
 
 #define  SUB2(x,xx,y,yy,z,zz,r,s)                    \
-           r=(x)-(y);  s=(ABS(x)>ABS(y)) ?           \
-                       (((((x)-r)-(y))-(yy))+(xx)) : \
-                       ((((x)-((y)+r))+(xx))-(yy));  \
-           z=r+s;  zz=(r-z)+s;
+	   r=(x)-(y);  s=(ABS(x)>ABS(y)) ?           \
+		       (((((x)-r)-(y))-(yy))+(xx)) : \
+		       ((((x)-((y)+r))+(xx))-(yy));  \
+	   z=r+s;  zz=(r-z)+s;
 
 
 /* Double-length multiplication, Dekker. The macro produces a double-length  */
@@ -110,8 +120,8 @@
 /* temporary storage variables of type double.                               */
 
 #define  MUL2(x,xx,y,yy,z,zz,p,hx,tx,hy,ty,q,c,cc)  \
-           MUL12(x,y,c,cc,p,hx,tx,hy,ty,q)          \
-           cc=((x)*(yy)+(xx)*(y))+cc;   z=c+cc;   zz=(c-z)+cc;
+	   MUL12(x,y,c,cc,p,hx,tx,hy,ty,q)          \
+	   cc=((x)*(yy)+(xx)*(y))+cc;   z=c+cc;   zz=(c-z)+cc;
 
 
 /* Double-length division, Dekker. The macro produces a double-length        */
@@ -121,8 +131,8 @@
 /* are temporary storage variables of type double.                           */
 
 #define  DIV2(x,xx,y,yy,z,zz,p,hx,tx,hy,ty,q,c,cc,u,uu)  \
-           c=(x)/(y);   MUL12(c,y,u,uu,p,hx,tx,hy,ty,q)  \
-           cc=(((((x)-u)-uu)+(xx))-c*(yy))/(y);   z=c+cc;   zz=(c-z)+cc;
+	   c=(x)/(y);   MUL12(c,y,u,uu,p,hx,tx,hy,ty,q)  \
+	   cc=(((((x)-u)-uu)+(xx))-c*(yy))/(y);   z=c+cc;   zz=(c-z)+cc;
 
 
 /* Double-length addition, slower but more accurate than ADD2.               */
@@ -133,17 +143,17 @@
 /* are temporary storage variables of type double.                           */
 
 #define  ADD2A(x,xx,y,yy,z,zz,r,rr,s,ss,u,uu,w)                        \
-           r=(x)+(y);                                                  \
-           if (ABS(x)>ABS(y)) { rr=((x)-r)+(y);  s=(rr+(yy))+(xx); }   \
-           else               { rr=((y)-r)+(x);  s=(rr+(xx))+(yy); }   \
-           if (rr!=0.0) {                                              \
-             z=r+s;  zz=(r-z)+s; }                                     \
-           else {                                                      \
-             ss=(ABS(xx)>ABS(yy)) ? (((xx)-s)+(yy)) : (((yy)-s)+(xx)); \
-             u=r+s;                                                    \
-             uu=(ABS(r)>ABS(s))   ? ((r-u)+s)   : ((s-u)+r)  ;         \
-             w=uu+ss;  z=u+w;                                          \
-             zz=(ABS(u)>ABS(w))   ? ((u-z)+w)   : ((w-z)+u)  ; }
+	   r=(x)+(y);                                                  \
+	   if (ABS(x)>ABS(y)) { rr=((x)-r)+(y);  s=(rr+(yy))+(xx); }   \
+	   else               { rr=((y)-r)+(x);  s=(rr+(xx))+(yy); }   \
+	   if (rr!=0.0) {                                              \
+	     z=r+s;  zz=(r-z)+s; }                                     \
+	   else {                                                      \
+	     ss=(ABS(xx)>ABS(yy)) ? (((xx)-s)+(yy)) : (((yy)-s)+(xx)); \
+	     u=r+s;                                                    \
+	     uu=(ABS(r)>ABS(s))   ? ((r-u)+s)   : ((s-u)+r)  ;         \
+	     w=uu+ss;  z=u+w;                                          \
+	     zz=(ABS(u)>ABS(w))   ? ((u-z)+w)   : ((w-z)+u)  ; }
 
 
 /* Double-length subtraction, slower but more accurate than SUB2.            */
@@ -154,21 +164,14 @@
 /* are temporary storage variables of type double.                           */
 
 #define  SUB2A(x,xx,y,yy,z,zz,r,rr,s,ss,u,uu,w)                        \
-           r=(x)-(y);                                                  \
-           if (ABS(x)>ABS(y)) { rr=((x)-r)-(y);  s=(rr-(yy))+(xx); }   \
-           else               { rr=(x)-((y)+r);  s=(rr+(xx))-(yy); }   \
-           if (rr!=0.0) {                                              \
-             z=r+s;  zz=(r-z)+s; }                                     \
-           else {                                                      \
-             ss=(ABS(xx)>ABS(yy)) ? (((xx)-s)-(yy)) : ((xx)-((yy)+s)); \
-             u=r+s;                                                    \
-             uu=(ABS(r)>ABS(s))   ? ((r-u)+s)   : ((s-u)+r)  ;         \
-             w=uu+ss;  z=u+w;                                          \
-             zz=(ABS(u)>ABS(w))   ? ((u-z)+w)   : ((w-z)+u)  ; }
-
-
-
-
-
-
-
+	   r=(x)-(y);                                                  \
+	   if (ABS(x)>ABS(y)) { rr=((x)-r)-(y);  s=(rr-(yy))+(xx); }   \
+	   else               { rr=(x)-((y)+r);  s=(rr+(xx))-(yy); }   \
+	   if (rr!=0.0) {                                              \
+	     z=r+s;  zz=(r-z)+s; }                                     \
+	   else {                                                      \
+	     ss=(ABS(xx)>ABS(yy)) ? (((xx)-s)-(yy)) : ((xx)-((yy)+s)); \
+	     u=r+s;                                                    \
+	     uu=(ABS(r)>ABS(s))   ? ((r-u)+s)   : ((s-u)+r)  ;         \
+	     w=uu+ss;  z=u+w;                                          \
+	     zz=(ABS(u)>ABS(w))   ? ((u-z)+w)   : ((w-z)+u)  ; }
diff --git a/libc/sysdeps/ieee754/dbl-64/doasin.c b/libc/sysdeps/ieee754/dbl-64/doasin.c
index 79f344af2..c21d4b7df 100644
--- a/libc/sysdeps/ieee754/dbl-64/doasin.c
+++ b/libc/sysdeps/ieee754/dbl-64/doasin.c
@@ -1,7 +1,7 @@
 /*
  * IBM Accurate Mathematical Library
  * written by International Business Machines Corp.
- * Copyright (C) 2001 Free Software Foundation
+ * Copyright (C) 2001, 2011 Free Software Foundation
  *
  * This program is free software; you can redistribute it and/or modify
  * it under the terms of the GNU Lesser General Public License as published by
@@ -31,7 +31,7 @@
 
 #include "endian.h"
 #include "mydefs.h"
-#include "dla.h"
+#include <dla.h>
 #include "math_private.h"
 
 /********************************************************************/
@@ -52,7 +52,10 @@ void __doasin(double x, double dx, double v[]) {
     d11 = 0.79470250400727425881446981833568758E-02;
 
   double xx,p,pp,u,uu,r,s;
-  double hx,tx,hy,ty,tp,tq,tc,tcc;
+  double tc,tcc;
+#ifndef DLA_FMA
+  double hx,tx,hy,ty,tp,tq;
+#endif
 
 
 /* Taylor series for arcsin for Double-Length numbers         */
diff --git a/libc/sysdeps/ieee754/dbl-64/dosincos.c b/libc/sysdeps/ieee754/dbl-64/dosincos.c
index 1d347a4bc..4ae88c31c 100644
--- a/libc/sysdeps/ieee754/dbl-64/dosincos.c
+++ b/libc/sysdeps/ieee754/dbl-64/dosincos.c
@@ -1,7 +1,7 @@
 /*
  * IBM Accurate Mathematical Library
  * written by International Business Machines Corp.
- * Copyright (C) 2001 Free Software Foundation
+ * Copyright (C) 2001, 2011 Free Software Foundation
  *
  * This program is free software; you can redistribute it and/or modify
  * it under the terms of the GNU Lesser General Public License as published by
@@ -36,7 +36,7 @@
 #include "endian.h"
 #include "mydefs.h"
 #include "sincos.tbl"
-#include "dla.h"
+#include <dla.h>
 #include "dosincos.h"
 #include "math_private.h"
 
@@ -48,8 +48,11 @@
 /***********************************************************************/
 
 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,
+  double r,s,c,cc,d,dd,d2,dd2,e,ee,
     sn,ssn,cs,ccs,ds,dss,dc,dcc;
+#ifndef DLA_FMA
+  double p,hx,tx,hy,ty,q;
+#endif
 #if 0
   double xx,y,yy,z,zz;
 #endif
@@ -61,7 +64,7 @@ void __dubsin(double x, double dx, double v[]) {
   x=x-(u.x-big.x);
   d=x+dx;
   dd=(x-d)+dx;
-         /* sin(x+dx)=sin(Xi+t)=sin(Xi)*cos(t) + cos(Xi)sin(t) where t ->0 */
+	 /* sin(x+dx)=sin(Xi+t)=sin(Xi)*cos(t) + cos(Xi)sin(t) where t ->0 */
   MUL2(d,dd,d,dd,d2,dd2,p,hx,tx,hy,ty,q,c,cc);
   sn=sincos.x[k];     /*                                  */
   ssn=sincos.x[k+1];  /*      sin(Xi) and cos(Xi)         */
@@ -99,8 +102,11 @@ void __dubsin(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,
+  double r,s,c,cc,d,dd,d2,dd2,e,ee,
     sn,ssn,cs,ccs,ds,dss,dc,dcc;
+#ifndef DLA_FMA
+  double p,hx,tx,hy,ty,q;
+#endif
 #if 0
   double xx,y,yy,z,zz;
 #endif
@@ -166,15 +172,15 @@ void __docos(double x, double dx, double v[]) {
   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];}
-                                       /* cos(x) = sin ( 90 -  x ) */
-         else {__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 { /* y>= 3/4 * PI */
diff --git a/libc/sysdeps/ieee754/dbl-64/e_acosh.c b/libc/sysdeps/ieee754/dbl-64/e_acosh.c
index 27c29cd8c..6ef10cb84 100644
--- a/libc/sysdeps/ieee754/dbl-64/e_acosh.c
+++ b/libc/sysdeps/ieee754/dbl-64/e_acosh.c
@@ -5,18 +5,14 @@
  *
  * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  */
 
-#if defined(LIBM_SCCS) && !defined(lint)
-static char rcsid[] = "$NetBSD: e_acosh.c,v 1.9 1995/05/12 04:57:18 jtc Exp $";
-#endif
-
 /* __ieee754_acosh(x)
  * Method :
- *	Based on 
+ *	Based on
  *		acosh(x) = log [ x + sqrt(x*x-1) ]
  *	we have
  *		acosh(x) := log(x)+ln2,	if x is large; else
@@ -31,21 +27,13 @@ static char rcsid[] = "$NetBSD: e_acosh.c,v 1.9 1995/05/12 04:57:18 jtc Exp $";
 #include "math.h"
 #include "math_private.h"
 
-#ifdef __STDC__
-static const double 
-#else
-static double 
-#endif
+static const double
 one	= 1.0,
 ln2	= 6.93147180559945286227e-01;  /* 0x3FE62E42, 0xFEFA39EF */
 
-#ifdef __STDC__
-	double __ieee754_acosh(double x)
-#else
-	double __ieee754_acosh(x)
-	double x;
-#endif
-{	
+double
+__ieee754_acosh(double x)
+{
 	double t;
 	int32_t hx;
 	u_int32_t lx;
@@ -54,8 +42,8 @@ ln2	= 6.93147180559945286227e-01;  /* 0x3FE62E42, 0xFEFA39EF */
 	    return (x-x)/(x-x);
 	} else if(hx >=0x41b00000) {	/* x > 2**28 */
 	    if(hx >=0x7ff00000) {	/* x is inf of NaN */
-	        return x+x;
-	    } else 
+		return x+x;
+	    } else
 		return __ieee754_log(x)+ln2;	/* acosh(huge)=log(2x) */
 	} else if(((hx-0x3ff00000)|lx)==0) {
 	    return 0.0;			/* acosh(1) = 0 */
@@ -64,6 +52,7 @@ ln2	= 6.93147180559945286227e-01;  /* 0x3FE62E42, 0xFEFA39EF */
 	    return __ieee754_log(2.0*x-one/(x+__ieee754_sqrt(t-one)));
 	} else {			/* 1<x<2 */
 	    t = x-one;
-	    return __log1p(t+__sqrt(2.0*t+t*t));
+	    return __log1p(t+__ieee754_sqrt(2.0*t+t*t));
 	}
 }
+strong_alias (__ieee754_acosh, __acosh_finite)
diff --git a/libc/sysdeps/ieee754/dbl-64/e_asin.c b/libc/sysdeps/ieee754/dbl-64/e_asin.c
index ce5d227b7..02efb7ad2 100644
--- a/libc/sysdeps/ieee754/dbl-64/e_asin.c
+++ b/libc/sysdeps/ieee754/dbl-64/e_asin.c
@@ -1,7 +1,7 @@
 /*
  * IBM Accurate Mathematical Library
  * written by International Business Machines Corp.
- * Copyright (C) 2001 Free Software Foundation
+ * Copyright (C) 2001, 2011 Free Software Foundation
  *
  * This program is free software; you can redistribute it and/or modify
  * it under the terms of the GNU Lesser General Public License as published by
@@ -209,9 +209,9 @@ double __ieee754_asin(double x){
     else xx = -x - asncs.x[n];
     t = asncs.x[n+1]*xx;
     p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+
-                      xx*(asncs.x[n+5]+xx*(asncs.x[n+6]
+		      xx*(asncs.x[n+5]+xx*(asncs.x[n+6]
 		      +xx*(asncs.x[n+7]+xx*(asncs.x[n+8]+
-                      xx*asncs.x[n+9])))))))+asncs.x[n+10];
+		      xx*asncs.x[n+9])))))))+asncs.x[n+10];
     t+=p;
     res =asncs.x[n+11] +t;
     cor = (asncs.x[n+11]-res)+t;
@@ -248,9 +248,9 @@ double __ieee754_asin(double x){
     else xx = -x - asncs.x[n];
     t = asncs.x[n+1]*xx;
     p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+
-                         xx*(asncs.x[n+5]+xx*(asncs.x[n+6]
+			 xx*(asncs.x[n+5]+xx*(asncs.x[n+6]
 			 +xx*(asncs.x[n+7]+xx*(asncs.x[n+8]+
-                    xx*(asncs.x[n+9]+xx*asncs.x[n+10]))))))))+asncs.x[n+11];
+		    xx*(asncs.x[n+9]+xx*asncs.x[n+10]))))))))+asncs.x[n+11];
     t+=p;
     res =asncs.x[n+12] +t;
     cor = (asncs.x[n+12]-res)+t;
@@ -324,6 +324,7 @@ double __ieee754_asin(double x){
     return u.x/v.x;  /* NaN */
  }
 }
+strong_alias (__ieee754_asin, __asin_finite)
 
 /*******************************************************************/
 /*                                                                 */
@@ -397,7 +398,7 @@ double __ieee754_acos(double x)
     else xx = -x - asncs.x[n];
     t = asncs.x[n+1]*xx;
     p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+
-                   xx*(asncs.x[n+5]+xx*asncs.x[n+6]))))+asncs.x[n+7];
+		   xx*(asncs.x[n+5]+xx*asncs.x[n+6]))))+asncs.x[n+7];
     t+=p;
     y = (m>0)?(hp0.x-asncs.x[n+8]):(hp0.x+asncs.x[n+8]);
     t = (m>0)?(hp1.x-t):(hp1.x+t);
@@ -433,8 +434,8 @@ double __ieee754_acos(double x)
     else {xx = -x - asncs.x[n]; eps=1.02; }
     t = asncs.x[n+1]*xx;
     p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+
-                   xx*(asncs.x[n+5]+xx*(asncs.x[n+6]+
-                   xx*asncs.x[n+7])))))+asncs.x[n+8];
+		   xx*(asncs.x[n+5]+xx*(asncs.x[n+6]+
+		   xx*asncs.x[n+7])))))+asncs.x[n+8];
     t+=p;
    y = (m>0)?(hp0.x-asncs.x[n+9]):(hp0.x+asncs.x[n+9]);
    t = (m>0)?(hp1.x-t):(hp1.x+t);
@@ -468,8 +469,8 @@ double __ieee754_acos(double x)
     else {xx = -x - asncs.x[n]; eps = 1.01; }
     t = asncs.x[n+1]*xx;
     p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+
-                      xx*(asncs.x[n+5]+xx*(asncs.x[n+6]+xx*(asncs.x[n+7]+
-                      xx*asncs.x[n+8]))))))+asncs.x[n+9];
+		      xx*(asncs.x[n+5]+xx*(asncs.x[n+6]+xx*(asncs.x[n+7]+
+		      xx*asncs.x[n+8]))))))+asncs.x[n+9];
     t+=p;
     y = (m>0)?(hp0.x-asncs.x[n+10]):(hp0.x+asncs.x[n+10]);
     t = (m>0)?(hp1.x-t):(hp1.x+t);
@@ -503,9 +504,9 @@ double __ieee754_acos(double x)
     else {xx = -x - asncs.x[n]; eps =1.005; }
     t = asncs.x[n+1]*xx;
     p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+
-                   xx*(asncs.x[n+5]+xx*(asncs.x[n+6]
+		   xx*(asncs.x[n+5]+xx*(asncs.x[n+6]
 		   +xx*(asncs.x[n+7]+xx*(asncs.x[n+8]+
-                   xx*asncs.x[n+9])))))))+asncs.x[n+10];
+		   xx*asncs.x[n+9])))))))+asncs.x[n+10];
     t+=p;
     y = (m>0)?(hp0.x-asncs.x[n+11]):(hp0.x+asncs.x[n+11]);
     t = (m>0)?(hp1.x-t):(hp1.x+t);
@@ -539,9 +540,9 @@ double __ieee754_acos(double x)
     else {xx = -x - asncs.x[n]; eps=1.005;}
     t = asncs.x[n+1]*xx;
     p=xx*xx*(asncs.x[n+2]+xx*(asncs.x[n+3]+xx*(asncs.x[n+4]+
-            xx*(asncs.x[n+5]+xx*(asncs.x[n+6]
+	    xx*(asncs.x[n+5]+xx*(asncs.x[n+6]
 	    +xx*(asncs.x[n+7]+xx*(asncs.x[n+8]+xx*(asncs.x[n+9]+
-            xx*asncs.x[n+10]))))))))+asncs.x[n+11];
+	    xx*asncs.x[n+10]))))))))+asncs.x[n+11];
     t+=p;
     y = (m>0)?(hp0.x-asncs.x[n+12]):(hp0.x+asncs.x[n+12]);
    t = (m>0)?(hp1.x-t):(hp1.x+t);
@@ -635,3 +636,4 @@ double __ieee754_acos(double x)
     return u.x/v.x;
   }
 }
+strong_alias (__ieee754_acos, __acos_finite)
diff --git a/libc/sysdeps/ieee754/dbl-64/e_atan2.c b/libc/sysdeps/ieee754/dbl-64/e_atan2.c
index 9e1a794ec..f8f678bc5 100644
--- a/libc/sysdeps/ieee754/dbl-64/e_atan2.c
+++ b/libc/sysdeps/ieee754/dbl-64/e_atan2.c
@@ -1,7 +1,7 @@
 /*
  * IBM Accurate Mathematical Library
  * written by International Business Machines Corp.
- * Copyright (C) 2001 Free Software Foundation
+ * Copyright (C) 2001, 2011 Free Software Foundation
  *
  * This program is free software; you can redistribute it and/or modify
  * it under the terms of the GNU Lesser General Public License as published by
@@ -37,7 +37,7 @@
 /*                                                                      */
 /************************************************************************/
 
-#include "dla.h"
+#include <dla.h>
 #include "mpa.h"
 #include "MathLib.h"
 #include "uatan.tbl"
@@ -62,8 +62,11 @@ double __ieee754_atan2(double y,double x) {
   int p;
 #endif
   static const int pr[MM]={6,8,10,20,32};
-  double ax,ay,u,du,u9,ua,v,vv,dv,t1,t2,t3,t4,t5,t6,t7,t8,
-         z,zz,cor,s1,ss1,s2,ss2;
+  double ax,ay,u,du,u9,ua,v,vv,dv,t1,t2,t3,t7,t8,
+	 z,zz,cor,s1,ss1,s2,ss2;
+#ifndef DLA_FMA
+  double t4,t5,t6;
+#endif
 #if 0
   double z1,z2;
 #endif
@@ -73,7 +76,7 @@ double __ieee754_atan2(double y,double x) {
 #endif
 
   static const int ep= 59768832,   /*  57*16**5   */
-                   em=-59768832;   /* -57*16**5   */
+		   em=-59768832;   /* -57*16**5   */
 
   /* x=NaN or y=NaN */
   num.d = x;  ux = num.i[HIGH_HALF];  dx = num.i[LOW_HALF];
@@ -102,23 +105,23 @@ double __ieee754_atan2(double y,double x) {
   if          (ux==0x7ff00000) {
     if        (dx==0x00000000) {
       if      (uy==0x7ff00000) {
-        if    (dy==0x00000000)  return qpi.d; }
+	if    (dy==0x00000000)  return qpi.d; }
       else if (uy==0xfff00000) {
-        if    (dy==0x00000000)  return mqpi.d; }
+	if    (dy==0x00000000)  return mqpi.d; }
       else {
-        if    ((uy&0x80000000)==0x00000000)  return ZERO;
-        else                                 return MZERO; }
+	if    ((uy&0x80000000)==0x00000000)  return ZERO;
+	else                                 return MZERO; }
     }
   }
   else if     (ux==0xfff00000) {
     if        (dx==0x00000000) {
       if      (uy==0x7ff00000) {
-        if    (dy==0x00000000)  return tqpi.d; }
+	if    (dy==0x00000000)  return tqpi.d; }
       else if (uy==0xfff00000) {
-        if    (dy==0x00000000)  return mtqpi.d; }
+	if    (dy==0x00000000)  return mtqpi.d; }
       else                     {
-        if    ((uy&0x80000000)==0x00000000)  return opi.d;
-        else                                 return mopi.d; }
+	if    ((uy&0x80000000)==0x00000000)  return opi.d;
+	else                                 return mopi.d; }
     }
   }
 
@@ -156,108 +159,108 @@ double __ieee754_atan2(double y,double x) {
     /* (i)   x>0, abs(y)< abs(x):  atan(ay/ax) */
     if (ay<ax) {
       if (u<inv16.d) {
-        v=u*u;  zz=du+u*v*(d3.d+v*(d5.d+v*(d7.d+v*(d9.d+v*(d11.d+v*d13.d)))));
-        if ((z=u+(zz-u1.d*u)) == u+(zz+u1.d*u))  return signArctan2(y,z);
+	v=u*u;  zz=du+u*v*(d3.d+v*(d5.d+v*(d7.d+v*(d9.d+v*(d11.d+v*d13.d)))));
+	if ((z=u+(zz-u1.d*u)) == u+(zz+u1.d*u))  return signArctan2(y,z);
 
-        MUL2(u,du,u,du,v,vv,t1,t2,t3,t4,t5,t6,t7,t8)
-        s1=v*(f11.d+v*(f13.d+v*(f15.d+v*(f17.d+v*f19.d))));
-        ADD2(f9.d,ff9.d,s1,ZERO,s2,ss2,t1,t2)
-        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
-        ADD2(f7.d,ff7.d,s1,ss1,s2,ss2,t1,t2)
-        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
-        ADD2(f5.d,ff5.d,s1,ss1,s2,ss2,t1,t2)
-        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
-        ADD2(f3.d,ff3.d,s1,ss1,s2,ss2,t1,t2)
-        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
-        MUL2(u,du,s1,ss1,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
-        ADD2(u,du,s2,ss2,s1,ss1,t1,t2)
-        if ((z=s1+(ss1-u5.d*s1)) == s1+(ss1+u5.d*s1))  return signArctan2(y,z);
-        return atan2Mp(x,y,pr);
+	MUL2(u,du,u,du,v,vv,t1,t2,t3,t4,t5,t6,t7,t8)
+	s1=v*(f11.d+v*(f13.d+v*(f15.d+v*(f17.d+v*f19.d))));
+	ADD2(f9.d,ff9.d,s1,ZERO,s2,ss2,t1,t2)
+	MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
+	ADD2(f7.d,ff7.d,s1,ss1,s2,ss2,t1,t2)
+	MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
+	ADD2(f5.d,ff5.d,s1,ss1,s2,ss2,t1,t2)
+	MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
+	ADD2(f3.d,ff3.d,s1,ss1,s2,ss2,t1,t2)
+	MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
+	MUL2(u,du,s1,ss1,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
+	ADD2(u,du,s2,ss2,s1,ss1,t1,t2)
+	if ((z=s1+(ss1-u5.d*s1)) == s1+(ss1+u5.d*s1))  return signArctan2(y,z);
+	return atan2Mp(x,y,pr);
       }
       else {
-        i=(TWO52+TWO8*u)-TWO52;  i-=16;
-        t3=u-cij[i][0].d;
-        EADD(t3,du,v,dv)
-        t1=cij[i][1].d;  t2=cij[i][2].d;
-        zz=v*t2+(dv*t2+v*v*(cij[i][3].d+v*(cij[i][4].d+
-                         v*(cij[i][5].d+v* cij[i][6].d))));
-        if (i<112) {
-          if (i<48)  u9=u91.d;    /* u < 1/4        */
-          else       u9=u92.d; }  /* 1/4 <= u < 1/2 */
-        else {
-          if (i<176) u9=u93.d;    /* 1/2 <= u < 3/4 */
-          else       u9=u94.d; }  /* 3/4 <= u <= 1  */
-        if ((z=t1+(zz-u9*t1)) == t1+(zz+u9*t1))  return signArctan2(y,z);
+	i=(TWO52+TWO8*u)-TWO52;  i-=16;
+	t3=u-cij[i][0].d;
+	EADD(t3,du,v,dv)
+	t1=cij[i][1].d;  t2=cij[i][2].d;
+	zz=v*t2+(dv*t2+v*v*(cij[i][3].d+v*(cij[i][4].d+
+			 v*(cij[i][5].d+v* cij[i][6].d))));
+	if (i<112) {
+	  if (i<48)  u9=u91.d;    /* u < 1/4        */
+	  else       u9=u92.d; }  /* 1/4 <= u < 1/2 */
+	else {
+	  if (i<176) u9=u93.d;    /* 1/2 <= u < 3/4 */
+	  else       u9=u94.d; }  /* 3/4 <= u <= 1  */
+	if ((z=t1+(zz-u9*t1)) == t1+(zz+u9*t1))  return signArctan2(y,z);
 
-        t1=u-hij[i][0].d;
-        EADD(t1,du,v,vv)
-        s1=v*(hij[i][11].d+v*(hij[i][12].d+v*(hij[i][13].d+
-           v*(hij[i][14].d+v* hij[i][15].d))));
-        ADD2(hij[i][9].d,hij[i][10].d,s1,ZERO,s2,ss2,t1,t2)
-        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
-        ADD2(hij[i][7].d,hij[i][8].d,s1,ss1,s2,ss2,t1,t2)
-        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
-        ADD2(hij[i][5].d,hij[i][6].d,s1,ss1,s2,ss2,t1,t2)
-        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
-        ADD2(hij[i][3].d,hij[i][4].d,s1,ss1,s2,ss2,t1,t2)
-        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
-        ADD2(hij[i][1].d,hij[i][2].d,s1,ss1,s2,ss2,t1,t2)
-        if ((z=s2+(ss2-ub.d*s2)) == s2+(ss2+ub.d*s2))  return signArctan2(y,z);
-        return atan2Mp(x,y,pr);
+	t1=u-hij[i][0].d;
+	EADD(t1,du,v,vv)
+	s1=v*(hij[i][11].d+v*(hij[i][12].d+v*(hij[i][13].d+
+	   v*(hij[i][14].d+v* hij[i][15].d))));
+	ADD2(hij[i][9].d,hij[i][10].d,s1,ZERO,s2,ss2,t1,t2)
+	MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
+	ADD2(hij[i][7].d,hij[i][8].d,s1,ss1,s2,ss2,t1,t2)
+	MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
+	ADD2(hij[i][5].d,hij[i][6].d,s1,ss1,s2,ss2,t1,t2)
+	MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
+	ADD2(hij[i][3].d,hij[i][4].d,s1,ss1,s2,ss2,t1,t2)
+	MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
+	ADD2(hij[i][1].d,hij[i][2].d,s1,ss1,s2,ss2,t1,t2)
+	if ((z=s2+(ss2-ub.d*s2)) == s2+(ss2+ub.d*s2))  return signArctan2(y,z);
+	return atan2Mp(x,y,pr);
       }
     }
 
     /* (ii)  x>0, abs(x)<=abs(y):  pi/2-atan(ax/ay) */
     else {
       if (u<inv16.d) {
-        v=u*u;
-        zz=u*v*(d3.d+v*(d5.d+v*(d7.d+v*(d9.d+v*(d11.d+v*d13.d)))));
-        ESUB(hpi.d,u,t2,cor)
-        t3=((hpi1.d+cor)-du)-zz;
-        if ((z=t2+(t3-u2.d)) == t2+(t3+u2.d))  return signArctan2(y,z);
+	v=u*u;
+	zz=u*v*(d3.d+v*(d5.d+v*(d7.d+v*(d9.d+v*(d11.d+v*d13.d)))));
+	ESUB(hpi.d,u,t2,cor)
+	t3=((hpi1.d+cor)-du)-zz;
+	if ((z=t2+(t3-u2.d)) == t2+(t3+u2.d))  return signArctan2(y,z);
 
-        MUL2(u,du,u,du,v,vv,t1,t2,t3,t4,t5,t6,t7,t8)
-        s1=v*(f11.d+v*(f13.d+v*(f15.d+v*(f17.d+v*f19.d))));
-        ADD2(f9.d,ff9.d,s1,ZERO,s2,ss2,t1,t2)
-        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
-        ADD2(f7.d,ff7.d,s1,ss1,s2,ss2,t1,t2)
-        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
-        ADD2(f5.d,ff5.d,s1,ss1,s2,ss2,t1,t2)
-        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
-        ADD2(f3.d,ff3.d,s1,ss1,s2,ss2,t1,t2)
-        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
-        MUL2(u,du,s1,ss1,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
-        ADD2(u,du,s2,ss2,s1,ss1,t1,t2)
-        SUB2(hpi.d,hpi1.d,s1,ss1,s2,ss2,t1,t2)
-        if ((z=s2+(ss2-u6.d)) == s2+(ss2+u6.d))  return signArctan2(y,z);
-        return atan2Mp(x,y,pr);
+	MUL2(u,du,u,du,v,vv,t1,t2,t3,t4,t5,t6,t7,t8)
+	s1=v*(f11.d+v*(f13.d+v*(f15.d+v*(f17.d+v*f19.d))));
+	ADD2(f9.d,ff9.d,s1,ZERO,s2,ss2,t1,t2)
+	MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
+	ADD2(f7.d,ff7.d,s1,ss1,s2,ss2,t1,t2)
+	MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
+	ADD2(f5.d,ff5.d,s1,ss1,s2,ss2,t1,t2)
+	MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
+	ADD2(f3.d,ff3.d,s1,ss1,s2,ss2,t1,t2)
+	MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
+	MUL2(u,du,s1,ss1,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
+	ADD2(u,du,s2,ss2,s1,ss1,t1,t2)
+	SUB2(hpi.d,hpi1.d,s1,ss1,s2,ss2,t1,t2)
+	if ((z=s2+(ss2-u6.d)) == s2+(ss2+u6.d))  return signArctan2(y,z);
+	return atan2Mp(x,y,pr);
       }
       else {
-        i=(TWO52+TWO8*u)-TWO52;  i-=16;
-        v=(u-cij[i][0].d)+du;
-        zz=hpi1.d-v*(cij[i][2].d+v*(cij[i][3].d+v*(cij[i][4].d+
-                                 v*(cij[i][5].d+v* cij[i][6].d))));
-        t1=hpi.d-cij[i][1].d;
-        if (i<112)  ua=ua1.d;  /* w <  1/2 */
-        else        ua=ua2.d;  /* w >= 1/2 */
-        if ((z=t1+(zz-ua)) == t1+(zz+ua))  return signArctan2(y,z);
+	i=(TWO52+TWO8*u)-TWO52;  i-=16;
+	v=(u-cij[i][0].d)+du;
+	zz=hpi1.d-v*(cij[i][2].d+v*(cij[i][3].d+v*(cij[i][4].d+
+				 v*(cij[i][5].d+v* cij[i][6].d))));
+	t1=hpi.d-cij[i][1].d;
+	if (i<112)  ua=ua1.d;  /* w <  1/2 */
+	else        ua=ua2.d;  /* w >= 1/2 */
+	if ((z=t1+(zz-ua)) == t1+(zz+ua))  return signArctan2(y,z);
 
-        t1=u-hij[i][0].d;
-        EADD(t1,du,v,vv)
-        s1=v*(hij[i][11].d+v*(hij[i][12].d+v*(hij[i][13].d+
-           v*(hij[i][14].d+v* hij[i][15].d))));
-        ADD2(hij[i][9].d,hij[i][10].d,s1,ZERO,s2,ss2,t1,t2)
-        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
-        ADD2(hij[i][7].d,hij[i][8].d,s1,ss1,s2,ss2,t1,t2)
-        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
-        ADD2(hij[i][5].d,hij[i][6].d,s1,ss1,s2,ss2,t1,t2)
-        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
-        ADD2(hij[i][3].d,hij[i][4].d,s1,ss1,s2,ss2,t1,t2)
-        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
-        ADD2(hij[i][1].d,hij[i][2].d,s1,ss1,s2,ss2,t1,t2)
-        SUB2(hpi.d,hpi1.d,s2,ss2,s1,ss1,t1,t2)
-        if ((z=s1+(ss1-uc.d)) == s1+(ss1+uc.d))  return signArctan2(y,z);
-        return atan2Mp(x,y,pr);
+	t1=u-hij[i][0].d;
+	EADD(t1,du,v,vv)
+	s1=v*(hij[i][11].d+v*(hij[i][12].d+v*(hij[i][13].d+
+	   v*(hij[i][14].d+v* hij[i][15].d))));
+	ADD2(hij[i][9].d,hij[i][10].d,s1,ZERO,s2,ss2,t1,t2)
+	MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
+	ADD2(hij[i][7].d,hij[i][8].d,s1,ss1,s2,ss2,t1,t2)
+	MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
+	ADD2(hij[i][5].d,hij[i][6].d,s1,ss1,s2,ss2,t1,t2)
+	MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
+	ADD2(hij[i][3].d,hij[i][4].d,s1,ss1,s2,ss2,t1,t2)
+	MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
+	ADD2(hij[i][1].d,hij[i][2].d,s1,ss1,s2,ss2,t1,t2)
+	SUB2(hpi.d,hpi1.d,s2,ss2,s1,ss1,t1,t2)
+	if ((z=s1+(ss1-uc.d)) == s1+(ss1+uc.d))  return signArctan2(y,z);
+	return atan2Mp(x,y,pr);
       }
     }
   }
@@ -266,112 +269,114 @@ double __ieee754_atan2(double y,double x) {
     /* (iii) x<0, abs(x)< abs(y):  pi/2+atan(ax/ay) */
     if (ax<ay) {
       if (u<inv16.d) {
-        v=u*u;
-        zz=u*v*(d3.d+v*(d5.d+v*(d7.d+v*(d9.d+v*(d11.d+v*d13.d)))));
-        EADD(hpi.d,u,t2,cor)
-        t3=((hpi1.d+cor)+du)+zz;
-        if ((z=t2+(t3-u3.d)) == t2+(t3+u3.d))  return signArctan2(y,z);
+	v=u*u;
+	zz=u*v*(d3.d+v*(d5.d+v*(d7.d+v*(d9.d+v*(d11.d+v*d13.d)))));
+	EADD(hpi.d,u,t2,cor)
+	t3=((hpi1.d+cor)+du)+zz;
+	if ((z=t2+(t3-u3.d)) == t2+(t3+u3.d))  return signArctan2(y,z);
 
-        MUL2(u,du,u,du,v,vv,t1,t2,t3,t4,t5,t6,t7,t8)
-        s1=v*(f11.d+v*(f13.d+v*(f15.d+v*(f17.d+v*f19.d))));
-        ADD2(f9.d,ff9.d,s1,ZERO,s2,ss2,t1,t2)
-        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
-        ADD2(f7.d,ff7.d,s1,ss1,s2,ss2,t1,t2)
-        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
-        ADD2(f5.d,ff5.d,s1,ss1,s2,ss2,t1,t2)
-        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
-        ADD2(f3.d,ff3.d,s1,ss1,s2,ss2,t1,t2)
-        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
-        MUL2(u,du,s1,ss1,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
-        ADD2(u,du,s2,ss2,s1,ss1,t1,t2)
-        ADD2(hpi.d,hpi1.d,s1,ss1,s2,ss2,t1,t2)
-        if ((z=s2+(ss2-u7.d)) == s2+(ss2+u7.d))  return signArctan2(y,z);
-        return atan2Mp(x,y,pr);
+	MUL2(u,du,u,du,v,vv,t1,t2,t3,t4,t5,t6,t7,t8)
+	s1=v*(f11.d+v*(f13.d+v*(f15.d+v*(f17.d+v*f19.d))));
+	ADD2(f9.d,ff9.d,s1,ZERO,s2,ss2,t1,t2)
+	MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
+	ADD2(f7.d,ff7.d,s1,ss1,s2,ss2,t1,t2)
+	MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
+	ADD2(f5.d,ff5.d,s1,ss1,s2,ss2,t1,t2)
+	MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
+	ADD2(f3.d,ff3.d,s1,ss1,s2,ss2,t1,t2)
+	MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
+	MUL2(u,du,s1,ss1,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
+	ADD2(u,du,s2,ss2,s1,ss1,t1,t2)
+	ADD2(hpi.d,hpi1.d,s1,ss1,s2,ss2,t1,t2)
+	if ((z=s2+(ss2-u7.d)) == s2+(ss2+u7.d))  return signArctan2(y,z);
+	return atan2Mp(x,y,pr);
       }
       else {
-        i=(TWO52+TWO8*u)-TWO52;  i-=16;
-        v=(u-cij[i][0].d)+du;
-        zz=hpi1.d+v*(cij[i][2].d+v*(cij[i][3].d+v*(cij[i][4].d+
-                                 v*(cij[i][5].d+v* cij[i][6].d))));
-        t1=hpi.d+cij[i][1].d;
-        if (i<112)  ua=ua1.d;  /* w <  1/2 */
-        else        ua=ua2.d;  /* w >= 1/2 */
-        if ((z=t1+(zz-ua)) == t1+(zz+ua))  return signArctan2(y,z);
+	i=(TWO52+TWO8*u)-TWO52;  i-=16;
+	v=(u-cij[i][0].d)+du;
+	zz=hpi1.d+v*(cij[i][2].d+v*(cij[i][3].d+v*(cij[i][4].d+
+				 v*(cij[i][5].d+v* cij[i][6].d))));
+	t1=hpi.d+cij[i][1].d;
+	if (i<112)  ua=ua1.d;  /* w <  1/2 */
+	else        ua=ua2.d;  /* w >= 1/2 */
+	if ((z=t1+(zz-ua)) == t1+(zz+ua))  return signArctan2(y,z);
 
-        t1=u-hij[i][0].d;
-        EADD(t1,du,v,vv)
-        s1=v*(hij[i][11].d+v*(hij[i][12].d+v*(hij[i][13].d+
-           v*(hij[i][14].d+v* hij[i][15].d))));
-        ADD2(hij[i][9].d,hij[i][10].d,s1,ZERO,s2,ss2,t1,t2)
-        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
-        ADD2(hij[i][7].d,hij[i][8].d,s1,ss1,s2,ss2,t1,t2)
-        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
-        ADD2(hij[i][5].d,hij[i][6].d,s1,ss1,s2,ss2,t1,t2)
-        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
-        ADD2(hij[i][3].d,hij[i][4].d,s1,ss1,s2,ss2,t1,t2)
-        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
-        ADD2(hij[i][1].d,hij[i][2].d,s1,ss1,s2,ss2,t1,t2)
-        ADD2(hpi.d,hpi1.d,s2,ss2,s1,ss1,t1,t2)
-        if ((z=s1+(ss1-uc.d)) == s1+(ss1+uc.d))  return signArctan2(y,z);
-        return atan2Mp(x,y,pr);
+	t1=u-hij[i][0].d;
+	EADD(t1,du,v,vv)
+	s1=v*(hij[i][11].d+v*(hij[i][12].d+v*(hij[i][13].d+
+	   v*(hij[i][14].d+v* hij[i][15].d))));
+	ADD2(hij[i][9].d,hij[i][10].d,s1,ZERO,s2,ss2,t1,t2)
+	MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
+	ADD2(hij[i][7].d,hij[i][8].d,s1,ss1,s2,ss2,t1,t2)
+	MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
+	ADD2(hij[i][5].d,hij[i][6].d,s1,ss1,s2,ss2,t1,t2)
+	MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
+	ADD2(hij[i][3].d,hij[i][4].d,s1,ss1,s2,ss2,t1,t2)
+	MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
+	ADD2(hij[i][1].d,hij[i][2].d,s1,ss1,s2,ss2,t1,t2)
+	ADD2(hpi.d,hpi1.d,s2,ss2,s1,ss1,t1,t2)
+	if ((z=s1+(ss1-uc.d)) == s1+(ss1+uc.d))  return signArctan2(y,z);
+	return atan2Mp(x,y,pr);
       }
     }
 
     /* (iv)  x<0, abs(y)<=abs(x):  pi-atan(ax/ay) */
     else {
       if (u<inv16.d) {
-        v=u*u;
-        zz=u*v*(d3.d+v*(d5.d+v*(d7.d+v*(d9.d+v*(d11.d+v*d13.d)))));
-        ESUB(opi.d,u,t2,cor)
-        t3=((opi1.d+cor)-du)-zz;
-        if ((z=t2+(t3-u4.d)) == t2+(t3+u4.d))  return signArctan2(y,z);
+	v=u*u;
+	zz=u*v*(d3.d+v*(d5.d+v*(d7.d+v*(d9.d+v*(d11.d+v*d13.d)))));
+	ESUB(opi.d,u,t2,cor)
+	t3=((opi1.d+cor)-du)-zz;
+	if ((z=t2+(t3-u4.d)) == t2+(t3+u4.d))  return signArctan2(y,z);
 
-        MUL2(u,du,u,du,v,vv,t1,t2,t3,t4,t5,t6,t7,t8)
-        s1=v*(f11.d+v*(f13.d+v*(f15.d+v*(f17.d+v*f19.d))));
-        ADD2(f9.d,ff9.d,s1,ZERO,s2,ss2,t1,t2)
-        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
-        ADD2(f7.d,ff7.d,s1,ss1,s2,ss2,t1,t2)
-        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
-        ADD2(f5.d,ff5.d,s1,ss1,s2,ss2,t1,t2)
-        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
-        ADD2(f3.d,ff3.d,s1,ss1,s2,ss2,t1,t2)
-        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
-        MUL2(u,du,s1,ss1,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
-        ADD2(u,du,s2,ss2,s1,ss1,t1,t2)
-        SUB2(opi.d,opi1.d,s1,ss1,s2,ss2,t1,t2)
-        if ((z=s2+(ss2-u8.d)) == s2+(ss2+u8.d))  return signArctan2(y,z);
-        return atan2Mp(x,y,pr);
+	MUL2(u,du,u,du,v,vv,t1,t2,t3,t4,t5,t6,t7,t8)
+	s1=v*(f11.d+v*(f13.d+v*(f15.d+v*(f17.d+v*f19.d))));
+	ADD2(f9.d,ff9.d,s1,ZERO,s2,ss2,t1,t2)
+	MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
+	ADD2(f7.d,ff7.d,s1,ss1,s2,ss2,t1,t2)
+	MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
+	ADD2(f5.d,ff5.d,s1,ss1,s2,ss2,t1,t2)
+	MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
+	ADD2(f3.d,ff3.d,s1,ss1,s2,ss2,t1,t2)
+	MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
+	MUL2(u,du,s1,ss1,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
+	ADD2(u,du,s2,ss2,s1,ss1,t1,t2)
+	SUB2(opi.d,opi1.d,s1,ss1,s2,ss2,t1,t2)
+	if ((z=s2+(ss2-u8.d)) == s2+(ss2+u8.d))  return signArctan2(y,z);
+	return atan2Mp(x,y,pr);
       }
       else {
-        i=(TWO52+TWO8*u)-TWO52;  i-=16;
-        v=(u-cij[i][0].d)+du;
-        zz=opi1.d-v*(cij[i][2].d+v*(cij[i][3].d+v*(cij[i][4].d+
-                                 v*(cij[i][5].d+v* cij[i][6].d))));
-        t1=opi.d-cij[i][1].d;
-        if (i<112)  ua=ua1.d;  /* w <  1/2 */
-        else        ua=ua2.d;  /* w >= 1/2 */
-        if ((z=t1+(zz-ua)) == t1+(zz+ua))  return signArctan2(y,z);
+	i=(TWO52+TWO8*u)-TWO52;  i-=16;
+	v=(u-cij[i][0].d)+du;
+	zz=opi1.d-v*(cij[i][2].d+v*(cij[i][3].d+v*(cij[i][4].d+
+				 v*(cij[i][5].d+v* cij[i][6].d))));
+	t1=opi.d-cij[i][1].d;
+	if (i<112)  ua=ua1.d;  /* w <  1/2 */
+	else        ua=ua2.d;  /* w >= 1/2 */
+	if ((z=t1+(zz-ua)) == t1+(zz+ua))  return signArctan2(y,z);
 
-        t1=u-hij[i][0].d;
-        EADD(t1,du,v,vv)
-        s1=v*(hij[i][11].d+v*(hij[i][12].d+v*(hij[i][13].d+
-           v*(hij[i][14].d+v* hij[i][15].d))));
-        ADD2(hij[i][9].d,hij[i][10].d,s1,ZERO,s2,ss2,t1,t2)
-        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
-        ADD2(hij[i][7].d,hij[i][8].d,s1,ss1,s2,ss2,t1,t2)
-        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
-        ADD2(hij[i][5].d,hij[i][6].d,s1,ss1,s2,ss2,t1,t2)
-        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
-        ADD2(hij[i][3].d,hij[i][4].d,s1,ss1,s2,ss2,t1,t2)
-        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
-        ADD2(hij[i][1].d,hij[i][2].d,s1,ss1,s2,ss2,t1,t2)
-        SUB2(opi.d,opi1.d,s2,ss2,s1,ss1,t1,t2)
-        if ((z=s1+(ss1-uc.d)) == s1+(ss1+uc.d))  return signArctan2(y,z);
-        return atan2Mp(x,y,pr);
+	t1=u-hij[i][0].d;
+	EADD(t1,du,v,vv)
+	s1=v*(hij[i][11].d+v*(hij[i][12].d+v*(hij[i][13].d+
+	   v*(hij[i][14].d+v* hij[i][15].d))));
+	ADD2(hij[i][9].d,hij[i][10].d,s1,ZERO,s2,ss2,t1,t2)
+	MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
+	ADD2(hij[i][7].d,hij[i][8].d,s1,ss1,s2,ss2,t1,t2)
+	MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
+	ADD2(hij[i][5].d,hij[i][6].d,s1,ss1,s2,ss2,t1,t2)
+	MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
+	ADD2(hij[i][3].d,hij[i][4].d,s1,ss1,s2,ss2,t1,t2)
+	MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
+	ADD2(hij[i][1].d,hij[i][2].d,s1,ss1,s2,ss2,t1,t2)
+	SUB2(opi.d,opi1.d,s2,ss2,s1,ss1,t1,t2)
+	if ((z=s1+(ss1-uc.d)) == s1+(ss1+uc.d))  return signArctan2(y,z);
+	return atan2Mp(x,y,pr);
       }
     }
   }
 }
+strong_alias (__ieee754_atan2, __atan2_finite)
+
   /* Treat the Denormalized case */
 static double  normalized(double ax,double ay,double y, double z)
     { int p;
diff --git a/libc/sysdeps/ieee754/dbl-64/e_atanh.c b/libc/sysdeps/ieee754/dbl-64/e_atanh.c
index fa4fe675c..de3bc8f14 100644
--- a/libc/sysdeps/ieee754/dbl-64/e_atanh.c
+++ b/libc/sysdeps/ieee754/dbl-64/e_atanh.c
@@ -1,74 +1,70 @@
-/* @(#)e_atanh.c 5.1 93/09/24 */
-/*
- * ====================================================
- * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
- *
- * Developed at SunPro, a Sun Microsystems, Inc. business.
- * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
- * is preserved.
- * ====================================================
- */
+/* Copyright (C) 2011 Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
+   Contributed by Ulrich Drepper <drepper@gmail.com>, 2011.
+
+   The GNU C Library 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.1 of the License, or (at your option) any later version.
+
+   The GNU C Library 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 GNU
+   Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public
+   License along with the GNU C Library; if not, write to the Free
+   Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
+   02111-1307 USA.  */
 
-#if defined(LIBM_SCCS) && !defined(lint)
-static char rcsid[] = "$NetBSD: e_atanh.c,v 1.8 1995/05/10 20:44:55 jtc Exp $";
-#endif
 
 /* __ieee754_atanh(x)
- * Method :
- *    1.Reduced x to positive by atanh(-x) = -atanh(x)
- *    2.For x>=0.5
- *                  1              2x                          x
- *	atanh(x) = --- * log(1 + -------) = 0.5 * log1p(2 * --------)
- *                  2             1 - x                      1 - x
- *	
- * 	For x<0.5
- *	atanh(x) = 0.5*log1p(2x+2x*x/(1-x))
- *
- * Special cases:
- *	atanh(x) is NaN if |x| > 1 with signal;
- *	atanh(NaN) is that NaN with no signal;
- *	atanh(+-1) is +-INF with signal.
- *
+   Method :
+      1.Reduced x to positive by atanh(-x) = -atanh(x)
+      2.For x>=0.5
+		    1              2x                          x
+	atanh(x) = --- * log(1 + -------) = 0.5 * log1p(2 * --------)
+		    2             1 - x                      1 - x
+
+	For x<0.5
+	atanh(x) = 0.5*log1p(2x+2x*x/(1-x))
+
+   Special cases:
+	atanh(x) is NaN if |x| > 1 with signal;
+	atanh(NaN) is that NaN with no signal;
+	atanh(+-1) is +-INF with signal.
+
  */
 
+#include <inttypes.h>
 #include "math.h"
 #include "math_private.h"
 
-#ifdef __STDC__
-static const double one = 1.0, huge = 1e300;
-#else
-static double one = 1.0, huge = 1e300;
-#endif
-
-#ifdef __STDC__
-static const double zero = 0.0;
-#else
-static double zero = 0.0;
-#endif
+static const double huge = 1e300;
 
-#ifdef __STDC__
-	double __ieee754_atanh(double x)
-#else
-	double __ieee754_atanh(x)
-	double x;
-#endif
+double
+__ieee754_atanh (double x)
 {
-	double t;
-	int32_t hx,ix;
-	u_int32_t lx;
-	EXTRACT_WORDS(hx,lx,x);
-	ix = hx&0x7fffffff;
-	if ((ix|((lx|(-lx))>>31))>0x3ff00000) /* |x|>1 */
-	    return (x-x)/(x-x);
-	if(ix==0x3ff00000) 
-	    return x/zero;
-	if(ix<0x3e300000&&(huge+x)>zero) return x;	/* x<2**-28 */
-	SET_HIGH_WORD(x,ix);
-	if(ix<0x3fe00000) {		/* x < 0.5 */
-	    t = x+x;
-	    t = 0.5*__log1p(t+t*x/(one-x));
-	} else 
-	    t = 0.5*__log1p((x+x)/(one-x));
-	if(hx>=0) return t; else return -t;
+  double xa = fabs (x);
+  double t;
+  if (xa < 0.5)
+    {
+      if (__builtin_expect (xa < 0x1.0p-28, 0) && (huge + x) > 0.0)
+	return x;
+
+      t = xa + xa;
+      t = 0.5 * __log1p (t + t * xa / (1.0 - xa));
+    }
+  else if (__builtin_expect (xa < 1.0, 1))
+    t = 0.5 * __log1p ((xa + xa) / (1.0 - xa));
+  else
+    {
+      if (xa > 1.0)
+	return (x - x) / (x - x);
+
+      return x / 0.0;
+    }
+
+  return __copysign (t, x);
 }
+strong_alias (__ieee754_atanh, __atanh_finite)
diff --git a/libc/sysdeps/ieee754/dbl-64/e_cosh.c b/libc/sysdeps/ieee754/dbl-64/e_cosh.c
index 65106b998..b9f79e47a 100644
--- a/libc/sysdeps/ieee754/dbl-64/e_cosh.c
+++ b/libc/sysdeps/ieee754/dbl-64/e_cosh.c
@@ -1,32 +1,28 @@
-/* @(#)e_cosh.c 5.1 93/09/24 */
+/* Optimized by Ulrich Drepper <drepper@gmail.com>, 2011 */
 /*
  * ====================================================
  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
  *
  * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  */
 
-#if defined(LIBM_SCCS) && !defined(lint)
-static char rcsid[] = "$NetBSD: e_cosh.c,v 1.7 1995/05/10 20:44:58 jtc Exp $";
-#endif
-
 /* __ieee754_cosh(x)
- * Method : 
+ * Method :
  * mathematically cosh(x) if defined to be (exp(x)+exp(-x))/2
- *	1. Replace x by |x| (cosh(x) = cosh(-x)). 
- *	2. 
- *		                                        [ exp(x) - 1 ]^2 
+ *	1. Replace x by |x| (cosh(x) = cosh(-x)).
+ *	2.
+ *							[ exp(x) - 1 ]^2
  *	    0        <= x <= ln2/2  :  cosh(x) := 1 + -------------------
- *			       			           2*exp(x)
+ *							   2*exp(x)
  *
- *		                                  exp(x) +  1/exp(x)
+ *						  exp(x) +  1/exp(x)
  *	    ln2/2    <= x <= 22     :  cosh(x) := -------------------
- *			       			          2
- *	    22       <= x <= lnovft :  cosh(x) := exp(x)/2 
+ *							  2
+ *	    22       <= x <= lnovft :  cosh(x) := exp(x)/2
  *	    lnovft   <= x <= ln2ovft:  cosh(x) := exp(x/2)/2 * exp(x/2)
  *	    ln2ovft  <  x	    :  cosh(x) := huge*huge (overflow)
  *
@@ -38,19 +34,11 @@ static char rcsid[] = "$NetBSD: e_cosh.c,v 1.7 1995/05/10 20:44:58 jtc Exp $";
 #include "math.h"
 #include "math_private.h"
 
-#ifdef __STDC__
 static const double one = 1.0, half=0.5, huge = 1.0e300;
-#else
-static double one = 1.0, half=0.5, huge = 1.0e300;
-#endif
 
-#ifdef __STDC__
-	double __ieee754_cosh(double x)
-#else
-	double __ieee754_cosh(x)
-	double x;
-#endif
-{	
+double
+__ieee754_cosh (double x)
+{
 	double t,w;
 	int32_t ix;
 	u_int32_t lx;
@@ -59,19 +47,17 @@ static double one = 1.0, half=0.5, huge = 1.0e300;
 	GET_HIGH_WORD(ix,x);
 	ix &= 0x7fffffff;
 
-    /* x is INF or NaN */
-	if(ix>=0x7ff00000) return x*x;	
-
-    /* |x| in [0,0.5*ln2], return 1+expm1(|x|)^2/(2*exp(|x|)) */
-	if(ix<0x3fd62e43) {
-	    t = __expm1(fabs(x));
-	    w = one+t;
-	    if (ix<0x3c800000) return w;	/* cosh(tiny) = 1 */
-	    return one+(t*t)/(w+w);
-	}
-
-    /* |x| in [0.5*ln2,22], return (exp(|x|)+1/exp(|x|)/2; */
+    /* |x| in [0,22] */
 	if (ix < 0x40360000) {
+	    /* |x| in [0,0.5*ln2], return 1+expm1(|x|)^2/(2*exp(|x|)) */
+		if(ix<0x3fd62e43) {
+		    t = __expm1(fabs(x));
+		    w = one+t;
+		    if (ix<0x3c800000) return w;	/* cosh(tiny) = 1 */
+		    return one+(t*t)/(w+w);
+		}
+
+	    /* |x| in [0.5*ln2,22], return (exp(|x|)+1/exp(|x|)/2; */
 		t = __ieee754_exp(fabs(x));
 		return half*t+half/t;
 	}
@@ -87,6 +73,10 @@ static double one = 1.0, half=0.5, huge = 1.0e300;
 	    return t*w;
 	}
 
+    /* x is INF or NaN */
+	if(ix>=0x7ff00000) return x*x;
+
     /* |x| > overflowthresold, cosh(x) overflow */
 	return huge*huge;
 }
+strong_alias (__ieee754_cosh, __cosh_finite)
diff --git a/libc/sysdeps/ieee754/dbl-64/e_exp.c b/libc/sysdeps/ieee754/dbl-64/e_exp.c
index 717469e25..f4b34a636 100644
--- a/libc/sysdeps/ieee754/dbl-64/e_exp.c
+++ b/libc/sysdeps/ieee754/dbl-64/e_exp.c
@@ -1,7 +1,7 @@
 /*
  * IBM Accurate Mathematical Library
  * written by International Business Machines Corp.
- * Copyright (C) 2001 Free Software Foundation
+ * Copyright (C) 2001, 2011 Free Software Foundation
  *
  * This program is free software; you can redistribute it and/or modify
  * it under the terms of the GNU Lesser General Public License as published by
@@ -145,6 +145,7 @@ double __ieee754_exp(double x) {
     else return __slowexp(x);
   }
 }
+strong_alias (__ieee754_exp, __exp_finite)
 
 /************************************************************************/
 /* Compute e^(x+xx)(Double-Length number) .The routine also receive     */
diff --git a/libc/sysdeps/ieee754/dbl-64/e_exp2.c b/libc/sysdeps/ieee754/dbl-64/e_exp2.c
index ce6368be4..0b7330aac 100644
--- a/libc/sysdeps/ieee754/dbl-64/e_exp2.c
+++ b/libc/sysdeps/ieee754/dbl-64/e_exp2.c
@@ -1,5 +1,6 @@
 /* Double-precision floating point 2^x.
-   Copyright (C) 1997,1998,2000,2001,2005,2006 Free Software Foundation, Inc.
+   Copyright (C) 1997,1998,2000,2001,2005,2006,2011
+   Free Software Foundation, Inc.
    This file is part of the GNU C Library.
    Contributed by Geoffrey Keating <geoffk@ozemail.com.au>
 
@@ -24,9 +25,6 @@
    17 (1), March 1991, pp. 26-45.
    It has been slightly modified to compute 2^x instead of e^x.
    */
-#ifndef _GNU_SOURCE
-#define _GNU_SOURCE
-#endif
 #include <stdlib.h>
 #include <float.h>
 #include <ieee754.h>
@@ -37,13 +35,8 @@
 
 #include "t_exp2.h"
 
-/* XXX I know the assembler generates a warning about incorrect section
-   attributes. But without the attribute here the compiler places the
-   constants in the .data section.  Ideally the constant is placed in
-   .rodata.cst8 so that it can be merged, but gcc sucks, it ICEs when
-   we try to force this section on it.  --drepper  */
-static const volatile double TWO1023 = 8.988465674311579539e+307;
-static const volatile double TWOM1000 = 9.3326361850321887899e-302;
+static const double TWO1023 = 8.988465674311579539e+307;
+static const double TWOM1000 = 9.3326361850321887899e-302;
 
 double
 __ieee754_exp2 (double x)
@@ -52,19 +45,26 @@ __ieee754_exp2 (double x)
   static const double lomark = (double) (DBL_MIN_EXP - DBL_MANT_DIG - 1);
 
   /* Check for usual case.  */
-  if (isless (x, himark) && isgreaterequal (x, lomark))
+  if (__builtin_expect (isless (x, himark), 1))
     {
+      /* Exceptional cases:  */
+      if (__builtin_expect (! isgreaterequal (x, lomark), 0))
+	{
+	  if (__isinf (x))
+	    /* e^-inf == 0, with no error.  */
+	    return 0;
+	  else
+	    /* Underflow */
+	    return TWOM1000 * TWOM1000;
+	}
+
       static const double THREEp42 = 13194139533312.0;
       int tval, unsafe;
       double rx, x22, result;
       union ieee754_double ex2_u, scale_u;
       fenv_t oldenv;
 
-      feholdexcept (&oldenv);
-#ifdef FE_TONEAREST
-      /* If we don't have this, it's too bad.  */
-      fesetround (FE_TONEAREST);
-#endif
+      libc_feholdexcept_setround (&oldenv, FE_TONEAREST);
 
       /* 1. Argument reduction.
 	 Choose integers ex, -256 <= t < 256, and some real
@@ -108,9 +108,10 @@ __ieee754_exp2 (double x)
 	       * x + .055504110254308625)
 	      * x + .240226506959100583)
 	     * x + .69314718055994495) * ex2_u.d;
+      math_opt_barrier (x22);
 
       /* 5. Return (2^x2-1) * 2^(t/512+e+ex) + 2^(t/512+e+ex).  */
-      fesetenv (&oldenv);
+      libc_fesetenv (&oldenv);
 
       result = x22 * x + ex2_u.d;
 
@@ -119,17 +120,8 @@ __ieee754_exp2 (double x)
       else
 	return result * scale_u.d;
     }
-  /* Exceptional cases:  */
-  else if (isless (x, himark))
-    {
-      if (__isinf (x))
-	/* e^-inf == 0, with no error.  */
-	return 0;
-      else
-	/* Underflow */
-	return TWOM1000 * TWOM1000;
-    }
   else
     /* Return x, if x is a NaN or Inf; or overflow, otherwise.  */
     return TWO1023*x;
 }
+strong_alias (__ieee754_exp2, __exp2_finite)
diff --git a/libc/sysdeps/ieee754/dbl-64/e_fmod.c b/libc/sysdeps/ieee754/dbl-64/e_fmod.c
index 2ce613574..a575f616b 100644
--- a/libc/sysdeps/ieee754/dbl-64/e_fmod.c
+++ b/libc/sysdeps/ieee754/dbl-64/e_fmod.c
@@ -5,16 +5,12 @@
  *
  * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  */
 
-#if defined(LIBM_SCCS) && !defined(lint)
-static char rcsid[] = "$NetBSD: e_fmod.c,v 1.8 1995/05/10 20:45:07 jtc Exp $";
-#endif
-
-/* 
+/*
  * __ieee754_fmod(x,y)
  * Return x mod y in exact arithmetic
  * Method: shift and subtract
@@ -23,18 +19,10 @@ static char rcsid[] = "$NetBSD: e_fmod.c,v 1.8 1995/05/10 20:45:07 jtc Exp $";
 #include "math.h"
 #include "math_private.h"
 
-#ifdef __STDC__
 static const double one = 1.0, Zero[] = {0.0, -0.0,};
-#else
-static double one = 1.0, Zero[] = {0.0, -0.0,};
-#endif
 
-#ifdef __STDC__
-	double __ieee754_fmod(double x, double y)
-#else
-	double __ieee754_fmod(x,y)
-	double x,y ;
-#endif
+double
+__ieee754_fmod (double x, double y)
 {
 	int32_t n,hx,hy,hz,ix,iy,sx,i;
 	u_int32_t lx,ly,lz;
@@ -51,7 +39,7 @@ static double one = 1.0, Zero[] = {0.0, -0.0,};
 	    return (x*y)/(x*y);
 	if(hx<=hy) {
 	    if((hx<hy)||(lx<ly)) return x;	/* |x|<|y| return x */
-	    if(lx==ly) 
+	    if(lx==ly)
 		return Zero[(u_int32_t)sx>>31];	/* |x|=|y| return x*0*/
 	}
 
@@ -74,25 +62,25 @@ static double one = 1.0, Zero[] = {0.0, -0.0,};
 	} else iy = (hy>>20)-1023;
 
     /* set up {hx,lx}, {hy,ly} and align y to x */
-	if(ix >= -1022) 
+	if(ix >= -1022)
 	    hx = 0x00100000|(0x000fffff&hx);
 	else {		/* subnormal x, shift x to normal */
 	    n = -1022-ix;
 	    if(n<=31) {
-	        hx = (hx<<n)|(lx>>(32-n));
-	        lx <<= n;
+		hx = (hx<<n)|(lx>>(32-n));
+		lx <<= n;
 	    } else {
 		hx = lx<<(n-32);
 		lx = 0;
 	    }
 	}
-	if(iy >= -1022) 
+	if(iy >= -1022)
 	    hy = 0x00100000|(0x000fffff&hy);
 	else {		/* subnormal y, shift y to normal */
 	    n = -1022-iy;
 	    if(n<=31) {
-	        hy = (hy<<n)|(ly>>(32-n));
-	        ly <<= n;
+		hy = (hy<<n)|(ly>>(32-n));
+		ly <<= n;
 	    } else {
 		hy = ly<<(n-32);
 		ly = 0;
@@ -105,17 +93,17 @@ static double one = 1.0, Zero[] = {0.0, -0.0,};
 	    hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
 	    if(hz<0){hx = hx+hx+(lx>>31); lx = lx+lx;}
 	    else {
-	    	if((hz|lz)==0) 		/* return sign(x)*0 */
+		if((hz|lz)==0)		/* return sign(x)*0 */
 		    return Zero[(u_int32_t)sx>>31];
-	    	hx = hz+hz+(lz>>31); lx = lz+lz;
+		hx = hz+hz+(lz>>31); lx = lz+lz;
 	    }
 	}
 	hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
 	if(hz>=0) {hx=hz;lx=lz;}
 
     /* convert back to floating value and restore the sign */
-	if((hx|lx)==0) 			/* return sign(x)*0 */
-	    return Zero[(u_int32_t)sx>>31];	
+	if((hx|lx)==0)			/* return sign(x)*0 */
+	    return Zero[(u_int32_t)sx>>31];
 	while(hx<0x00100000) {		/* normalize x */
 	    hx = hx+hx+(lx>>31); lx = lx+lx;
 	    iy -= 1;
@@ -138,3 +126,4 @@ static double one = 1.0, Zero[] = {0.0, -0.0,};
 	}
 	return x;		/* exact output */
 }
+strong_alias (__ieee754_fmod, __fmod_finite)
diff --git a/libc/sysdeps/ieee754/dbl-64/e_gamma_r.c b/libc/sysdeps/ieee754/dbl-64/e_gamma_r.c
index f32309350..c4b7470e5 100644
--- a/libc/sysdeps/ieee754/dbl-64/e_gamma_r.c
+++ b/libc/sysdeps/ieee754/dbl-64/e_gamma_r.c
@@ -1,5 +1,5 @@
 /* Implementation of gamma function according to ISO C.
-   Copyright (C) 1997, 1999, 2001, 2004 Free Software Foundation, Inc.
+   Copyright (C) 1997, 1999, 2001, 2004, 2011 Free Software Foundation, Inc.
    This file is part of the GNU C Library.
    Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
 
@@ -33,19 +33,20 @@ __ieee754_gamma_r (double x, int *signgamp)
 
   EXTRACT_WORDS (hx, lx, x);
 
-  if (((hx & 0x7fffffff) | lx) == 0)
+  if (__builtin_expect (((hx & 0x7fffffff) | lx) == 0, 0))
     {
       /* Return value for x == 0 is Inf with divide by zero exception.  */
       *signgamp = 0;
       return 1.0 / x;
     }
-  if (hx < 0 && (u_int32_t) hx < 0xfff00000 && __rint (x) == x)
+  if (__builtin_expect (hx < 0, 0)
+      && (u_int32_t) hx < 0xfff00000 && __rint (x) == x)
     {
       /* Return value for integer x < 0 is NaN with invalid exception.  */
       *signgamp = 0;
       return (x - x) / (x - x);
     }
-  if ((unsigned int) hx == 0xfff00000 && lx==0)
+  if (__builtin_expect ((unsigned int) hx == 0xfff00000 && lx==0, 0))
     {
       /* x == -Inf.  According to ISO this is NaN.  */
       *signgamp = 0;
@@ -55,3 +56,4 @@ __ieee754_gamma_r (double x, int *signgamp)
   /* XXX FIXME.  */
   return __ieee754_exp (__ieee754_lgamma_r (x, signgamp));
 }
+strong_alias (__ieee754_gamma_r, __gamma_r_finite)
diff --git a/libc/sysdeps/ieee754/dbl-64/e_hypot.c b/libc/sysdeps/ieee754/dbl-64/e_hypot.c
index 76a77ec33..a89ccaa35 100644
--- a/libc/sysdeps/ieee754/dbl-64/e_hypot.c
+++ b/libc/sysdeps/ieee754/dbl-64/e_hypot.c
@@ -10,10 +10,6 @@
  * ====================================================
  */
 
-#if defined(LIBM_SCCS) && !defined(lint)
-static char rcsid[] = "$NetBSD: e_hypot.c,v 1.9 1995/05/12 04:57:27 jtc Exp $";
-#endif
-
 /* __ieee754_hypot(x,y)
  *
  * Method :
@@ -42,19 +38,15 @@ static char rcsid[] = "$NetBSD: e_hypot.c,v 1.9 1995/05/12 04:57:27 jtc Exp $";
  *	hypot(x,y) is NAN if x or y is NAN.
  *
  * Accuracy:
- * 	hypot(x,y) returns sqrt(x^2+y^2) with error less
- * 	than 1 ulps (units in the last place)
+ *	hypot(x,y) returns sqrt(x^2+y^2) with error less
+ *	than 1 ulps (units in the last place)
  */
 
 #include "math.h"
 #include "math_private.h"
 
-#ifdef __STDC__
-	double __ieee754_hypot(double x, double y)
-#else
-	double __ieee754_hypot(x,y)
-	double x, y;
-#endif
+double
+__ieee754_hypot(double x, double y)
 {
 	double a,b,t1,t2,y1,y2,w;
 	int32_t j,k,ha,hb;
@@ -68,7 +60,7 @@ static char rcsid[] = "$NetBSD: e_hypot.c,v 1.9 1995/05/12 04:57:27 jtc Exp $";
 	SET_HIGH_WORD(b,hb);	/* b <- |b| */
 	if((ha-hb)>0x3c00000) {return a+b;} /* x/y > 2**60 */
 	k=0;
-	if(ha > 0x5f300000) {	/* a>2**500 */
+	if(__builtin_expect(ha > 0x5f300000, 0)) {	/* a>2**500 */
 	   if(ha >= 0x7ff00000) {	/* Inf or NaN */
 	       u_int32_t low;
 	       w = a+b;			/* for sNaN */
@@ -83,9 +75,9 @@ static char rcsid[] = "$NetBSD: e_hypot.c,v 1.9 1995/05/12 04:57:27 jtc Exp $";
 	   SET_HIGH_WORD(a,ha);
 	   SET_HIGH_WORD(b,hb);
 	}
-	if(hb < 0x20b00000) {	/* b < 2**-500 */
+	if(__builtin_expect(hb < 0x20b00000, 0)) {	/* b < 2**-500 */
 	    if(hb <= 0x000fffff) {	/* subnormal b or 0 */
-	        u_int32_t low;
+		u_int32_t low;
 		GET_LOW_WORD(low,b);
 		if((hb|low)==0) return a;
 		t1=0;
@@ -94,7 +86,7 @@ static char rcsid[] = "$NetBSD: e_hypot.c,v 1.9 1995/05/12 04:57:27 jtc Exp $";
 		a *= t1;
 		k -= 1022;
 	    } else {		/* scale a and b by 2^600 */
-	        ha += 0x25800000; 	/* a *= 2^600 */
+		ha += 0x25800000;	/* a *= 2^600 */
 		hb += 0x25800000;	/* b *= 2^600 */
 		k -= 600;
 		SET_HIGH_WORD(a,ha);
@@ -126,3 +118,4 @@ static char rcsid[] = "$NetBSD: e_hypot.c,v 1.9 1995/05/12 04:57:27 jtc Exp $";
 	    return t1*w;
 	} else return w;
 }
+strong_alias (__ieee754_hypot, __hypot_finite)
diff --git a/libc/sysdeps/ieee754/dbl-64/e_j0.c b/libc/sysdeps/ieee754/dbl-64/e_j0.c
index 302df49d6..5ebf2056b 100644
--- a/libc/sysdeps/ieee754/dbl-64/e_j0.c
+++ b/libc/sysdeps/ieee754/dbl-64/e_j0.c
@@ -13,10 +13,6 @@
    for performance improvement on pipelined processors.
 */
 
-#if defined(LIBM_SCCS) && !defined(lint)
-static char rcsid[] = "$NetBSD: e_j0.c,v 1.8 1995/05/10 20:45:23 jtc Exp $";
-#endif
-
 /* __ieee754_j0(x), __ieee754_y0(x)
  * Bessel function of the first and second kinds of order zero.
  * Method -- j0(x):
@@ -26,16 +22,16 @@ static char rcsid[] = "$NetBSD: e_j0.c,v 1.8 1995/05/10 20:45:23 jtc Exp $";
  *		j0(x) = 1-z/4+ z^2*R0/S0,  where z = x*x;
  *	   (precision:  |j0-1+z/4-z^2R0/S0 |<2**-63.67 )
  *	   for x in (2,inf)
- * 		j0(x) = sqrt(2/(pi*x))*(p0(x)*cos(x0)-q0(x)*sin(x0))
- * 	   where x0 = x-pi/4. It is better to compute sin(x0),cos(x0)
+ *		j0(x) = sqrt(2/(pi*x))*(p0(x)*cos(x0)-q0(x)*sin(x0))
+ *	   where x0 = x-pi/4. It is better to compute sin(x0),cos(x0)
  *	   as follow:
  *		cos(x0) = cos(x)cos(pi/4)+sin(x)sin(pi/4)
  *			= 1/sqrt(2) * (cos(x) + sin(x))
  *		sin(x0) = sin(x)cos(pi/4)-cos(x)sin(pi/4)
  *			= 1/sqrt(2) * (sin(x) - cos(x))
- * 	   (To avoid cancellation, use
+ *	   (To avoid cancellation, use
  *		sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
- * 	    to compute the worse one.)
+ *	    to compute the worse one.)
  *
  *	3 Special cases
  *		j0(nan)= nan
@@ -56,8 +52,8 @@ static char rcsid[] = "$NetBSD: e_j0.c,v 1.8 1995/05/10 20:45:23 jtc Exp $";
  *	   Note: For tiny x, U/V = u0 and j0(x)~1, hence
  *		y0(tiny) = u0 + (2/pi)*ln(tiny), (choose tiny<2**-27)
  *	2. For x>=2.
- * 		y0(x) = sqrt(2/(pi*x))*(p0(x)*cos(x0)+q0(x)*sin(x0))
- * 	   where x0 = x-pi/4. It is better to compute sin(x0),cos(x0)
+ *		y0(x) = sqrt(2/(pi*x))*(p0(x)*cos(x0)+q0(x)*sin(x0))
+ *	   where x0 = x-pi/4. It is better to compute sin(x0),cos(x0)
  *	   by the method mentioned above.
  *	3. Special cases: y0(0)=-inf, y0(x<0)=NaN, y0(inf)=0.
  */
@@ -65,22 +61,14 @@ static char rcsid[] = "$NetBSD: e_j0.c,v 1.8 1995/05/10 20:45:23 jtc Exp $";
 #include "math.h"
 #include "math_private.h"
 
-#ifdef __STDC__
 static double pzero(double), qzero(double);
-#else
-static double pzero(), qzero();
-#endif
 
-#ifdef __STDC__
 static const double
-#else
-static double
-#endif
-huge 	= 1e300,
+huge	= 1e300,
 one	= 1.0,
 invsqrtpi=  5.64189583547756279280e-01, /* 0x3FE20DD7, 0x50429B6D */
 tpi      =  6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */
- 		/* R0/S0 on [0, 2.00] */
+		/* R0/S0 on [0, 2.00] */
 R[]  =  {0.0, 0.0, 1.56249999999999947958e-02, /* 0x3F8FFFFF, 0xFFFFFFFD */
  -1.89979294238854721751e-04, /* 0xBF28E6A5, 0xB61AC6E9 */
   1.82954049532700665670e-06, /* 0x3EBEB1D1, 0x0C503919 */
@@ -90,18 +78,10 @@ S[]  =  {0.0, 1.56191029464890010492e-02, /* 0x3F8FFCE8, 0x82C8C2A4 */
   5.13546550207318111446e-07, /* 0x3EA13B54, 0xCE84D5A9 */
   1.16614003333790000205e-09}; /* 0x3E1408BC, 0xF4745D8F */
 
-#ifdef __STDC__
 static const double zero = 0.0;
-#else
-static double zero = 0.0;
-#endif
 
-#ifdef __STDC__
-	double __ieee754_j0(double x)
-#else
-	double __ieee754_j0(x)
-	double x;
-#endif
+double
+__ieee754_j0(double x)
 {
 	double z, s,c,ss,cc,r,u,v,r1,r2,s1,s2,z2,z4;
 	int32_t hx,ix;
@@ -117,7 +97,7 @@ static double zero = 0.0;
 		if(ix<0x7fe00000) {  /* make sure x+x not overflow */
 		    z = -__cos(x+x);
 		    if ((s*c)<zero) cc = z/ss;
-		    else 	    ss = z/cc;
+		    else	    ss = z/cc;
 		}
 	/*
 	 * j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x)
@@ -132,8 +112,8 @@ static double zero = 0.0;
 	}
 	if(ix<0x3f200000) {	/* |x| < 2**-13 */
 	    if(huge+x>one) {	/* raise inexact if x != 0 */
-	        if(ix<0x3e400000) return one;	/* |x|<2**-27 */
-	        else 	      return one - 0.25*x*x;
+		if(ix<0x3e400000) return one;	/* |x|<2**-27 */
+		else	      return one - 0.25*x*x;
 	    }
 	}
 	z = x*x;
@@ -155,12 +135,9 @@ static double zero = 0.0;
 	    return((one+u)*(one-u)+z*(r/s));
 	}
 }
+strong_alias (__ieee754_j0, __j0_finite)
 
-#ifdef __STDC__
 static const double
-#else
-static double
-#endif
 U[]  = {-7.38042951086872317523e-02, /* 0xBFB2E4D6, 0x99CBD01F */
   1.76666452509181115538e-01, /* 0x3FC69D01, 0x9DE9E3FC */
  -1.38185671945596898896e-02, /* 0xBF8C4CE8, 0xB16CFA97 */
@@ -173,52 +150,48 @@ V[]  =  {1.27304834834123699328e-02, /* 0x3F8A1270, 0x91C9C71A */
   2.59150851840457805467e-07, /* 0x3E91642D, 0x7FF202FD */
   4.41110311332675467403e-10}; /* 0x3DFE5018, 0x3BD6D9EF */
 
-#ifdef __STDC__
-	double __ieee754_y0(double x)
-#else
-	double __ieee754_y0(x)
-	double x;
-#endif
+double
+__ieee754_y0(double x)
 {
 	double z, s,c,ss,cc,u,v,z2,z4,z6,u1,u2,u3,v1,v2;
 	int32_t hx,ix,lx;
 
 	EXTRACT_WORDS(hx,lx,x);
-        ix = 0x7fffffff&hx;
+	ix = 0x7fffffff&hx;
     /* Y0(NaN) is NaN, y0(-inf) is Nan, y0(inf) is 0, y0(0) is -inf.  */
 	if(ix>=0x7ff00000) return  one/(x+x*x);
-        if((ix|lx)==0) return -HUGE_VAL+x; /* -inf and overflow exception.  */
-        if(hx<0) return zero/(zero*x);
-        if(ix >= 0x40000000) {  /* |x| >= 2.0 */
-        /* y0(x) = sqrt(2/(pi*x))*(p0(x)*sin(x0)+q0(x)*cos(x0))
-         * where x0 = x-pi/4
-         *      Better formula:
-         *              cos(x0) = cos(x)cos(pi/4)+sin(x)sin(pi/4)
-         *                      =  1/sqrt(2) * (sin(x) + cos(x))
-         *              sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4)
-         *                      =  1/sqrt(2) * (sin(x) - cos(x))
-         * To avoid cancellation, use
-         *              sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
-         * to compute the worse one.
-         */
+	if((ix|lx)==0) return -HUGE_VAL+x; /* -inf and overflow exception.  */
+	if(hx<0) return zero/(zero*x);
+	if(ix >= 0x40000000) {  /* |x| >= 2.0 */
+	/* y0(x) = sqrt(2/(pi*x))*(p0(x)*sin(x0)+q0(x)*cos(x0))
+	 * where x0 = x-pi/4
+	 *      Better formula:
+	 *              cos(x0) = cos(x)cos(pi/4)+sin(x)sin(pi/4)
+	 *                      =  1/sqrt(2) * (sin(x) + cos(x))
+	 *              sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4)
+	 *                      =  1/sqrt(2) * (sin(x) - cos(x))
+	 * To avoid cancellation, use
+	 *              sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
+	 * to compute the worse one.
+	 */
 		__sincos (x, &s, &c);
-                ss = s-c;
-                cc = s+c;
+		ss = s-c;
+		cc = s+c;
 	/*
 	 * j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x)
 	 * y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x)
 	 */
-                if(ix<0x7fe00000) {  /* make sure x+x not overflow */
-                    z = -__cos(x+x);
-                    if ((s*c)<zero) cc = z/ss;
-                    else            ss = z/cc;
-                }
-                if(ix>0x48000000) z = (invsqrtpi*ss)/__ieee754_sqrt(x);
-                else {
-                    u = pzero(x); v = qzero(x);
-                    z = invsqrtpi*(u*ss+v*cc)/__ieee754_sqrt(x);
-                }
-                return z;
+		if(ix<0x7fe00000) {  /* make sure x+x not overflow */
+		    z = -__cos(x+x);
+		    if ((s*c)<zero) cc = z/ss;
+		    else            ss = z/cc;
+		}
+		if(ix>0x48000000) z = (invsqrtpi*ss)/__ieee754_sqrt(x);
+		else {
+		    u = pzero(x); v = qzero(x);
+		    z = invsqrtpi*(u*ss+v*cc)/__ieee754_sqrt(x);
+		}
+		return z;
 	}
 	if(ix<=0x3e400000) {	/* x < 2**-27 */
 	    return(U[0] + tpi*__ieee754_log(x));
@@ -238,21 +211,18 @@ V[]  =  {1.27304834834123699328e-02, /* 0x3F8A1270, 0x91C9C71A */
 #endif
 	return(u/v + tpi*(__ieee754_j0(x)*__ieee754_log(x)));
 }
+strong_alias (__ieee754_y0, __y0_finite)
 
 /* The asymptotic expansions of pzero is
  *	1 - 9/128 s^2 + 11025/98304 s^4 - ...,	where s = 1/x.
  * For x >= 2, We approximate pzero by
- * 	pzero(x) = 1 + (R/S)
+ *	pzero(x) = 1 + (R/S)
  * where  R = pR0 + pR1*s^2 + pR2*s^4 + ... + pR5*s^10
- * 	  S = 1 + pS0*s^2 + ... + pS4*s^10
+ *	  S = 1 + pS0*s^2 + ... + pS4*s^10
  * and
  *	| pzero(x)-1-R/S | <= 2  ** ( -60.26)
  */
-#ifdef __STDC__
 static const double pR8[6] = { /* for x in [inf, 8]=1/[0,0.125] */
-#else
-static double pR8[6] = { /* for x in [inf, 8]=1/[0,0.125] */
-#endif
   0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */
  -7.03124999999900357484e-02, /* 0xBFB1FFFF, 0xFFFFFD32 */
  -8.08167041275349795626e+00, /* 0xC02029D0, 0xB44FA779 */
@@ -260,11 +230,7 @@ static double pR8[6] = { /* for x in [inf, 8]=1/[0,0.125] */
  -2.48521641009428822144e+03, /* 0xC0A36A6E, 0xCD4DCAFC */
  -5.25304380490729545272e+03, /* 0xC0B4850B, 0x36CC643D */
 };
-#ifdef __STDC__
 static const double pS8[5] = {
-#else
-static double pS8[5] = {
-#endif
   1.16534364619668181717e+02, /* 0x405D2233, 0x07A96751 */
   3.83374475364121826715e+03, /* 0x40ADF37D, 0x50596938 */
   4.05978572648472545552e+04, /* 0x40E3D2BB, 0x6EB6B05F */
@@ -272,11 +238,7 @@ static double pS8[5] = {
   4.76277284146730962675e+04, /* 0x40E74177, 0x4F2C49DC */
 };
 
-#ifdef __STDC__
 static const double pR5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */
-#else
-static double pR5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */
-#endif
  -1.14125464691894502584e-11, /* 0xBDA918B1, 0x47E495CC */
  -7.03124940873599280078e-02, /* 0xBFB1FFFF, 0xE69AFBC6 */
  -4.15961064470587782438e+00, /* 0xC010A370, 0xF90C6BBF */
@@ -284,11 +246,7 @@ static double pR5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */
  -3.31231299649172967747e+02, /* 0xC074B3B3, 0x6742CC63 */
  -3.46433388365604912451e+02, /* 0xC075A6EF, 0x28A38BD7 */
 };
-#ifdef __STDC__
 static const double pS5[5] = {
-#else
-static double pS5[5] = {
-#endif
   6.07539382692300335975e+01, /* 0x404E6081, 0x0C98C5DE */
   1.05125230595704579173e+03, /* 0x40906D02, 0x5C7E2864 */
   5.97897094333855784498e+03, /* 0x40B75AF8, 0x8FBE1D60 */
@@ -296,11 +254,7 @@ static double pS5[5] = {
   2.40605815922939109441e+03, /* 0x40A2CC1D, 0xC70BE864 */
 };
 
-#ifdef __STDC__
 static const double pR3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */
-#else
-static double pR3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */
-#endif
  -2.54704601771951915620e-09, /* 0xBE25E103, 0x6FE1AA86 */
  -7.03119616381481654654e-02, /* 0xBFB1FFF6, 0xF7C0E24B */
  -2.40903221549529611423e+00, /* 0xC00345B2, 0xAEA48074 */
@@ -308,11 +262,7 @@ static double pR3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */
  -5.80791704701737572236e+01, /* 0xC04D0A22, 0x420A1A45 */
  -3.14479470594888503854e+01, /* 0xC03F72AC, 0xA892D80F */
 };
-#ifdef __STDC__
 static const double pS3[5] = {
-#else
-static double pS3[5] = {
-#endif
   3.58560338055209726349e+01, /* 0x4041ED92, 0x84077DD3 */
   3.61513983050303863820e+02, /* 0x40769839, 0x464A7C0E */
   1.19360783792111533330e+03, /* 0x4092A66E, 0x6D1061D6 */
@@ -320,11 +270,7 @@ static double pS3[5] = {
   1.73580930813335754692e+02, /* 0x4065B296, 0xFC379081 */
 };
 
-#ifdef __STDC__
 static const double pR2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */
-#else
-static double pR2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */
-#endif
  -8.87534333032526411254e-08, /* 0xBE77D316, 0xE927026D */
  -7.03030995483624743247e-02, /* 0xBFB1FF62, 0x495E1E42 */
  -1.45073846780952986357e+00, /* 0xBFF73639, 0x8A24A843 */
@@ -332,11 +278,7 @@ static double pR2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */
  -1.11931668860356747786e+01, /* 0xC02662E6, 0xC5246303 */
  -3.23364579351335335033e+00, /* 0xC009DE81, 0xAF8FE70F */
 };
-#ifdef __STDC__
 static const double pS2[5] = {
-#else
-static double pS2[5] = {
-#endif
   2.22202997532088808441e+01, /* 0x40363865, 0x908B5959 */
   1.36206794218215208048e+02, /* 0x4061069E, 0x0EE8878F */
   2.70470278658083486789e+02, /* 0x4070E786, 0x42EA079B */
@@ -344,18 +286,10 @@ static double pS2[5] = {
   1.46576176948256193810e+01, /* 0x402D50B3, 0x44391809 */
 };
 
-#ifdef __STDC__
-	static double pzero(double x)
-#else
-	static double pzero(x)
-	double x;
-#endif
+static double
+pzero(double x)
 {
-#ifdef __STDC__
 	const double *p,*q;
-#else
-	double *p,*q;
-#endif
 	double z,r,s,z2,z4,r1,r2,r3,s1,s2,s3;
 	int32_t ix;
 	GET_HIGH_WORD(ix,x);
@@ -385,17 +319,13 @@ static double pS2[5] = {
 /* For x >= 8, the asymptotic expansions of qzero is
  *	-1/8 s + 75/1024 s^3 - ..., where s = 1/x.
  * We approximate pzero by
- * 	qzero(x) = s*(-1.25 + (R/S))
+ *	qzero(x) = s*(-1.25 + (R/S))
  * where  R = qR0 + qR1*s^2 + qR2*s^4 + ... + qR5*s^10
- * 	  S = 1 + qS0*s^2 + ... + qS5*s^12
+ *	  S = 1 + qS0*s^2 + ... + qS5*s^12
  * and
  *	| qzero(x)/s +1.25-R/S | <= 2  ** ( -61.22)
  */
-#ifdef __STDC__
 static const double qR8[6] = { /* for x in [inf, 8]=1/[0,0.125] */
-#else
-static double qR8[6] = { /* for x in [inf, 8]=1/[0,0.125] */
-#endif
   0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */
   7.32421874999935051953e-02, /* 0x3FB2BFFF, 0xFFFFFE2C */
   1.17682064682252693899e+01, /* 0x40278952, 0x5BB334D6 */
@@ -403,11 +333,7 @@ static double qR8[6] = { /* for x in [inf, 8]=1/[0,0.125] */
   8.85919720756468632317e+03, /* 0x40C14D99, 0x3E18F46D */
   3.70146267776887834771e+04, /* 0x40E212D4, 0x0E901566 */
 };
-#ifdef __STDC__
 static const double qS8[6] = {
-#else
-static double qS8[6] = {
-#endif
   1.63776026895689824414e+02, /* 0x406478D5, 0x365B39BC */
   8.09834494656449805916e+03, /* 0x40BFA258, 0x4E6B0563 */
   1.42538291419120476348e+05, /* 0x41016652, 0x54D38C3F */
@@ -416,11 +342,7 @@ static double qS8[6] = {
  -3.43899293537866615225e+05, /* 0xC114FD6D, 0x2C9530C5 */
 };
 
-#ifdef __STDC__
 static const double qR5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */
-#else
-static double qR5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */
-#endif
   1.84085963594515531381e-11, /* 0x3DB43D8F, 0x29CC8CD9 */
   7.32421766612684765896e-02, /* 0x3FB2BFFF, 0xD172B04C */
   5.83563508962056953777e+00, /* 0x401757B0, 0xB9953DD3 */
@@ -428,11 +350,7 @@ static double qR5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */
   1.02724376596164097464e+03, /* 0x40900CF9, 0x9DC8C481 */
   1.98997785864605384631e+03, /* 0x409F17E9, 0x53C6E3A6 */
 };
-#ifdef __STDC__
 static const double qS5[6] = {
-#else
-static double qS5[6] = {
-#endif
   8.27766102236537761883e+01, /* 0x4054B1B3, 0xFB5E1543 */
   2.07781416421392987104e+03, /* 0x40A03BA0, 0xDA21C0CE */
   1.88472887785718085070e+04, /* 0x40D267D2, 0x7B591E6D */
@@ -441,11 +359,7 @@ static double qS5[6] = {
  -5.35434275601944773371e+03, /* 0xC0B4EA57, 0xBEDBC609 */
 };
 
-#ifdef __STDC__
 static const double qR3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */
-#else
-static double qR3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */
-#endif
   4.37741014089738620906e-09, /* 0x3E32CD03, 0x6ADECB82 */
   7.32411180042911447163e-02, /* 0x3FB2BFEE, 0x0E8D0842 */
   3.34423137516170720929e+00, /* 0x400AC0FC, 0x61149CF5 */
@@ -453,11 +367,7 @@ static double qR3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */
   1.70808091340565596283e+02, /* 0x406559DB, 0xE25EFD1F */
   1.66733948696651168575e+02, /* 0x4064D77C, 0x81FA21E0 */
 };
-#ifdef __STDC__
 static const double qS3[6] = {
-#else
-static double qS3[6] = {
-#endif
   4.87588729724587182091e+01, /* 0x40486122, 0xBFE343A6 */
   7.09689221056606015736e+02, /* 0x40862D83, 0x86544EB3 */
   3.70414822620111362994e+03, /* 0x40ACF04B, 0xE44DFC63 */
@@ -466,11 +376,7 @@ static double qS3[6] = {
  -1.49247451836156386662e+02, /* 0xC062A7EB, 0x201CF40F */
 };
 
-#ifdef __STDC__
 static const double qR2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */
-#else
-static double qR2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */
-#endif
   1.50444444886983272379e-07, /* 0x3E84313B, 0x54F76BDB */
   7.32234265963079278272e-02, /* 0x3FB2BEC5, 0x3E883E34 */
   1.99819174093815998816e+00, /* 0x3FFFF897, 0xE727779C */
@@ -478,11 +384,7 @@ static double qR2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */
   3.16662317504781540833e+01, /* 0x403FAA8E, 0x29FBDC4A */
   1.62527075710929267416e+01, /* 0x403040B1, 0x71814BB4 */
 };
-#ifdef __STDC__
 static const double qS2[6] = {
-#else
-static double qS2[6] = {
-#endif
   3.03655848355219184498e+01, /* 0x403E5D96, 0xF7C07AED */
   2.69348118608049844624e+02, /* 0x4070D591, 0xE4D14B40 */
   8.44783757595320139444e+02, /* 0x408A6645, 0x22B3BF22 */
@@ -491,18 +393,10 @@ static double qS2[6] = {
  -5.31095493882666946917e+00, /* 0xC0153E6A, 0xF8B32931 */
 };
 
-#ifdef __STDC__
-	static double qzero(double x)
-#else
-	static double qzero(x)
-	double x;
-#endif
+static double
+qzero(double x)
 {
-#ifdef __STDC__
 	const double *p,*q;
-#else
-	double *p,*q;
-#endif
 	double s,r,z,z2,z4,z6,r1,r2,r3,s1,s2,s3;
 	int32_t ix;
 	GET_HIGH_WORD(ix,x);
diff --git a/libc/sysdeps/ieee754/dbl-64/e_j1.c b/libc/sysdeps/ieee754/dbl-64/e_j1.c
index 8a3b2ffd1..fdc6b5b89 100644
--- a/libc/sysdeps/ieee754/dbl-64/e_j1.c
+++ b/libc/sysdeps/ieee754/dbl-64/e_j1.c
@@ -13,10 +13,6 @@
    for performance improvement on pipelined processors.
 */
 
-#if defined(LIBM_SCCS) && !defined(lint)
-static char rcsid[] = "$NetBSD: e_j1.c,v 1.8 1995/05/10 20:45:27 jtc Exp $";
-#endif
-
 /* __ieee754_j1(x), __ieee754_y1(x)
  * Bessel function of the first and second kinds of order zero.
  * Method -- j1(x):
@@ -26,17 +22,17 @@ static char rcsid[] = "$NetBSD: e_j1.c,v 1.8 1995/05/10 20:45:27 jtc Exp $";
  *		j1(x) = x/2 + x*z*R0/S0,  where z = x*x;
  *	   (precision:  |j1/x - 1/2 - R0/S0 |<2**-61.51 )
  *	   for x in (2,inf)
- * 		j1(x) = sqrt(2/(pi*x))*(p1(x)*cos(x1)-q1(x)*sin(x1))
- * 		y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x1)+q1(x)*cos(x1))
- * 	   where x1 = x-3*pi/4. It is better to compute sin(x1),cos(x1)
+ *		j1(x) = sqrt(2/(pi*x))*(p1(x)*cos(x1)-q1(x)*sin(x1))
+ *		y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x1)+q1(x)*cos(x1))
+ *	   where x1 = x-3*pi/4. It is better to compute sin(x1),cos(x1)
  *	   as follow:
  *		cos(x1) =  cos(x)cos(3pi/4)+sin(x)sin(3pi/4)
  *			=  1/sqrt(2) * (sin(x) - cos(x))
  *		sin(x1) =  sin(x)cos(3pi/4)-cos(x)sin(3pi/4)
  *			= -1/sqrt(2) * (sin(x) + cos(x))
- * 	   (To avoid cancellation, use
+ *	   (To avoid cancellation, use
  *		sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
- * 	    to compute the worse one.)
+ *	    to compute the worse one.)
  *
  *	3 Special cases
  *		j1(nan)= nan
@@ -57,25 +53,17 @@ static char rcsid[] = "$NetBSD: e_j1.c,v 1.8 1995/05/10 20:45:27 jtc Exp $";
  *	   Note: For tiny x, 1/x dominate y1 and hence
  *		y1(tiny) = -2/pi/tiny, (choose tiny<2**-54)
  *	3. For x>=2.
- * 		y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x1)+q1(x)*cos(x1))
- * 	   where x1 = x-3*pi/4. It is better to compute sin(x1),cos(x1)
+ *		y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x1)+q1(x)*cos(x1))
+ *	   where x1 = x-3*pi/4. It is better to compute sin(x1),cos(x1)
  *	   by method mentioned above.
  */
 
 #include "math.h"
 #include "math_private.h"
 
-#ifdef __STDC__
 static double pone(double), qone(double);
-#else
-static double pone(), qone();
-#endif
 
-#ifdef __STDC__
 static const double
-#else
-static double
-#endif
 huge    = 1e300,
 one	= 1.0,
 invsqrtpi=  5.64189583547756279280e-01, /* 0x3FE20DD7, 0x50429B6D */
@@ -91,25 +79,17 @@ S[]  =  {0.0, 1.91537599538363460805e-02, /* 0x3F939D0B, 0x12637E53 */
   5.04636257076217042715e-09, /* 0x3E35AC88, 0xC97DFF2C */
   1.23542274426137913908e-11}; /* 0x3DAB2ACF, 0xCFB97ED8 */
 
-#ifdef __STDC__
 static const double zero    = 0.0;
-#else
-static double zero    = 0.0;
-#endif
 
-#ifdef __STDC__
-	double __ieee754_j1(double x)
-#else
-	double __ieee754_j1(x)
-	double x;
-#endif
+double
+__ieee754_j1(double x)
 {
 	double z, s,c,ss,cc,r,u,v,y,r1,r2,s1,s2,s3,z2,z4;
 	int32_t hx,ix;
 
 	GET_HIGH_WORD(hx,x);
 	ix = hx&0x7fffffff;
-	if(ix>=0x7ff00000) return one/x;
+	if(__builtin_expect(ix>=0x7ff00000, 0)) return one/x;
 	y = fabs(x);
 	if(ix >= 0x40000000) {	/* |x| >= 2.0 */
 		__sincos (y, &s, &c);
@@ -118,7 +98,7 @@ static double zero    = 0.0;
 		if(ix<0x7fe00000) {  /* make sure y+y not overflow */
 		    z = __cos(y+y);
 		    if ((s*c)>zero) cc = z/ss;
-		    else 	    ss = z/cc;
+		    else	    ss = z/cc;
 		}
 	/*
 	 * j1(x) = 1/sqrt(pi) * (P(1,x)*cc - Q(1,x)*ss) / sqrt(x)
@@ -130,9 +110,9 @@ static double zero    = 0.0;
 		    z = invsqrtpi*(u*cc-v*ss)/__ieee754_sqrt(y);
 		}
 		if(hx<0) return -z;
-		else  	 return  z;
+		else	 return  z;
 	}
-	if(ix<0x3e400000) {	/* |x|<2**-27 */
+	if(__builtin_expect(ix<0x3e400000, 0)) {	/* |x|<2**-27 */
 	    if(huge+x>one) return 0.5*x;/* inexact if x!=0 necessary */
 	}
 	z = x*x;
@@ -144,7 +124,7 @@ static double zero    = 0.0;
 	r1 = z*R[0]; z2=z*z;
 	r2 = R[1]+z*R[2]; z4=z2*z2;
 	r = r1 + z2*r2 + z4*R[3];
-  	r *= x;
+	r *= x;
 	s1 = one+z*S[1];
 	s2 = S[2]+z*S[3];
 	s3 = S[4]+z*S[5];
@@ -152,23 +132,16 @@ static double zero    = 0.0;
 #endif
 	return(x*0.5+r/s);
 }
+strong_alias (__ieee754_j1, __j1_finite)
 
-#ifdef __STDC__
 static const double U0[5] = {
-#else
-static double U0[5] = {
-#endif
  -1.96057090646238940668e-01, /* 0xBFC91866, 0x143CBC8A */
   5.04438716639811282616e-02, /* 0x3FA9D3C7, 0x76292CD1 */
  -1.91256895875763547298e-03, /* 0xBF5F55E5, 0x4844F50F */
   2.35252600561610495928e-05, /* 0x3EF8AB03, 0x8FA6B88E */
  -9.19099158039878874504e-08, /* 0xBE78AC00, 0x569105B8 */
 };
-#ifdef __STDC__
 static const double V0[5] = {
-#else
-static double V0[5] = {
-#endif
   1.99167318236649903973e-02, /* 0x3F94650D, 0x3F4DA9F0 */
   2.02552581025135171496e-04, /* 0x3F2A8C89, 0x6C257764 */
   1.35608801097516229404e-06, /* 0x3EB6C05A, 0x894E8CA6 */
@@ -176,56 +149,53 @@ static double V0[5] = {
   1.66559246207992079114e-11, /* 0x3DB25039, 0xDACA772A */
 };
 
-#ifdef __STDC__
-	double __ieee754_y1(double x)
-#else
-	double __ieee754_y1(x)
-	double x;
-#endif
+double
+__ieee754_y1(double x)
 {
 	double z, s,c,ss,cc,u,v,u1,u2,v1,v2,v3,z2,z4;
 	int32_t hx,ix,lx;
 
 	EXTRACT_WORDS(hx,lx,x);
-        ix = 0x7fffffff&hx;
+	ix = 0x7fffffff&hx;
     /* if Y1(NaN) is NaN, Y1(-inf) is NaN, Y1(inf) is 0 */
-	if(ix>=0x7ff00000) return  one/(x+x*x);
-        if((ix|lx)==0) return -HUGE_VAL+x; /* -inf and overflow exception.  */;
-        if(hx<0) return zero/(zero*x);
-        if(ix >= 0x40000000) {  /* |x| >= 2.0 */
+	if(__builtin_expect(ix>=0x7ff00000, 0)) return  one/(x+x*x);
+	if(__builtin_expect((ix|lx)==0, 0))
+		return -HUGE_VAL+x; /* -inf and overflow exception.  */;
+	if(__builtin_expect(hx<0, 0)) return zero/(zero*x);
+	if(ix >= 0x40000000) {  /* |x| >= 2.0 */
 		__sincos (x, &s, &c);
-                ss = -s-c;
-                cc = s-c;
-                if(ix<0x7fe00000) {  /* make sure x+x not overflow */
-                    z = __cos(x+x);
-                    if ((s*c)>zero) cc = z/ss;
-                    else            ss = z/cc;
-                }
-        /* y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x0)+q1(x)*cos(x0))
-         * where x0 = x-3pi/4
-         *      Better formula:
-         *              cos(x0) = cos(x)cos(3pi/4)+sin(x)sin(3pi/4)
-         *                      =  1/sqrt(2) * (sin(x) - cos(x))
-         *              sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4)
-         *                      = -1/sqrt(2) * (cos(x) + sin(x))
-         * To avoid cancellation, use
-         *              sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
-         * to compute the worse one.
-         */
-                if(ix>0x48000000) z = (invsqrtpi*ss)/__ieee754_sqrt(x);
-                else {
-                    u = pone(x); v = qone(x);
-                    z = invsqrtpi*(u*ss+v*cc)/__ieee754_sqrt(x);
-                }
-                return z;
-        }
-        if(ix<=0x3c900000) {    /* x < 2**-54 */
-            return(-tpi/x);
-        }
-        z = x*x;
+		ss = -s-c;
+		cc = s-c;
+		if(ix<0x7fe00000) {  /* make sure x+x not overflow */
+		    z = __cos(x+x);
+		    if ((s*c)>zero) cc = z/ss;
+		    else            ss = z/cc;
+		}
+	/* y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x0)+q1(x)*cos(x0))
+	 * where x0 = x-3pi/4
+	 *      Better formula:
+	 *              cos(x0) = cos(x)cos(3pi/4)+sin(x)sin(3pi/4)
+	 *                      =  1/sqrt(2) * (sin(x) - cos(x))
+	 *              sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4)
+	 *                      = -1/sqrt(2) * (cos(x) + sin(x))
+	 * To avoid cancellation, use
+	 *              sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
+	 * to compute the worse one.
+	 */
+		if(ix>0x48000000) z = (invsqrtpi*ss)/__ieee754_sqrt(x);
+		else {
+		    u = pone(x); v = qone(x);
+		    z = invsqrtpi*(u*ss+v*cc)/__ieee754_sqrt(x);
+		}
+		return z;
+	}
+	if(__builtin_expect(ix<=0x3c900000, 0)) {    /* x < 2**-54 */
+	    return(-tpi/x);
+	}
+	z = x*x;
 #ifdef DO_NOT_USE_THIS
-        u = U0[0]+z*(U0[1]+z*(U0[2]+z*(U0[3]+z*U0[4])));
-        v = one+z*(V0[0]+z*(V0[1]+z*(V0[2]+z*(V0[3]+z*V0[4]))));
+	u = U0[0]+z*(U0[1]+z*(U0[2]+z*(U0[3]+z*U0[4])));
+	v = one+z*(V0[0]+z*(V0[1]+z*(V0[2]+z*(V0[3]+z*V0[4]))));
 #else
 	u1 = U0[0]+z*U0[1];z2=z*z;
 	u2 = U0[2]+z*U0[3];z4=z2*z2;
@@ -235,24 +205,21 @@ static double V0[5] = {
 	v3 = V0[3]+z*V0[4];
 	v = v1 + z2*v2 + z4*v3;
 #endif
-        return(x*(u/v) + tpi*(__ieee754_j1(x)*__ieee754_log(x)-one/x));
+	return(x*(u/v) + tpi*(__ieee754_j1(x)*__ieee754_log(x)-one/x));
 }
+strong_alias (__ieee754_y1, __y1_finite)
 
 /* For x >= 8, the asymptotic expansions of pone is
  *	1 + 15/128 s^2 - 4725/2^15 s^4 - ...,	where s = 1/x.
  * We approximate pone by
- * 	pone(x) = 1 + (R/S)
+ *	pone(x) = 1 + (R/S)
  * where  R = pr0 + pr1*s^2 + pr2*s^4 + ... + pr5*s^10
- * 	  S = 1 + ps0*s^2 + ... + ps4*s^10
+ *	  S = 1 + ps0*s^2 + ... + ps4*s^10
  * and
  *	| pone(x)-1-R/S | <= 2  ** ( -60.06)
  */
 
-#ifdef __STDC__
 static const double pr8[6] = { /* for x in [inf, 8]=1/[0,0.125] */
-#else
-static double pr8[6] = { /* for x in [inf, 8]=1/[0,0.125] */
-#endif
   0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */
   1.17187499999988647970e-01, /* 0x3FBDFFFF, 0xFFFFFCCE */
   1.32394806593073575129e+01, /* 0x402A7A9D, 0x357F7FCE */
@@ -260,11 +227,7 @@ static double pr8[6] = { /* for x in [inf, 8]=1/[0,0.125] */
   3.87474538913960532227e+03, /* 0x40AE457D, 0xA3A532CC */
   7.91447954031891731574e+03, /* 0x40BEEA7A, 0xC32782DD */
 };
-#ifdef __STDC__
 static const double ps8[5] = {
-#else
-static double ps8[5] = {
-#endif
   1.14207370375678408436e+02, /* 0x405C8D45, 0x8E656CAC */
   3.65093083420853463394e+03, /* 0x40AC85DC, 0x964D274F */
   3.69562060269033463555e+04, /* 0x40E20B86, 0x97C5BB7F */
@@ -272,11 +235,7 @@ static double ps8[5] = {
   3.08042720627888811578e+04, /* 0x40DE1511, 0x697A0B2D */
 };
 
-#ifdef __STDC__
 static const double pr5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */
-#else
-static double pr5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */
-#endif
   1.31990519556243522749e-11, /* 0x3DAD0667, 0xDAE1CA7D */
   1.17187493190614097638e-01, /* 0x3FBDFFFF, 0xE2C10043 */
   6.80275127868432871736e+00, /* 0x401B3604, 0x6E6315E3 */
@@ -284,11 +243,7 @@ static double pr5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */
   5.17636139533199752805e+02, /* 0x40802D16, 0xD052D649 */
   5.28715201363337541807e+02, /* 0x408085B8, 0xBB7E0CB7 */
 };
-#ifdef __STDC__
 static const double ps5[5] = {
-#else
-static double ps5[5] = {
-#endif
   5.92805987221131331921e+01, /* 0x404DA3EA, 0xA8AF633D */
   9.91401418733614377743e+02, /* 0x408EFB36, 0x1B066701 */
   5.35326695291487976647e+03, /* 0x40B4E944, 0x5706B6FB */
@@ -296,11 +251,7 @@ static double ps5[5] = {
   1.50404688810361062679e+03, /* 0x40978030, 0x036F5E51 */
 };
 
-#ifdef __STDC__
 static const double pr3[6] = {
-#else
-static double pr3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */
-#endif
   3.02503916137373618024e-09, /* 0x3E29FC21, 0xA7AD9EDD */
   1.17186865567253592491e-01, /* 0x3FBDFFF5, 0x5B21D17B */
   3.93297750033315640650e+00, /* 0x400F76BC, 0xE85EAD8A */
@@ -308,11 +259,7 @@ static double pr3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */
   9.10550110750781271918e+01, /* 0x4056C385, 0x4D2C1837 */
   4.85590685197364919645e+01, /* 0x4048478F, 0x8EA83EE5 */
 };
-#ifdef __STDC__
 static const double ps3[5] = {
-#else
-static double ps3[5] = {
-#endif
   3.47913095001251519989e+01, /* 0x40416549, 0xA134069C */
   3.36762458747825746741e+02, /* 0x40750C33, 0x07F1A75F */
   1.04687139975775130551e+03, /* 0x40905B7C, 0x5037D523 */
@@ -320,11 +267,7 @@ static double ps3[5] = {
   1.03787932439639277504e+02, /* 0x4059F26D, 0x7C2EED53 */
 };
 
-#ifdef __STDC__
 static const double pr2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */
-#else
-static double pr2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */
-#endif
   1.07710830106873743082e-07, /* 0x3E7CE9D4, 0xF65544F4 */
   1.17176219462683348094e-01, /* 0x3FBDFF42, 0xBE760D83 */
   2.36851496667608785174e+00, /* 0x4002F2B7, 0xF98FAEC0 */
@@ -332,11 +275,7 @@ static double pr2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */
   1.76939711271687727390e+01, /* 0x4031B1A8, 0x177F8EE2 */
   5.07352312588818499250e+00, /* 0x40144B49, 0xA574C1FE */
 };
-#ifdef __STDC__
 static const double ps2[5] = {
-#else
-static double ps2[5] = {
-#endif
   2.14364859363821409488e+01, /* 0x40356FBD, 0x8AD5ECDC */
   1.25290227168402751090e+02, /* 0x405F5293, 0x14F92CD5 */
   2.32276469057162813669e+02, /* 0x406D08D8, 0xD5A2DBD9 */
@@ -344,30 +283,22 @@ static double ps2[5] = {
   8.36463893371618283368e+00, /* 0x4020BAB1, 0xF44E5192 */
 };
 
-#ifdef __STDC__
-	static double pone(double x)
-#else
-	static double pone(x)
-	double x;
-#endif
+static double
+pone(double x)
 {
-#ifdef __STDC__
 	const double *p,*q;
-#else
-	double *p,*q;
-#endif
 	double z,r,s,r1,r2,r3,s1,s2,s3,z2,z4;
-        int32_t ix;
+	int32_t ix;
 	GET_HIGH_WORD(ix,x);
 	ix &= 0x7fffffff;
-        if(ix>=0x40200000)     {p = pr8; q= ps8;}
-        else if(ix>=0x40122E8B){p = pr5; q= ps5;}
-        else if(ix>=0x4006DB6D){p = pr3; q= ps3;}
-        else if(ix>=0x40000000){p = pr2; q= ps2;}
-        z = one/(x*x);
+	if(ix>=0x40200000)     {p = pr8; q= ps8;}
+	else if(ix>=0x40122E8B){p = pr5; q= ps5;}
+	else if(ix>=0x4006DB6D){p = pr3; q= ps3;}
+	else if(ix>=0x40000000){p = pr2; q= ps2;}
+	z = one/(x*x);
 #ifdef DO_NOT_USE_THIS
-        r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
-        s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4]))));
+	r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
+	s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4]))));
 #else
 	r1 = p[0]+z*p[1]; z2=z*z;
 	r2 = p[2]+z*p[3]; z4=z2*z2;
@@ -378,25 +309,21 @@ static double ps2[5] = {
 	s3 = q[3]+z*q[4];
 	s = s1 + z2*s2 + z4*s3;
 #endif
-        return one+ r/s;
+	return one+ r/s;
 }
 
 
 /* For x >= 8, the asymptotic expansions of qone is
  *	3/8 s - 105/1024 s^3 - ..., where s = 1/x.
  * We approximate pone by
- * 	qone(x) = s*(0.375 + (R/S))
+ *	qone(x) = s*(0.375 + (R/S))
  * where  R = qr1*s^2 + qr2*s^4 + ... + qr5*s^10
- * 	  S = 1 + qs1*s^2 + ... + qs6*s^12
+ *	  S = 1 + qs1*s^2 + ... + qs6*s^12
  * and
  *	| qone(x)/s -0.375-R/S | <= 2  ** ( -61.13)
  */
 
-#ifdef __STDC__
 static const double qr8[6] = { /* for x in [inf, 8]=1/[0,0.125] */
-#else
-static double qr8[6] = { /* for x in [inf, 8]=1/[0,0.125] */
-#endif
   0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */
  -1.02539062499992714161e-01, /* 0xBFBA3FFF, 0xFFFFFDF3 */
  -1.62717534544589987888e+01, /* 0xC0304591, 0xA26779F7 */
@@ -404,11 +331,7 @@ static double qr8[6] = { /* for x in [inf, 8]=1/[0,0.125] */
  -1.18498066702429587167e+04, /* 0xC0C724E7, 0x40F87415 */
  -4.84385124285750353010e+04, /* 0xC0E7A6D0, 0x65D09C6A */
 };
-#ifdef __STDC__
 static const double qs8[6] = {
-#else
-static double qs8[6] = {
-#endif
   1.61395369700722909556e+02, /* 0x40642CA6, 0xDE5BCDE5 */
   7.82538599923348465381e+03, /* 0x40BE9162, 0xD0D88419 */
   1.33875336287249578163e+05, /* 0x4100579A, 0xB0B75E98 */
@@ -417,11 +340,7 @@ static double qs8[6] = {
  -2.94490264303834643215e+05, /* 0xC111F969, 0x0EA5AA18 */
 };
 
-#ifdef __STDC__
 static const double qr5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */
-#else
-static double qr5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */
-#endif
  -2.08979931141764104297e-11, /* 0xBDB6FA43, 0x1AA1A098 */
  -1.02539050241375426231e-01, /* 0xBFBA3FFF, 0xCB597FEF */
  -8.05644828123936029840e+00, /* 0xC0201CE6, 0xCA03AD4B */
@@ -429,11 +348,7 @@ static double qr5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */
  -1.37319376065508163265e+03, /* 0xC09574C6, 0x6931734F */
  -2.61244440453215656817e+03, /* 0xC0A468E3, 0x88FDA79D */
 };
-#ifdef __STDC__
 static const double qs5[6] = {
-#else
-static double qs5[6] = {
-#endif
   8.12765501384335777857e+01, /* 0x405451B2, 0xFF5A11B2 */
   1.99179873460485964642e+03, /* 0x409F1F31, 0xE77BF839 */
   1.74684851924908907677e+04, /* 0x40D10F1F, 0x0D64CE29 */
@@ -442,11 +357,7 @@ static double qs5[6] = {
  -4.71918354795128470869e+03, /* 0xC0B26F2E, 0xFCFFA004 */
 };
 
-#ifdef __STDC__
 static const double qr3[6] = {
-#else
-static double qr3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */
-#endif
  -5.07831226461766561369e-09, /* 0xBE35CFA9, 0xD38FC84F */
  -1.02537829820837089745e-01, /* 0xBFBA3FEB, 0x51AEED54 */
  -4.61011581139473403113e+00, /* 0xC01270C2, 0x3302D9FF */
@@ -454,11 +365,7 @@ static double qr3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */
  -2.28244540737631695038e+02, /* 0xC06C87D3, 0x4718D55F */
  -2.19210128478909325622e+02, /* 0xC06B66B9, 0x5F5C1BF6 */
 };
-#ifdef __STDC__
 static const double qs3[6] = {
-#else
-static double qs3[6] = {
-#endif
   4.76651550323729509273e+01, /* 0x4047D523, 0xCCD367E4 */
   6.73865112676699709482e+02, /* 0x40850EEB, 0xC031EE3E */
   3.38015286679526343505e+03, /* 0x40AA684E, 0x448E7C9A */
@@ -467,11 +374,7 @@ static double qs3[6] = {
  -1.35201191444307340817e+02, /* 0xC060E670, 0x290A311F */
 };
 
-#ifdef __STDC__
 static const double qr2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */
-#else
-static double qr2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */
-#endif
  -1.78381727510958865572e-07, /* 0xBE87F126, 0x44C626D2 */
  -1.02517042607985553460e-01, /* 0xBFBA3E8E, 0x9148B010 */
  -2.75220568278187460720e+00, /* 0xC0060484, 0x69BB4EDA */
@@ -479,11 +382,7 @@ static double qr2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */
  -4.23253133372830490089e+01, /* 0xC04529A3, 0xDE104AAA */
  -2.13719211703704061733e+01, /* 0xC0355F36, 0x39CF6E52 */
 };
-#ifdef __STDC__
 static const double qs2[6] = {
-#else
-static double qs2[6] = {
-#endif
   2.95333629060523854548e+01, /* 0x403D888A, 0x78AE64FF */
   2.52981549982190529136e+02, /* 0x406F9F68, 0xDB821CBA */
   7.57502834868645436472e+02, /* 0x4087AC05, 0xCE49A0F7 */
@@ -492,18 +391,10 @@ static double qs2[6] = {
  -4.95949898822628210127e+00, /* 0xC013D686, 0xE71BE86B */
 };
 
-#ifdef __STDC__
-	static double qone(double x)
-#else
-	static double qone(x)
-	double x;
-#endif
+static double
+qone(double x)
 {
-#ifdef __STDC__
 	const double *p,*q;
-#else
-	double *p,*q;
-#endif
 	double  s,r,z,r1,r2,r3,s1,s2,s3,z2,z4,z6;
 	int32_t ix;
 	GET_HIGH_WORD(ix,x);
diff --git a/libc/sysdeps/ieee754/dbl-64/e_jn.c b/libc/sysdeps/ieee754/dbl-64/e_jn.c
index bf4a13d97..f8b8e2ee6 100644
--- a/libc/sysdeps/ieee754/dbl-64/e_jn.c
+++ b/libc/sysdeps/ieee754/dbl-64/e_jn.c
@@ -10,10 +10,6 @@
  * ====================================================
  */
 
-#if defined(LIBM_SCCS) && !defined(lint)
-static char rcsid[] = "$NetBSD: e_jn.c,v 1.9 1995/05/10 20:45:34 jtc Exp $";
-#endif
-
 /*
  * __ieee754_jn(n, x), __ieee754_yn(n, x)
  * floating point Bessel's function of the 1st and 2nd kind
@@ -43,27 +39,15 @@ static char rcsid[] = "$NetBSD: e_jn.c,v 1.9 1995/05/10 20:45:34 jtc Exp $";
 #include "math.h"
 #include "math_private.h"
 
-#ifdef __STDC__
 static const double
-#else
-static double
-#endif
 invsqrtpi=  5.64189583547756279280e-01, /* 0x3FE20DD7, 0x50429B6D */
 two   =  2.00000000000000000000e+00, /* 0x40000000, 0x00000000 */
 one   =  1.00000000000000000000e+00; /* 0x3FF00000, 0x00000000 */
 
-#ifdef __STDC__
 static const double zero  =  0.00000000000000000000e+00;
-#else
-static double zero  =  0.00000000000000000000e+00;
-#endif
 
-#ifdef __STDC__
-	double __ieee754_jn(int n, double x)
-#else
-	double __ieee754_jn(n,x)
-	int n; double x;
-#endif
+double
+__ieee754_jn(int n, double x)
 {
 	int32_t i,hx,ix,lx, sgn;
 	double a, b, temp, di;
@@ -75,7 +59,8 @@ static double zero  =  0.00000000000000000000e+00;
 	EXTRACT_WORDS(hx,lx,x);
 	ix = 0x7fffffff&hx;
     /* if J(n,NaN) is NaN */
-	if((ix|((u_int32_t)(lx|-lx))>>31)>0x7ff00000) return x+x;
+	if(__builtin_expect((ix|((u_int32_t)(lx|-lx))>>31)>0x7ff00000, 0))
+		return x+x;
 	if(n<0){
 		n = -n;
 		x = -x;
@@ -85,8 +70,9 @@ static double zero  =  0.00000000000000000000e+00;
 	if(n==1) return(__ieee754_j1(x));
 	sgn = (n&1)&(hx>>31);	/* even n -- 0, odd n -- sign(x) */
 	x = fabs(x);
-	if((ix|lx)==0||ix>=0x7ff00000) 	/* if x is 0 or inf */
-	    b = zero;
+	if(__builtin_expect((ix|lx)==0||ix>=0x7ff00000,0))
+		/* if x is 0 or inf */
+		b = zero;
 	else if((double)n<=x) {
 		/* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */
 	    if(ix>=0x52D00000) { /* x > 2**302 */
@@ -99,7 +85,7 @@ static double zero  =  0.00000000000000000000e+00;
      *		   n	sin(xn)*sqt2	cos(xn)*sqt2
      *		----------------------------------
      *		   0	 s-c		 c+s
-     *		   1	-s-c 		-c+s
+     *		   1	-s-c		-c+s
      *		   2	-s+c		-c-s
      *		   3	 s+c		 c-s
      */
@@ -114,13 +100,13 @@ static double zero  =  0.00000000000000000000e+00;
 		}
 		b = invsqrtpi*temp/__ieee754_sqrt(x);
 	    } else {
-	        a = __ieee754_j0(x);
-	        b = __ieee754_j1(x);
-	        for(i=1;i<n;i++){
+		a = __ieee754_j0(x);
+		b = __ieee754_j1(x);
+		for(i=1;i<n;i++){
 		    temp = b;
 		    b = b*((double)(i+i)/x) - a; /* avoid underflow */
 		    a = temp;
-	        }
+		}
 	    }
 	} else {
 	    if(ix<0x3e100000) {	/* x < 2**-29 */
@@ -139,11 +125,11 @@ static double zero  =  0.00000000000000000000e+00;
 		}
 	    } else {
 		/* use backward recurrence */
-		/* 			x      x^2      x^2
+		/*			x      x^2      x^2
 		 *  J(n,x)/J(n-1,x) =  ----   ------   ------   .....
 		 *			2n  - 2(n+1) - 2(n+2)
 		 *
-		 * 			1      1        1
+		 *			1      1        1
 		 *  (for large x)   =  ----  ------   ------   .....
 		 *			2n   2(n+1)   2(n+2)
 		 *			-- - ------ - ------ -
@@ -156,7 +142,7 @@ static double zero  =  0.00000000000000000000e+00;
 		 *		       1
 		 *	   w - -----------------
 		 *			  1
-		 * 	        w+h - ---------
+		 *		w+h - ---------
 		 *		       w+2h - ...
 		 *
 		 * To determine how many terms needed, let
@@ -193,19 +179,19 @@ static double zero  =  0.00000000000000000000e+00;
 		v = two/x;
 		tmp = tmp*__ieee754_log(fabs(v*tmp));
 		if(tmp<7.09782712893383973096e+02) {
-	    	    for(i=n-1,di=(double)(i+i);i>0;i--){
-		        temp = b;
+		    for(i=n-1,di=(double)(i+i);i>0;i--){
+			temp = b;
 			b *= di;
 			b  = b/x - a;
-		        a = temp;
+			a = temp;
 			di -= two;
-	     	    }
+		    }
 		} else {
-	    	    for(i=n-1,di=(double)(i+i);i>0;i--){
-		        temp = b;
+		    for(i=n-1,di=(double)(i+i);i>0;i--){
+			temp = b;
 			b *= di;
 			b  = b/x - a;
-		        a = temp;
+			a = temp;
 			di -= two;
 		    /* scale b to avoid spurious overflow */
 			if(b>1e100) {
@@ -213,20 +199,26 @@ static double zero  =  0.00000000000000000000e+00;
 			    t /= b;
 			    b  = one;
 			}
-	     	    }
+		    }
 		}
-	    	b = (t*__ieee754_j0(x)/b);
+		/* j0() and j1() suffer enormous loss of precision at and
+		 * near zero; however, we know that their zero points never
+		 * coincide, so just choose the one further away from zero.
+		 */
+		z = __ieee754_j0 (x);
+		w = __ieee754_j1 (x);
+		if (fabs (z) >= fabs (w))
+		  b = (t * z / b);
+		else
+		  b = (t * w / a);
 	    }
 	}
 	if(sgn==1) return -b; else return b;
 }
+strong_alias (__ieee754_jn, __jn_finite)
 
-#ifdef __STDC__
-	double __ieee754_yn(int n, double x)
-#else
-	double __ieee754_yn(n,x)
-	int n; double x;
-#endif
+double
+__ieee754_yn(int n, double x)
 {
 	int32_t i,hx,ix,lx;
 	int32_t sign;
@@ -235,9 +227,11 @@ static double zero  =  0.00000000000000000000e+00;
 	EXTRACT_WORDS(hx,lx,x);
 	ix = 0x7fffffff&hx;
     /* if Y(n,NaN) is NaN */
-	if((ix|((u_int32_t)(lx|-lx))>>31)>0x7ff00000) return x+x;
-	if((ix|lx)==0) return -HUGE_VAL+x; /* -inf and overflow exception.  */;
-	if(hx<0) return zero/(zero*x);
+	if(__builtin_expect((ix|((u_int32_t)(lx|-lx))>>31)>0x7ff00000,0))
+		return x+x;
+	if(__builtin_expect((ix|lx)==0, 0))
+		return -HUGE_VAL+x; /* -inf and overflow exception.  */;
+	if(__builtin_expect(hx<0, 0)) return zero/(zero*x);
 	sign = 1;
 	if(n<0){
 		n = -n;
@@ -245,7 +239,7 @@ static double zero  =  0.00000000000000000000e+00;
 	}
 	if(n==0) return(__ieee754_y0(x));
 	if(n==1) return(sign*__ieee754_y1(x));
-	if(ix==0x7ff00000) return zero;
+	if(__builtin_expect(ix==0x7ff00000, 0)) return zero;
 	if(ix>=0x52D00000) { /* x > 2**302 */
     /* (x >> n**2)
      *	    Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
@@ -256,7 +250,7 @@ static double zero  =  0.00000000000000000000e+00;
      *		   n	sin(xn)*sqt2	cos(xn)*sqt2
      *		----------------------------------
      *		   0	 s-c		 c+s
-     *		   1	-s-c 		-c+s
+     *		   1	-s-c		-c+s
      *		   2	-s+c		-c-s
      *		   3	 s+c		 c-s
      */
@@ -285,3 +279,4 @@ static double zero  =  0.00000000000000000000e+00;
 	}
 	if(sign>0) return b; else return -b;
 }
+strong_alias (__ieee754_yn, __yn_finite)
diff --git a/libc/sysdeps/ieee754/dbl-64/e_lgamma_r.c b/libc/sysdeps/ieee754/dbl-64/e_lgamma_r.c
index a298a5a2a..e26ce8a24 100644
--- a/libc/sysdeps/ieee754/dbl-64/e_lgamma_r.c
+++ b/libc/sysdeps/ieee754/dbl-64/e_lgamma_r.c
@@ -10,19 +10,15 @@
  * ====================================================
  */
 
-#if defined(LIBM_SCCS) && !defined(lint)
-static char rcsid[] = "$NetBSD: e_lgamma_r.c,v 1.7 1995/05/10 20:45:42 jtc Exp $";
-#endif
-
 /* __ieee754_lgamma_r(x, signgamp)
  * Reentrant version of the logarithm of the Gamma function
  * with user provide pointer for the sign of Gamma(x).
  *
  * Method:
  *   1. Argument Reduction for 0 < x <= 8
- * 	Since gamma(1+s)=s*gamma(s), for x in [0,8], we may
- * 	reduce x to a number in [1.5,2.5] by
- * 		lgamma(1+s) = log(s) + lgamma(s)
+ *	Since gamma(1+s)=s*gamma(s), for x in [0,8], we may
+ *	reduce x to a number in [1.5,2.5] by
+ *		lgamma(1+s) = log(s) + lgamma(s)
  *	for example,
  *		lgamma(7.3) = log(6.3) + lgamma(6.3)
  *			    = log(6.3*5.3) + lgamma(5.3)
@@ -56,15 +52,15 @@ static char rcsid[] = "$NetBSD: e_lgamma_r.c,v 1.7 1995/05/10 20:45:42 jtc Exp $
  *	Let z = 1/x, then we approximation
  *		f(z) = lgamma(x) - (x-0.5)(log(x)-1)
  *	by
- *	  			    3       5             11
+ *				    3       5             11
  *		w = w0 + w1*z + w2*z  + w3*z  + ... + w6*z
  *	where
  *		|w - f(z)| < 2**-58.74
  *
  *   4. For negative x, since (G is gamma function)
  *		-x*G(-x)*G(x) = pi/sin(pi*x),
- * 	we have
- * 		G(x) = pi/(sin(pi*x)*(-x)*G(-x))
+ *	we have
+ *		G(x) = pi/(sin(pi*x)*(-x)*G(-x))
  *	since G(-x) is positive, sign(G(x)) = sign(sin(pi*x)) for x<0
  *	Hence, for x<0, signgam = sign(sin(pi*x)) and
  *		lgamma(x) = log(|Gamma(x)|)
@@ -77,18 +73,14 @@ static char rcsid[] = "$NetBSD: e_lgamma_r.c,v 1.7 1995/05/10 20:45:42 jtc Exp $
  *		lgamma(1)=lgamma(2)=0
  *		lgamma(x) ~ -log(x) for tiny x
  *		lgamma(0) = lgamma(inf) = inf
- *	 	lgamma(-integer) = +-inf
+ *		lgamma(-integer) = +-inf
  *
  */
 
 #include "math.h"
 #include "math_private.h"
 
-#ifdef __STDC__
 static const double
-#else
-static double
-#endif
 two52=  4.50359962737049600000e+15, /* 0x43300000, 0x00000000 */
 half=  5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */
 one =  1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
@@ -156,18 +148,10 @@ w4  = -5.95187557450339963135e-04, /* 0xBF4380CB, 0x8C0FE741 */
 w5  =  8.36339918996282139126e-04, /* 0x3F4B67BA, 0x4CDAD5D1 */
 w6  = -1.63092934096575273989e-03; /* 0xBF5AB89D, 0x0B9E43E4 */
 
-#ifdef __STDC__
 static const double zero=  0.00000000000000000000e+00;
-#else
-static double zero=  0.00000000000000000000e+00;
-#endif
 
-#ifdef __STDC__
-	static double sin_pi(double x)
-#else
-	static double sin_pi(x)
-	double x;
-#endif
+static double
+sin_pi(double x)
 {
 	double y,z;
 	int n,ix;
@@ -188,16 +172,16 @@ static double zero=  0.00000000000000000000e+00;
 	    y   = 2.0*(y - __floor(y));		/* y = |x| mod 2.0 */
 	    n   = (int) (y*4.0);
 	} else {
-            if(ix>=0x43400000) {
-                y = zero; n = 0;                 /* y must be even */
-            } else {
-                if(ix<0x43300000) z = y+two52;	/* exact */
+	    if(ix>=0x43400000) {
+		y = zero; n = 0;                 /* y must be even */
+	    } else {
+		if(ix<0x43300000) z = y+two52;	/* exact */
 		GET_LOW_WORD(n,z);
 		n &= 1;
-                y  = n;
-                n<<= 2;
-            }
-        }
+		y  = n;
+		n<<= 2;
+	    }
+	}
 	switch (n) {
 	    case 0:   y =  __sin(pi*y); break;
 	    case 1:
@@ -212,12 +196,8 @@ static double zero=  0.00000000000000000000e+00;
 }
 
 
-#ifdef __STDC__
-	double __ieee754_lgamma_r(double x, int *signgamp)
-#else
-	double __ieee754_lgamma_r(x,signgamp)
-	double x; int *signgamp;
-#endif
+double
+__ieee754_lgamma_r(double x, int *signgamp)
 {
 	double t,y,z,nadj,p,p1,p2,p3,q,r,w;
 	int i,hx,lx,ix;
@@ -227,21 +207,23 @@ static double zero=  0.00000000000000000000e+00;
     /* purge off +-inf, NaN, +-0, and negative arguments */
 	*signgamp = 1;
 	ix = hx&0x7fffffff;
-	if(ix>=0x7ff00000) return x*x;
-	if((ix|lx)==0)
+	if(__builtin_expect(ix>=0x7ff00000, 0)) return x*x;
+	if(__builtin_expect((ix|lx)==0, 0))
 	  {
 	    if (hx < 0)
 	      *signgamp = -1;
 	    return one/fabs(x);
 	  }
-	if(ix<0x3b900000) {	/* |x|<2**-70, return -log(|x|) */
+	if(__builtin_expect(ix<0x3b900000, 0)) {
+	    /* |x|<2**-70, return -log(|x|) */
 	    if(hx<0) {
-	        *signgamp = -1;
-	        return -__ieee754_log(-x);
+		*signgamp = -1;
+		return -__ieee754_log(-x);
 	    } else return -__ieee754_log(x);
 	}
 	if(hx<0) {
-	    if(ix>=0x43300000) 	/* |x|>=2**52, must be -integer */
+	    if(__builtin_expect(ix>=0x43300000, 0))
+		/* |x|>=2**52, must be -integer */
 		return x/zero;
 	    t = sin_pi(x);
 	    if(t==zero) return one/fabsf(t); /* -integer */
@@ -254,15 +236,15 @@ static double zero=  0.00000000000000000000e+00;
 	if((((ix-0x3ff00000)|lx)==0)||(((ix-0x40000000)|lx)==0)) r = 0;
     /* for x < 2.0 */
 	else if(ix<0x40000000) {
-	    if(ix<=0x3feccccc) { 	/* lgamma(x) = lgamma(x+1)-log(x) */
+	    if(ix<=0x3feccccc) {	/* lgamma(x) = lgamma(x+1)-log(x) */
 		r = -__ieee754_log(x);
 		if(ix>=0x3FE76944) {y = one-x; i= 0;}
 		else if(ix>=0x3FCDA661) {y= x-(tc-one); i=1;}
-	  	else {y = x; i=2;}
+		else {y = x; i=2;}
 	    } else {
-	  	r = zero;
-	        if(ix>=0x3FFBB4C3) {y=2.0-x;i=0;} /* [1.7316,2] */
-	        else if(ix>=0x3FF3B4C4) {y=x-tc;i=1;} /* [1.23,1.73] */
+		r = zero;
+		if(ix>=0x3FFBB4C3) {y=2.0-x;i=0;} /* [1.7316,2] */
+		else if(ix>=0x3FF3B4C4) {y=x-tc;i=1;} /* [1.23,1.73] */
 		else {y=x-one;i=2;}
 	    }
 	    switch(i) {
@@ -286,7 +268,7 @@ static double zero=  0.00000000000000000000e+00;
 		r += (-0.5*y + p1/p2);
 	    }
 	}
-	else if(ix<0x40200000) { 			/* x < 8.0 */
+	else if(ix<0x40200000) {			/* x < 8.0 */
 	    i = (int)x;
 	    t = zero;
 	    y = x-(double)i;
@@ -315,3 +297,4 @@ static double zero=  0.00000000000000000000e+00;
 	if(hx<0) r = nadj - r;
 	return r;
 }
+strong_alias (__ieee754_lgamma_r, __lgamma_r_finite)
diff --git a/libc/sysdeps/ieee754/dbl-64/e_log.c b/libc/sysdeps/ieee754/dbl-64/e_log.c
index 1a9967b54..14851638a 100644
--- a/libc/sysdeps/ieee754/dbl-64/e_log.c
+++ b/libc/sysdeps/ieee754/dbl-64/e_log.c
@@ -1,7 +1,7 @@
 /*
  * IBM Accurate Mathematical Library
  * written by International Business Machines Corp.
- * Copyright (C) 2001 Free Software Foundation
+ * Copyright (C) 2001, 2011 Free Software Foundation
  *
  * This program is free software; you can redistribute it and/or modify
  * it under the terms of the GNU Lesser General Public License as published by
@@ -36,7 +36,7 @@
 
 
 #include "endian.h"
-#include "dla.h"
+#include <dla.h>
 #include "mpa.h"
 #include "MathLib.h"
 #include "math_private.h"
@@ -55,9 +55,12 @@ double __ieee754_log(double x) {
   int k;
 #endif
   double dbl_n,u,p0,q,r0,w,nln2a,luai,lubi,lvaj,lvbj,
-         sij,ssij,ttij,A,B,B0,y,y1,y2,polI,polII,sa,sb,
-         t1,t2,t3,t4,t5,t6,t7,t8,t,ra,rb,ww,
-         a0,aa0,s1,s2,ss2,s3,ss3,a1,aa1,a,aa,b,bb,c;
+	 sij,ssij,ttij,A,B,B0,y,y1,y2,polI,polII,sa,sb,
+	 t1,t2,t7,t8,t,ra,rb,ww,
+	 a0,aa0,s1,s2,ss2,s3,ss3,a1,aa1,a,aa,b,bb,c;
+#ifndef DLA_FMA
+  double t3,t4,t5,t6;
+#endif
   number num;
   mp_no mpx,mpy,mpy1,mpy2,mperr;
 
@@ -68,18 +71,21 @@ double __ieee754_log(double x) {
 
   num.d = x;  ux = num.i[HIGH_HALF];  dx = num.i[LOW_HALF];
   n=0;
-  if (ux < 0x00100000) {
-    if (((ux & 0x7fffffff) | dx) == 0)  return MHALF/ZERO; /* return -INF */
-    if (ux < 0) return (x-x)/ZERO;                         /* return NaN  */
+  if (__builtin_expect(ux < 0x00100000, 0)) {
+    if (__builtin_expect(((ux & 0x7fffffff) | dx) == 0, 0))
+      return MHALF/ZERO; /* return -INF */
+    if (__builtin_expect(ux < 0, 0))
+      return (x-x)/ZERO;                         /* return NaN  */
     n -= 54;    x *= two54.d;                              /* scale x     */
     num.d = x;
   }
-  if (ux >= 0x7ff00000) return x+x;                        /* INF or NaN  */
+  if (__builtin_expect(ux >= 0x7ff00000, 0))
+    return x+x;                        /* INF or NaN  */
 
   /* Regular values of x */
 
   w = x-ONE;
-  if (ABS(w) > U03) { goto case_03; }
+  if (__builtin_expect(ABS(w) > U03, 1)) { goto case_03; }
 
 
   /*--- Stage I, the case abs(x-1) < 0.03 */
@@ -90,7 +96,7 @@ double __ieee754_log(double x) {
 
   /* Evaluate polynomial II */
   polII = (b0.d+w*(b1.d+w*(b2.d+w*(b3.d+w*(b4.d+
-          w*(b5.d+w*(b6.d+w*(b7.d+w*b8.d))))))))*w*w*w;
+	  w*(b5.d+w*(b6.d+w*(b7.d+w*b8.d))))))))*w*w*w;
   c = (aa+bb)+polII;
 
   /* End stage I, case abs(x-1) < 0.03 */
@@ -99,7 +105,7 @@ double __ieee754_log(double x) {
   /*--- Stage II, the case abs(x-1) < 0.03 */
 
   a = d11.d+w*(d12.d+w*(d13.d+w*(d14.d+w*(d15.d+w*(d16.d+
-            w*(d17.d+w*(d18.d+w*(d19.d+w*d20.d))))))));
+	    w*(d17.d+w*(d18.d+w*(d19.d+w*d20.d))))))));
   EMULV(w,a,s2,ss2,t1,t2,t3,t4,t5)
   ADD2(d10.d,dd10.d,s2,ss2,s3,ss3,t1,t2)
   MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
@@ -201,3 +207,4 @@ double __ieee754_log(double x) {
   }
   return y1;
 }
+strong_alias (__ieee754_log, __log_finite)
diff --git a/libc/sysdeps/ieee754/dbl-64/e_log10.c b/libc/sysdeps/ieee754/dbl-64/e_log10.c
index e8a3278ea..6a630bcef 100644
--- a/libc/sysdeps/ieee754/dbl-64/e_log10.c
+++ b/libc/sysdeps/ieee754/dbl-64/e_log10.c
@@ -10,10 +10,6 @@
  * ====================================================
  */
 
-#if defined(LIBM_SCCS) && !defined(lint)
-static char rcsid[] = "$NetBSD: e_log10.c,v 1.9 1995/05/10 20:45:51 jtc Exp $";
-#endif
-
 /* __ieee754_log10(x)
  * Return the base 10 logarithm of x
  *
@@ -50,28 +46,16 @@ static char rcsid[] = "$NetBSD: e_log10.c,v 1.9 1995/05/10 20:45:51 jtc Exp $";
 #include "math.h"
 #include "math_private.h"
 
-#ifdef __STDC__
 static const double
-#else
-static double
-#endif
 two54      =  1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */
 ivln10     =  4.34294481903251816668e-01, /* 0x3FDBCB7B, 0x1526E50E */
 log10_2hi  =  3.01029995663611771306e-01, /* 0x3FD34413, 0x509F6000 */
 log10_2lo  =  3.69423907715893078616e-13; /* 0x3D59FEF3, 0x11F12B36 */
 
-#ifdef __STDC__
 static const double zero   =  0.0;
-#else
-static double zero   =  0.0;
-#endif
 
-#ifdef __STDC__
-	double __ieee754_log10(double x)
-#else
-	double __ieee754_log10(x)
-	double x;
-#endif
+double
+__ieee754_log10(double x)
 {
 	double y,z;
 	int32_t i,k,hx;
@@ -79,20 +63,22 @@ static double zero   =  0.0;
 
 	EXTRACT_WORDS(hx,lx,x);
 
-        k=0;
-        if (hx < 0x00100000) {			/* x < 2**-1022  */
-            if (((hx&0x7fffffff)|lx)==0)
-                return -two54/(x-x);		/* log(+-0)=-inf */
-            if (hx<0) return (x-x)/(x-x);	/* log(-#) = NaN */
-            k -= 54; x *= two54; /* subnormal number, scale up x */
+	k=0;
+	if (hx < 0x00100000) {			/* x < 2**-1022  */
+	    if (__builtin_expect(((hx&0x7fffffff)|lx)==0, 0))
+		return -two54/(x-x);		/* log(+-0)=-inf */
+	    if (__builtin_expect(hx<0, 0))
+		return (x-x)/(x-x);	/* log(-#) = NaN */
+	    k -= 54; x *= two54; /* subnormal number, scale up x */
 	    GET_HIGH_WORD(hx,x);
-        }
-	if (hx >= 0x7ff00000) return x+x;
+	}
+	if (__builtin_expect(hx >= 0x7ff00000, 0)) return x+x;
 	k += (hx>>20)-1023;
 	i  = ((u_int32_t)k&0x80000000)>>31;
-        hx = (hx&0x000fffff)|((0x3ff-i)<<20);
-        y  = (double)(k+i);
+	hx = (hx&0x000fffff)|((0x3ff-i)<<20);
+	y  = (double)(k+i);
 	SET_HIGH_WORD(x,hx);
 	z  = y*log10_2lo + ivln10*__ieee754_log(x);
 	return  z+y*log10_2hi;
 }
+strong_alias (__ieee754_log10, __log10_finite)
diff --git a/libc/sysdeps/ieee754/dbl-64/e_log2.c b/libc/sysdeps/ieee754/dbl-64/e_log2.c
index f05d0ce96..be41cb4e6 100644
--- a/libc/sysdeps/ieee754/dbl-64/e_log2.c
+++ b/libc/sysdeps/ieee754/dbl-64/e_log2.c
@@ -21,14 +21,14 @@
  *   2. Approximation of log(1+f).
  *	Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s)
  *		 = 2s + 2/3 s**3 + 2/5 s**5 + .....,
- *	     	 = 2s + s*R
+ *		 = 2s + s*R
  *      We use a special Reme algorithm on [0,0.1716] to generate
- * 	a polynomial of degree 14 to approximate R The maximum error
+ *	a polynomial of degree 14 to approximate R The maximum error
  *	of this polynomial approximation is bounded by 2**-58.45. In
  *	other words,
- *		        2      4      6      8      10      12      14
+ *			2      4      6      8      10      12      14
  *	    R(z) ~ Lg1*s +Lg2*s +Lg3*s +Lg4*s +Lg5*s  +Lg6*s  +Lg7*s
- *  	(the values of Lg1 to Lg7 are listed in the program)
+ *	(the values of Lg1 to Lg7 are listed in the program)
  *	and
  *	    |      2          14          |     -58.45
  *	    | Lg1*s +...+Lg7*s    -  R(z) | <= 2
@@ -57,11 +57,7 @@
 #include "math.h"
 #include "math_private.h"
 
-#ifdef __STDC__
 static const double
-#else
-static double
-#endif
 ln2 = 0.69314718055994530942,
 two54   =  1.80143985094819840000e+16,  /* 43500000 00000000 */
 Lg1 = 6.666666666666735130e-01,  /* 3FE55555 55555593 */
@@ -72,18 +68,10 @@ Lg5 = 1.818357216161805012e-01,  /* 3FC74664 96CB03DE */
 Lg6 = 1.531383769920937332e-01,  /* 3FC39A09 D078C69F */
 Lg7 = 1.479819860511658591e-01;  /* 3FC2F112 DF3E5244 */
 
-#ifdef __STDC__
 static const double zero   =  0.0;
-#else
-static double zero   =  0.0;
-#endif
 
-#ifdef __STDC__
-	double __ieee754_log2(double x)
-#else
-	double __ieee754_log2(x)
-	double x;
-#endif
+double
+__ieee754_log2(double x)
 {
 	double hfsq,f,s,z,R,w,t1,t2,dk;
 	int32_t k,hx,i,j;
@@ -93,13 +81,14 @@ static double zero   =  0.0;
 
 	k=0;
 	if (hx < 0x00100000) {			/* x < 2**-1022  */
-	    if (((hx&0x7fffffff)|lx)==0)
+	    if (__builtin_expect(((hx&0x7fffffff)|lx)==0, 0))
 		return -two54/(x-x);		/* log(+-0)=-inf */
-	    if (hx<0) return (x-x)/(x-x);	/* log(-#) = NaN */
+	    if (__builtin_expect(hx<0, 0))
+		return (x-x)/(x-x);	/* log(-#) = NaN */
 	    k -= 54; x *= two54; /* subnormal number, scale up x */
 	    GET_HIGH_WORD(hx,x);
 	}
-	if (hx >= 0x7ff00000) return x+x;
+	if (__builtin_expect(hx >= 0x7ff00000, 0)) return x+x;
 	k += (hx>>20)-1023;
 	hx &= 0x000fffff;
 	i = (hx+0x95f64)&0x100000;
@@ -112,7 +101,7 @@ static double zero   =  0.0;
 	    R = f*f*(0.5-0.33333333333333333*f);
 	    return dk-(R-f)/ln2;
 	}
- 	s = f/(2.0+f);
+	s = f/(2.0+f);
 	z = s*s;
 	i = hx-0x6147a;
 	w = z*z;
@@ -128,3 +117,4 @@ static double zero   =  0.0;
 	    return dk-((s*(f-R))-f)/ln2;
 	}
 }
+strong_alias (__ieee754_log2, __log2_finite)
diff --git a/libc/sysdeps/ieee754/dbl-64/e_pow.c b/libc/sysdeps/ieee754/dbl-64/e_pow.c
index 1e159f2c0..789054015 100644
--- a/libc/sysdeps/ieee754/dbl-64/e_pow.c
+++ b/libc/sysdeps/ieee754/dbl-64/e_pow.c
@@ -1,7 +1,7 @@
 /*
  * IBM Accurate Mathematical Library
  * written by International Business Machines Corp.
- * Copyright (C) 2001, 2002, 2004 Free Software Foundation
+ * Copyright (C) 2001, 2002, 2004, 2011 Free Software Foundation
  *
  * This program is free software; you can redistribute it and/or modify
  * it under the terms of the GNU Lesser General Public License as published by
@@ -22,12 +22,12 @@
 /*                                                                         */
 /*  FUNCTIONS: upow                                                        */
 /*             power1                                                      */
-/*             my_log2                                                        */
+/*             my_log2                                                     */
 /*             log1                                                        */
 /*             checkint                                                    */
 /* FILES NEEDED: dla.h endian.h mpa.h mydefs.h                             */
 /*               halfulp.c mpexp.c mplog.c slowexp.c slowpow.c mpa.c       */
-/*                          uexp.c  upow.c			           */
+/*                          uexp.c  upow.c				   */
 /*               root.tbl uexp.tbl upow.tbl                                */
 /* An ultimate power routine. Given two IEEE double machine numbers y,x    */
 /* it computes the correctly rounded (to nearest) value of x^y.            */
@@ -37,7 +37,7 @@
 /***************************************************************************/
 #include "endian.h"
 #include "upow.h"
-#include "dla.h"
+#include <dla.h>
 #include "mydefs.h"
 #include "MathLib.h"
 #include "upow.tbl"
@@ -77,7 +77,7 @@ double __ieee754_pow(double x, double y) {
   /* else */
   if(((u.i[HIGH_HALF]>0 && u.i[HIGH_HALF]<0x7ff00000)||        /* x>0 and not x->0 */
        (u.i[HIGH_HALF]==0 && u.i[LOW_HALF]!=0))  &&
-                                      /*   2^-1023< x<= 2^-1023 * 0x1.0000ffffffff */
+				      /*   2^-1023< x<= 2^-1023 * 0x1.0000ffffffff */
       (v.i[HIGH_HALF]&0x7fffffff) < 0x4ff00000) {              /* if y<-1 or y>1   */
     z = log1(x,&aa,&error);                                 /* x^y  =e^(y log (X)) */
     t = y*134217729.0;
@@ -153,6 +153,7 @@ double __ieee754_pow(double x, double y) {
   if (y<0) return (x<1.0)?INF.x:0;
   return 0;     /* unreachable, to make the compiler happy */
 }
+strong_alias (__ieee754_pow, __pow_finite)
 
 /**************************************************************************/
 /* Computing x^y using more accurate but more slow log routine            */
@@ -282,7 +283,10 @@ static double my_log2(double x, double *delta, double *error) {
   double cor;
 #endif
   double ou1,ou2,lu1,lu2,ov,lv1,lv2,a,a1,a2;
-  double y,yy,z,zz,j1,j2,j3,j4,j5,j6,j7,j8;
+  double y,yy,z,zz,j1,j2,j7,j8;
+#ifndef DLA_FMA
+  double j3,j4,j5,j6;
+#endif
   mynumber u,v;
 #ifdef BIG_ENDI
   mynumber
diff --git a/libc/sysdeps/ieee754/dbl-64/e_remainder.c b/libc/sysdeps/ieee754/dbl-64/e_remainder.c
index cc06e18ce..d1782a15c 100644
--- a/libc/sysdeps/ieee754/dbl-64/e_remainder.c
+++ b/libc/sysdeps/ieee754/dbl-64/e_remainder.c
@@ -1,8 +1,8 @@
 /*
  * IBM Accurate Mathematical Library
  * written by International Business Machines Corp.
- * Copyright (C) 2001 Free Software Foundation
- * 
+ * Copyright (C) 2001, 2011 Free Software Foundation
+ *
  * This program is free software; you can redistribute it and/or modify
  * it under the terms of the GNU Lesser General Public License as published by
  * the Free Software Foundation; either version 2.1 of the License, or
@@ -96,9 +96,9 @@ double __ieee754_remainder(double x, double y)
       u.x=(u.x-d*w.x)-d*ww.x;
       if (ABS(u.x)<0.5*t.x) return (u.x!=0)?u.x:((x>0)?ZERO.x:nZERO.x);
       else
-        if (ABS(u.x)>0.5*t.x) return (d>z)?u.x+t.x:u.x-t.x;
-        else
-        {z=u.x/t.x; d=(z+big.x)-big.x; return ((u.x-d*w.x)-d*ww.x);}
+	if (ABS(u.x)>0.5*t.x) return (d>z)?u.x+t.x:u.x-t.x;
+	else
+	{z=u.x/t.x; d=(z+big.x)-big.x; return ((u.x-d*w.x)-d*ww.x);}
     }
 
   }   /*   (kx<0x7fe00000&&ky<0x7ff00000&&ky>=0x03500000)     */
@@ -128,3 +128,4 @@ double __ieee754_remainder(double x, double y)
    }
   }
 }
+strong_alias (__ieee754_remainder, __remainder_finite)
diff --git a/libc/sysdeps/ieee754/dbl-64/e_sinh.c b/libc/sysdeps/ieee754/dbl-64/e_sinh.c
index 1701b9bb6..50463c304 100644
--- a/libc/sysdeps/ieee754/dbl-64/e_sinh.c
+++ b/libc/sysdeps/ieee754/dbl-64/e_sinh.c
@@ -5,7 +5,7 @@
  *
  * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  */
@@ -15,15 +15,15 @@ static char rcsid[] = "$NetBSD: e_sinh.c,v 1.7 1995/05/10 20:46:13 jtc Exp $";
 #endif
 
 /* __ieee754_sinh(x)
- * Method : 
+ * Method :
  * mathematically sinh(x) if defined to be (exp(x)-exp(-x))/2
- *	1. Replace x by |x| (sinh(-x) = -sinh(x)). 
- *	2. 
- *		                                    E + E/(E+1)
+ *	1. Replace x by |x| (sinh(-x) = -sinh(x)).
+ *	2.
+ *						    E + E/(E+1)
  *	    0        <= x <= 22     :  sinh(x) := --------------, E=expm1(x)
- *			       			        2
+ *							2
  *
- *	    22       <= x <= lnovft :  sinh(x) := exp(x)/2 
+ *	    22       <= x <= lnovft :  sinh(x) := exp(x)/2
  *	    lnovft   <= x <= ln2ovft:  sinh(x) := exp(x/2)/2 * exp(x/2)
  *	    ln2ovft  <  x	    :  sinh(x) := x*shuge (overflow)
  *
@@ -35,19 +35,11 @@ static char rcsid[] = "$NetBSD: e_sinh.c,v 1.7 1995/05/10 20:46:13 jtc Exp $";
 #include "math.h"
 #include "math_private.h"
 
-#ifdef __STDC__
 static const double one = 1.0, shuge = 1.0e307;
-#else
-static double one = 1.0, shuge = 1.0e307;
-#endif
 
-#ifdef __STDC__
-	double __ieee754_sinh(double x)
-#else
-	double __ieee754_sinh(x)
-	double x;
-#endif
-{	
+double
+__ieee754_sinh(double x)
+{
 	double t,w,h;
 	int32_t ix,jx;
 	u_int32_t lx;
@@ -57,14 +49,15 @@ static double one = 1.0, shuge = 1.0e307;
 	ix = jx&0x7fffffff;
 
     /* x is INF or NaN */
-	if(ix>=0x7ff00000) return x+x;	
+	if(__builtin_expect(ix>=0x7ff00000, 0)) return x+x;
 
 	h = 0.5;
 	if (jx<0) h = -h;
     /* |x| in [0,22], return sign(x)*0.5*(E+E/(E+1))) */
 	if (ix < 0x40360000) {		/* |x|<22 */
-	    if (ix<0x3e300000) 		/* |x|<2**-28 */
-		if(shuge+x>one) return x;/* sinh(tiny) = tiny with inexact */
+	    if (__builtin_expect(ix<0x3e300000, 0))	/* |x|<2**-28 */
+		if(shuge+x>one)
+		    return x;/* sinh(tiny) = tiny with inexact */
 	    t = __expm1(fabs(x));
 	    if(ix<0x3ff00000) return h*(2.0*t-t*t/(t+one));
 	    return h*(t+t/(t+one));
@@ -84,3 +77,4 @@ static double one = 1.0, shuge = 1.0e307;
     /* |x| > overflowthresold, sinh(x) overflow */
 	return x*shuge;
 }
+strong_alias (__ieee754_sinh, __sinh_finite)
diff --git a/libc/sysdeps/ieee754/dbl-64/e_sqrt.c b/libc/sysdeps/ieee754/dbl-64/e_sqrt.c
index f7e805549..c507c598d 100644
--- a/libc/sysdeps/ieee754/dbl-64/e_sqrt.c
+++ b/libc/sysdeps/ieee754/dbl-64/e_sqrt.c
@@ -1,7 +1,7 @@
 /*
  * IBM Accurate Mathematical Library
  * written by International Business Machines Corp.
- * Copyright (C) 2001 Free Software Foundation
+ * Copyright (C) 2001, 2011 Free Software Foundation
  *
  * This program is free software; you can redistribute it and/or modify
  * it under the terms of the GNU Lesser General Public License as published by
@@ -35,7 +35,7 @@
 
 #include "endian.h"
 #include "mydefs.h"
-#include "dla.h"
+#include <dla.h>
 #include "MathLib.h"
 #include "root.tbl"
 #include "math_private.h"
@@ -86,3 +86,4 @@ double __ieee754_sqrt(double x) {
     return tm256.x*__ieee754_sqrt(x*t512.x);
   }
 }
+strong_alias (__ieee754_sqrt, __sqrt_finite)
diff --git a/libc/sysdeps/ieee754/dbl-64/halfulp.c b/libc/sysdeps/ieee754/dbl-64/halfulp.c
index 478a4bacf..373d40522 100644
--- a/libc/sysdeps/ieee754/dbl-64/halfulp.c
+++ b/libc/sysdeps/ieee754/dbl-64/halfulp.c
@@ -1,7 +1,7 @@
 /*
  * IBM Accurate Mathematical Library
  * written by International Business Machines Corp.
- * Copyright (C) 2001, 2005 Free Software Foundation
+ * Copyright (C) 2001, 2005, 2011 Free Software Foundation
  *
  * This program is free software; you can redistribute it and/or modify
  * it under the terms of the GNU Lesser General Public License as published by
@@ -37,22 +37,23 @@
 
 #include "endian.h"
 #include "mydefs.h"
-#include "dla.h"
+#include <dla.h>
 #include "math_private.h"
 
-double __ieee754_sqrt(double x);
-
 static const int4 tab54[32] = {
    262143, 11585, 1782, 511, 210, 107, 63, 42,
        30,    22,   17,  14,  12,  10,  9,  7,
-        7,     6,    5,   5,   5,   4,  4,  4,
-        3,     3,    3,   3,   3,   3,  3,  3 };
+	7,     6,    5,   5,   5,   4,  4,  4,
+	3,     3,    3,   3,   3,   3,  3,  3 };
 
 
 double __halfulp(double x, double y)
 {
   mynumber v;
-  double z,u,uu,j1,j2,j3,j4,j5;
+  double z,u,uu;
+#ifndef DLA_FMA
+  double j1,j2,j3,j4,j5;
+#endif
   int4 k,l,m,n;
   if (y <= 0) {               /*if power is negative or zero */
     v.x = y;
@@ -64,12 +65,12 @@ double __halfulp(double x, double y)
     z = (double) k;
     return (z*y == -1075.0)?0: -10.0;
   }
-                              /* if y > 0  */
+			      /* 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;
@@ -90,7 +91,7 @@ double __halfulp(double x, double y)
   k = -k;
   if (k>5) return -10.0;
 
-                            /*   now treat x        */
+			    /*   now treat x        */
   while (k>0) {
     z = __ieee754_sqrt(x);
     EMULV(z,z,u,uu,j1,j2,j3,j4,j5);
@@ -111,11 +112,11 @@ double __halfulp(double x, double y)
   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 */
+	    /*   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  */
+	     /* yes, it is - now compute x**n by simple multiplications  */
 
   u = x;
   for (k=1;k<n;k++) u = u*x;
diff --git a/libc/sysdeps/ieee754/dbl-64/s_asinh.c b/libc/sysdeps/ieee754/dbl-64/s_asinh.c
index 985cfe32e..93789fb41 100644
--- a/libc/sysdeps/ieee754/dbl-64/s_asinh.c
+++ b/libc/sysdeps/ieee754/dbl-64/s_asinh.c
@@ -10,10 +10,6 @@
  * ====================================================
  */
 
-#if defined(LIBM_SCCS) && !defined(lint)
-static char rcsid[] = "$NetBSD: s_asinh.c,v 1.9 1995/05/12 04:57:37 jtc Exp $";
-#endif
-
 /* asinh(x)
  * Method :
  *	Based on
@@ -28,40 +24,34 @@ static char rcsid[] = "$NetBSD: s_asinh.c,v 1.9 1995/05/12 04:57:37 jtc Exp $";
 #include "math.h"
 #include "math_private.h"
 
-#ifdef __STDC__
 static const double
-#else
-static double
-#endif
 one =  1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
 ln2 =  6.93147180559945286227e-01, /* 0x3FE62E42, 0xFEFA39EF */
 huge=  1.00000000000000000000e+300;
 
-#ifdef __STDC__
-	double __asinh(double x)
-#else
-	double __asinh(x)
-	double x;
-#endif
+double
+__asinh(double x)
 {
-	double t,w;
+	double w;
 	int32_t hx,ix;
 	GET_HIGH_WORD(hx,x);
 	ix = hx&0x7fffffff;
-	if(ix>=0x7ff00000) return x+x;	/* x is inf or NaN */
-	if(ix< 0x3e300000) {	/* |x|<2**-28 */
+	if(__builtin_expect(ix< 0x3e300000, 0)) {	/* |x|<2**-28 */
 	    if(huge+x>one) return x;	/* return x inexact except 0 */
 	}
-	if(ix>0x41b00000) {	/* |x| > 2**28 */
+	if(__builtin_expect(ix>0x41b00000, 0)) {	/* |x| > 2**28 */
+	    if(ix>=0x7ff00000) return x+x;	/* x is inf or NaN */
 	    w = __ieee754_log(fabs(x))+ln2;
-	} else if (ix>0x40000000) {	/* 2**28 > |x| > 2.0 */
-	    t = fabs(x);
-	    w = __ieee754_log(2.0*t+one/(__ieee754_sqrt(x*x+one)+t));
-	} else {		/* 2.0 > |x| > 2**-28 */
-	    t = x*x;
-	    w =__log1p(fabs(x)+t/(one+__ieee754_sqrt(one+t)));
+	} else {
+	    double xa = fabs(x);
+	    if (ix>0x40000000) {	/* 2**28 > |x| > 2.0 */
+		w = __ieee754_log(2.0*xa+one/(__ieee754_sqrt(xa*xa+one)+xa));
+	    } else {		/* 2.0 > |x| > 2**-28 */
+		double t = xa*xa;
+		w =__log1p(xa+t/(one+__ieee754_sqrt(one+t)));
+	    }
 	}
-	if(hx>0) return w; else return -w;
+	return __copysign(w, x);
 }
 weak_alias (__asinh, asinh)
 #ifdef NO_LONG_DOUBLE
diff --git a/libc/sysdeps/ieee754/dbl-64/s_atan.c b/libc/sysdeps/ieee754/dbl-64/s_atan.c
index 556f5b216..5ea83261a 100644
--- a/libc/sysdeps/ieee754/dbl-64/s_atan.c
+++ b/libc/sysdeps/ieee754/dbl-64/s_atan.c
@@ -1,7 +1,7 @@
 /*
  * IBM Accurate Mathematical Library
  * written by International Business Machines Corp.
- * Copyright (C) 2001 Free Software Foundation
+ * Copyright (C) 2001, 2011 Free Software Foundation
  *
  * This program is free software; you can redistribute it and/or modify
  * it under the terms of the GNU Lesser General Public License as published by
@@ -37,7 +37,7 @@
 /*                                                                      */
 /************************************************************************/
 
-#include "dla.h"
+#include <dla.h>
 #include "mpa.h"
 #include "MathLib.h"
 #include "uatan.tbl"
@@ -52,8 +52,11 @@ double __signArctan(double,double);
 double atan(double x) {
 
 
-  double cor,s1,ss1,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,u,u2,u3,
-         v,vv,w,ww,y,yy,z,zz;
+  double cor,s1,ss1,s2,ss2,t1,t2,t3,t7,t8,t9,t10,u,u2,u3,
+	 v,vv,w,ww,y,yy,z,zz;
+#ifndef DLA_FMA
+  double t4,t5,t6;
+#endif
 #if 0
   double y1,y2;
 #endif
@@ -78,44 +81,44 @@ double atan(double x) {
   if (u<C) {
     if (u<B) {
       if (u<A) {                                           /* u < A */
-         return x; }
+	 return x; }
       else {                                               /* A <= u < B */
-        v=x*x;  yy=x*v*(d3.d+v*(d5.d+v*(d7.d+v*(d9.d+v*(d11.d+v*d13.d)))));
-        if ((y=x+(yy-U1*x)) == x+(yy+U1*x))  return y;
-
-        EMULV(x,x,v,vv,t1,t2,t3,t4,t5)                       /* v+vv=x^2 */
-        s1=v*(f11.d+v*(f13.d+v*(f15.d+v*(f17.d+v*f19.d))));
-        ADD2(f9.d,ff9.d,s1,ZERO,s2,ss2,t1,t2)
-        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
-        ADD2(f7.d,ff7.d,s1,ss1,s2,ss2,t1,t2)
-        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
-        ADD2(f5.d,ff5.d,s1,ss1,s2,ss2,t1,t2)
-        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
-        ADD2(f3.d,ff3.d,s1,ss1,s2,ss2,t1,t2)
-        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
-        MUL2(x,ZERO,s1,ss1,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
-        ADD2(x,ZERO,s2,ss2,s1,ss1,t1,t2)
-        if ((y=s1+(ss1-U5*s1)) == s1+(ss1+U5*s1))  return y;
-
-        return atanMp(x,pr);
+	v=x*x;  yy=x*v*(d3.d+v*(d5.d+v*(d7.d+v*(d9.d+v*(d11.d+v*d13.d)))));
+	if ((y=x+(yy-U1*x)) == x+(yy+U1*x))  return y;
+
+	EMULV(x,x,v,vv,t1,t2,t3,t4,t5)                       /* v+vv=x^2 */
+	s1=v*(f11.d+v*(f13.d+v*(f15.d+v*(f17.d+v*f19.d))));
+	ADD2(f9.d,ff9.d,s1,ZERO,s2,ss2,t1,t2)
+	MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
+	ADD2(f7.d,ff7.d,s1,ss1,s2,ss2,t1,t2)
+	MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
+	ADD2(f5.d,ff5.d,s1,ss1,s2,ss2,t1,t2)
+	MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
+	ADD2(f3.d,ff3.d,s1,ss1,s2,ss2,t1,t2)
+	MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
+	MUL2(x,ZERO,s1,ss1,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
+	ADD2(x,ZERO,s2,ss2,s1,ss1,t1,t2)
+	if ((y=s1+(ss1-U5*s1)) == s1+(ss1+U5*s1))  return y;
+
+	return atanMp(x,pr);
       } }
     else {  /* B <= u < C */
       i=(TWO52+TWO8*u)-TWO52;  i-=16;
       z=u-cij[i][0].d;
       yy=z*(cij[i][2].d+z*(cij[i][3].d+z*(cij[i][4].d+
-                        z*(cij[i][5].d+z* cij[i][6].d))));
+			z*(cij[i][5].d+z* cij[i][6].d))));
       t1=cij[i][1].d;
       if (i<112) {
-        if (i<48)  u2=U21;    /* u < 1/4        */
-        else       u2=U22; }  /* 1/4 <= u < 1/2 */
+	if (i<48)  u2=U21;    /* u < 1/4        */
+	else       u2=U22; }  /* 1/4 <= u < 1/2 */
       else {
-        if (i<176) u2=U23;    /* 1/2 <= u < 3/4 */
-        else       u2=U24; }  /* 3/4 <= u <= 1  */
+	if (i<176) u2=U23;    /* 1/2 <= u < 3/4 */
+	else       u2=U24; }  /* 3/4 <= u <= 1  */
       if ((y=t1+(yy-u2*t1)) == t1+(yy+u2*t1))  return __signArctan(x,y);
 
       z=u-hij[i][0].d;
       s1=z*(hij[i][11].d+z*(hij[i][12].d+z*(hij[i][13].d+
-         z*(hij[i][14].d+z* hij[i][15].d))));
+	 z*(hij[i][14].d+z* hij[i][15].d))));
       ADD2(hij[i][9].d,hij[i][10].d,s1,ZERO,s2,ss2,t1,t2)
       MUL2(z,ZERO,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
       ADD2(hij[i][7].d,hij[i][8].d,s1,ss1,s2,ss2,t1,t2)
@@ -138,7 +141,7 @@ double atan(double x) {
       i=(TWO52+TWO8*w)-TWO52;  i-=16;
       z=(w-cij[i][0].d)+ww;
       yy=HPI1-z*(cij[i][2].d+z*(cij[i][3].d+z*(cij[i][4].d+
-                             z*(cij[i][5].d+z* cij[i][6].d))));
+			     z*(cij[i][5].d+z* cij[i][6].d))));
       t1=HPI-cij[i][1].d;
       if (i<112)  u3=U31;  /* w <  1/2 */
       else        u3=U32;  /* w >= 1/2 */
@@ -148,7 +151,7 @@ double atan(double x) {
       t1=w-hij[i][0].d;
       EADD(t1,ww,z,zz)
       s1=z*(hij[i][11].d+z*(hij[i][12].d+z*(hij[i][13].d+
-         z*(hij[i][14].d+z* hij[i][15].d))));
+	 z*(hij[i][14].d+z* hij[i][15].d))));
       ADD2(hij[i][9].d,hij[i][10].d,s1,ZERO,s2,ss2,t1,t2)
       MUL2(z,zz,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
       ADD2(hij[i][7].d,hij[i][8].d,s1,ss1,s2,ss2,t1,t2)
@@ -165,36 +168,36 @@ double atan(double x) {
     }
     else {
       if (u<E) { /* D <= u < E */
-        w=ONE/u;   v=w*w;
-        EMULV(w,u,t1,t2,t3,t4,t5,t6,t7)
-        yy=w*v*(d3.d+v*(d5.d+v*(d7.d+v*(d9.d+v*(d11.d+v*d13.d)))));
-        ww=w*((ONE-t1)-t2);
-        ESUB(HPI,w,t3,cor)
-        yy=((HPI1+cor)-ww)-yy;
-        if ((y=t3+(yy-U4)) == t3+(yy+U4))  return __signArctan(x,y);
-
-        DIV2(ONE,ZERO,u,ZERO,w,ww,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
-        MUL2(w,ww,w,ww,v,vv,t1,t2,t3,t4,t5,t6,t7,t8)
-        s1=v*(f11.d+v*(f13.d+v*(f15.d+v*(f17.d+v*f19.d))));
-        ADD2(f9.d,ff9.d,s1,ZERO,s2,ss2,t1,t2)
-        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
-        ADD2(f7.d,ff7.d,s1,ss1,s2,ss2,t1,t2)
-        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
-        ADD2(f5.d,ff5.d,s1,ss1,s2,ss2,t1,t2)
-        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
-        ADD2(f3.d,ff3.d,s1,ss1,s2,ss2,t1,t2)
-        MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
-        MUL2(w,ww,s1,ss1,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
-        ADD2(w,ww,s2,ss2,s1,ss1,t1,t2)
-        SUB2(HPI,HPI1,s1,ss1,s2,ss2,t1,t2)
-        if ((y=s2+(ss2-U8)) == s2+(ss2+U8))  return __signArctan(x,y);
+	w=ONE/u;   v=w*w;
+	EMULV(w,u,t1,t2,t3,t4,t5,t6,t7)
+	yy=w*v*(d3.d+v*(d5.d+v*(d7.d+v*(d9.d+v*(d11.d+v*d13.d)))));
+	ww=w*((ONE-t1)-t2);
+	ESUB(HPI,w,t3,cor)
+	yy=((HPI1+cor)-ww)-yy;
+	if ((y=t3+(yy-U4)) == t3+(yy+U4))  return __signArctan(x,y);
+
+	DIV2(ONE,ZERO,u,ZERO,w,ww,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
+	MUL2(w,ww,w,ww,v,vv,t1,t2,t3,t4,t5,t6,t7,t8)
+	s1=v*(f11.d+v*(f13.d+v*(f15.d+v*(f17.d+v*f19.d))));
+	ADD2(f9.d,ff9.d,s1,ZERO,s2,ss2,t1,t2)
+	MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
+	ADD2(f7.d,ff7.d,s1,ss1,s2,ss2,t1,t2)
+	MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
+	ADD2(f5.d,ff5.d,s1,ss1,s2,ss2,t1,t2)
+	MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
+	ADD2(f3.d,ff3.d,s1,ss1,s2,ss2,t1,t2)
+	MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
+	MUL2(w,ww,s1,ss1,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
+	ADD2(w,ww,s2,ss2,s1,ss1,t1,t2)
+	SUB2(HPI,HPI1,s1,ss1,s2,ss2,t1,t2)
+	if ((y=s2+(ss2-U8)) == s2+(ss2+U8))  return __signArctan(x,y);
 
       return atanMp(x,pr);
       }
       else {
-        /* u >= E */
-        if (x>0) return  HPI;
-        else     return MHPI; }
+	/* u >= E */
+	if (x>0) return  HPI;
+	else     return MHPI; }
     }
   }
 
diff --git a/libc/sysdeps/ieee754/dbl-64/s_ceil.c b/libc/sysdeps/ieee754/dbl-64/s_ceil.c
index 1b352a679..695cae5d5 100644
--- a/libc/sysdeps/ieee754/dbl-64/s_ceil.c
+++ b/libc/sysdeps/ieee754/dbl-64/s_ceil.c
@@ -10,10 +10,6 @@
  * ====================================================
  */
 
-#if defined(LIBM_SCCS) && !defined(lint)
-static char rcsid[] = "$NetBSD: s_ceil.c,v 1.8 1995/05/10 20:46:53 jtc Exp $";
-#endif
-
 /*
  * ceil(x)
  * Return x rounded toward -inf to integral value
@@ -26,18 +22,10 @@ static char rcsid[] = "$NetBSD: s_ceil.c,v 1.8 1995/05/10 20:46:53 jtc Exp $";
 #include "math.h"
 #include "math_private.h"
 
-#ifdef __STDC__
 static const double huge = 1.0e300;
-#else
-static double huge = 1.0e300;
-#endif
 
-#ifdef __STDC__
-	double __ceil(double x)
-#else
-	double __ceil(x)
-	double x;
-#endif
+double
+__ceil(double x)
 {
 	int32_t i0,i1,j0;
 	u_int32_t i,j;
@@ -78,8 +66,10 @@ static double huge = 1.0e300;
 	INSERT_WORDS(x,i0,i1);
 	return x;
 }
+#ifndef __ceil
 weak_alias (__ceil, ceil)
-#ifdef NO_LONG_DOUBLE
+# ifdef NO_LONG_DOUBLE
 strong_alias (__ceil, __ceill)
 weak_alias (__ceil, ceill)
+# endif
 #endif
diff --git a/libc/sysdeps/ieee754/dbl-64/s_finite.c b/libc/sysdeps/ieee754/dbl-64/s_finite.c
index 90de0f9d1..2ca3bf245 100644
--- a/libc/sysdeps/ieee754/dbl-64/s_finite.c
+++ b/libc/sysdeps/ieee754/dbl-64/s_finite.c
@@ -22,6 +22,7 @@ static char rcsid[] = "$NetBSD: s_finite.c,v 1.8 1995/05/10 20:47:17 jtc Exp $";
 #include "math.h"
 #include "math_private.h"
 
+#undef __finite
 #ifdef __STDC__
 	int __finite(double x)
 #else
diff --git a/libc/sysdeps/ieee754/dbl-64/s_floor.c b/libc/sysdeps/ieee754/dbl-64/s_floor.c
index 77db9ef39..5b593ca31 100644
--- a/libc/sysdeps/ieee754/dbl-64/s_floor.c
+++ b/libc/sysdeps/ieee754/dbl-64/s_floor.c
@@ -10,10 +10,6 @@
  * ====================================================
  */
 
-#if defined(LIBM_SCCS) && !defined(lint)
-static char rcsid[] = "$NetBSD: s_floor.c,v 1.8 1995/05/10 20:47:20 jtc Exp $";
-#endif
-
 /*
  * floor(x)
  * Return x rounded toward -inf to integral value
@@ -44,7 +40,7 @@ static double huge = 1.0e300;
 	EXTRACT_WORDS(i0,i1,x);
 	j0 = ((i0>>20)&0x7ff)-0x3ff;
 	if(j0<20) {
-	    if(j0<0) { 	/* raise inexact if x != 0 */
+	    if(j0<0) {	/* raise inexact if x != 0 */
 		if(huge+x>0.0) {/* return 0*sign(x) if |x|<1 */
 		    if(i0>=0) {i0=i1=0;}
 		    else if(((i0&0x7fffffff)|i1)!=0)
@@ -64,12 +60,12 @@ static double huge = 1.0e300;
 	} else {
 	    i = ((u_int32_t)(0xffffffff))>>(j0-20);
 	    if((i1&i)==0) return x;	/* x is integral */
-	    if(huge+x>0.0) { 		/* raise inexact flag */
+	    if(huge+x>0.0) {		/* raise inexact flag */
 		if(i0<0) {
 		    if(j0==20) i0+=1;
 		    else {
 			j = i1+(1<<(52-j0));
-			if(j<i1) i0 +=1 ; 	/* got a carry */
+			if(j<i1) i0 +=1 ;	/* got a carry */
 			i1=j;
 		    }
 		}
@@ -79,8 +75,10 @@ static double huge = 1.0e300;
 	INSERT_WORDS(x,i0,i1);
 	return x;
 }
+#ifndef __floor
 weak_alias (__floor, floor)
-#ifdef NO_LONG_DOUBLE
+# ifdef NO_LONG_DOUBLE
 strong_alias (__floor, __floorl)
 weak_alias (__floor, floorl)
+# endif
 #endif
diff --git a/libc/sysdeps/ieee754/dbl-64/s_fma.c b/libc/sysdeps/ieee754/dbl-64/s_fma.c
index 3b0bfd5ce..14c6503de 100644
--- a/libc/sysdeps/ieee754/dbl-64/s_fma.c
+++ b/libc/sysdeps/ieee754/dbl-64/s_fma.c
@@ -1,5 +1,5 @@
 /* Compute x * y + z as ternary operation.
-   Copyright (C) 2010 Free Software Foundation, Inc.
+   Copyright (C) 2010, 2011 Free Software Foundation, Inc.
    This file is part of the GNU C Library.
    Contributed by Jakub Jelinek <jakub@redhat.com>, 2010.
 
@@ -22,6 +22,7 @@
 #include <math.h>
 #include <fenv.h>
 #include <ieee754.h>
+#include <math_private.h>
 
 /* This implementation uses rounding to odd to avoid problems with
    double rounding.  See a paper by Boldo and Melquiond:
@@ -47,7 +48,7 @@ __fma (double x, double y, double z)
 	 z rather than NaN.  */
       if (w.ieee.exponent == 0x7ff
 	  && u.ieee.exponent != 0x7ff
-          && v.ieee.exponent != 0x7ff)
+	  && v.ieee.exponent != 0x7ff)
 	return (z + x) + y;
       /* If x or y or z is Inf/NaN, or if fma will certainly overflow,
 	 or if x * y is less than half of DBL_DENORM_MIN,
@@ -148,34 +149,33 @@ __fma (double x, double y, double z)
   double a2 = t1 + t2;
 
   fenv_t env;
-  feholdexcept (&env);
-  fesetround (FE_TOWARDZERO);
+  libc_feholdexcept_setround (&env, FE_TOWARDZERO);
   /* Perform m2 + a2 addition with round to odd.  */
   u.d = a2 + m2;
 
   if (__builtin_expect (adjust == 0, 1))
     {
       if ((u.ieee.mantissa1 & 1) == 0 && u.ieee.exponent != 0x7ff)
-	u.ieee.mantissa1 |= fetestexcept (FE_INEXACT) != 0;
-      feupdateenv (&env);
+	u.ieee.mantissa1 |= libc_fetestexcept (FE_INEXACT) != 0;
+      libc_feupdateenv (&env);
       /* Result is a1 + u.d.  */
       return a1 + u.d;
     }
   else if (__builtin_expect (adjust > 0, 1))
     {
       if ((u.ieee.mantissa1 & 1) == 0 && u.ieee.exponent != 0x7ff)
-	u.ieee.mantissa1 |= fetestexcept (FE_INEXACT) != 0;
-      feupdateenv (&env);
+	u.ieee.mantissa1 |= libc_fetestexcept (FE_INEXACT) != 0;
+      libc_feupdateenv (&env);
       /* Result is a1 + u.d, scaled up.  */
       return (a1 + u.d) * 0x1p53;
     }
   else
     {
       if ((u.ieee.mantissa1 & 1) == 0)
-	u.ieee.mantissa1 |= fetestexcept (FE_INEXACT) != 0;
+	u.ieee.mantissa1 |= libc_fetestexcept (FE_INEXACT) != 0;
       v.d = a1 + u.d;
-      int j = fetestexcept (FE_INEXACT) != 0;
-      feupdateenv (&env);
+      int j = libc_fetestexcept (FE_INEXACT) != 0;
+      libc_feupdateenv (&env);
       /* Ensure the following computations are performed in default rounding
 	 mode instead of just reusing the round to zero computation.  */
       asm volatile ("" : "=m" (u) : "m" (u));
diff --git a/libc/sysdeps/ieee754/dbl-64/s_fmaf.c b/libc/sysdeps/ieee754/dbl-64/s_fmaf.c
index cd16cd1dc..dc748e554 100644
--- a/libc/sysdeps/ieee754/dbl-64/s_fmaf.c
+++ b/libc/sysdeps/ieee754/dbl-64/s_fmaf.c
@@ -1,5 +1,5 @@
 /* Compute x * y + z as ternary operation.
-   Copyright (C) 2010 Free Software Foundation, Inc.
+   Copyright (C) 2010, 2011 Free Software Foundation, Inc.
    This file is part of the GNU C Library.
    Contributed by Jakub Jelinek <jakub@redhat.com>, 2010.
 
@@ -21,6 +21,7 @@
 #include <math.h>
 #include <fenv.h>
 #include <ieee754.h>
+#include <math_private.h>
 
 /* This implementation relies on double being more than twice as
    precise as float and uses rounding to odd in order to avoid problems
@@ -35,13 +36,12 @@ __fmaf (float x, float y, float z)
   /* Multiplication is always exact.  */
   double temp = (double) x * (double) y;
   union ieee754_double u;
-  feholdexcept (&env);
-  fesetround (FE_TOWARDZERO);
+  libc_feholdexcept_setroundf (&env, FE_TOWARDZERO);
   /* Perform addition with round to odd.  */
   u.d = temp + (double) z;
   if ((u.ieee.mantissa1 & 1) == 0 && u.ieee.exponent != 0x7ff)
-    u.ieee.mantissa1 |= fetestexcept (FE_INEXACT) != 0;
-  feupdateenv (&env);
+    u.ieee.mantissa1 |= libc_fetestexcept (FE_INEXACT) != 0;
+  libc_feupdateenv (&env);
   /* And finally truncation with round to nearest.  */
   return (float) u.d;
 }
diff --git a/libc/sysdeps/ieee754/dbl-64/s_isinf_ns.c b/libc/sysdeps/ieee754/dbl-64/s_isinf_ns.c
new file mode 100644
index 000000000..065522ed8
--- /dev/null
+++ b/libc/sysdeps/ieee754/dbl-64/s_isinf_ns.c
@@ -0,0 +1,20 @@
+/*
+ * Written by Ulrich Drepper <drepper@gmail.com>.
+ */
+
+/*
+ * __isinf_ns(x) returns != 0 if x is ±inf, else 0;
+ * no branching!
+ */
+
+#include "math.h"
+#include "math_private.h"
+
+#undef __isinf_ns
+int
+__isinf_ns (double x)
+{
+	int32_t hx,lx;
+	EXTRACT_WORDS(hx,lx,x);
+	return !(lx | ((hx & 0x7fffffff) ^ 0x7ff00000));
+}
diff --git a/libc/sysdeps/ieee754/dbl-64/s_isnan.c b/libc/sysdeps/ieee754/dbl-64/s_isnan.c
index 3a089e99a..74e829160 100644
--- a/libc/sysdeps/ieee754/dbl-64/s_isnan.c
+++ b/libc/sysdeps/ieee754/dbl-64/s_isnan.c
@@ -22,6 +22,7 @@ static char rcsid[] = "$NetBSD: s_isnan.c,v 1.8 1995/05/10 20:47:36 jtc Exp $";
 #include "math.h"
 #include "math_private.h"
 
+#undef __isnan
 #ifdef __STDC__
 	int __isnan(double x)
 #else
diff --git a/libc/sysdeps/ieee754/dbl-64/s_lround.c b/libc/sysdeps/ieee754/dbl-64/s_lround.c
index 4e1302ad4..a849997a1 100644
--- a/libc/sysdeps/ieee754/dbl-64/s_lround.c
+++ b/libc/sysdeps/ieee754/dbl-64/s_lround.c
@@ -51,7 +51,7 @@ __lround (double x)
   else if (j0 < (int32_t) (8 * sizeof (long int)) - 1)
     {
       if (j0 >= 52)
-	result = ((long int) i0 << (j0 - 20)) | (i1 << (j0 - 52));
+	result = ((long int) i0 << (j0 - 20)) | ((long int) i1 << (j0 - 52));
       else
 	{
 	  u_int32_t j = i1 + (0x80000000 >> (j0 - 20));
diff --git a/libc/sysdeps/ieee754/dbl-64/s_rint.c b/libc/sysdeps/ieee754/dbl-64/s_rint.c
index 4e6381efb..a671a6277 100644
--- a/libc/sysdeps/ieee754/dbl-64/s_rint.c
+++ b/libc/sysdeps/ieee754/dbl-64/s_rint.c
@@ -10,10 +10,6 @@
  * ====================================================
  */
 
-#if defined(LIBM_SCCS) && !defined(lint)
-static char rcsid[] = "$NetBSD: s_rint.c,v 1.8 1995/05/10 20:48:04 jtc Exp $";
-#endif
-
 /*
  * rint(x)
  * Return x rounded to integral value according to the prevailing
@@ -27,22 +23,14 @@ static char rcsid[] = "$NetBSD: s_rint.c,v 1.8 1995/05/10 20:48:04 jtc Exp $";
 #include "math.h"
 #include "math_private.h"
 
-#ifdef __STDC__
 static const double
-#else
-static double
-#endif
 TWO52[2]={
   4.50359962737049600000e+15, /* 0x43300000, 0x00000000 */
  -4.50359962737049600000e+15, /* 0xC3300000, 0x00000000 */
 };
 
-#ifdef __STDC__
-	double __rint(double x)
-#else
-	double __rint(x)
-	double x;
-#endif
+double
+__rint(double x)
 {
 	int32_t i0,j0,sx;
 	u_int32_t i,i1;
@@ -57,11 +45,11 @@ TWO52[2]={
 		i0 &= 0xfffe0000;
 		i0 |= ((i1|-i1)>>12)&0x80000;
 		SET_HIGH_WORD(x,i0);
-	        w = TWO52[sx]+x;
-	        t =  w-TWO52[sx];
+		w = TWO52[sx]+x;
+		t =  w-TWO52[sx];
 		GET_HIGH_WORD(i0,t);
 		SET_HIGH_WORD(t,(i0&0x7fffffff)|(sx<<31));
-	        return t;
+		return t;
 	    } else {
 		i = (0x000fffff)>>j0;
 		if(((i0&i)|i1)==0) return x; /* x is integral */
@@ -91,8 +79,10 @@ TWO52[2]={
 	w = TWO52[sx]+x;
 	return w-TWO52[sx];
 }
+#ifndef __rint
 weak_alias (__rint, rint)
-#ifdef NO_LONG_DOUBLE
+# ifdef NO_LONG_DOUBLE
 strong_alias (__rint, __rintl)
 weak_alias (__rint, rintl)
+# endif
 #endif
diff --git a/libc/sysdeps/ieee754/dbl-64/s_tan.c b/libc/sysdeps/ieee754/dbl-64/s_tan.c
index 4e26d90ae..df8eedd92 100644
--- a/libc/sysdeps/ieee754/dbl-64/s_tan.c
+++ b/libc/sysdeps/ieee754/dbl-64/s_tan.c
@@ -1,7 +1,7 @@
 /*
  * IBM Accurate Mathematical Library
  * written by International Business Machines Corp.
- * Copyright (C) 2001, 2009 Free Software Foundation
+ * Copyright (C) 2001, 2009, 2011 Free Software Foundation
  *
  * This program is free software; you can redistribute it and/or modify
  * it under the terms of the GNU Lesser General Public License as published by
@@ -36,7 +36,7 @@
 
 #include <errno.h>
 #include "endian.h"
-#include "dla.h"
+#include <dla.h>
 #include "mpa.h"
 #include "MathLib.h"
 #include "math.h"
@@ -50,7 +50,10 @@ double tan(double x) {
 
   int ux,i,n;
   double a,da,a2,b,db,c,dc,c1,cc1,c2,cc2,c3,cc3,fi,ffi,gi,pz,s,sy,
-  t,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,w,x2,xn,xx2,y,ya,yya,z0,z,zz,z2,zz2;
+  t,t1,t2,t3,t4,t7,t8,t9,t10,w,x2,xn,xx2,y,ya,yya,z0,z,zz,z2,zz2;
+#ifndef DLA_FMA
+  double t5,t6;
+#endif
   int p;
   number num,v;
   mp_no mpa,mpt1,mpt2;
@@ -84,7 +87,7 @@ double tan(double x) {
 
     /* Second stage */
     c1 = x2*(a15.d+x2*(a17.d+x2*(a19.d+x2*(a21.d+x2*(a23.d+x2*(a25.d+
-         x2*a27.d))))));
+	 x2*a27.d))))));
     EMULV(x,x,x2,xx2,t1,t2,t3,t4,t5)
     ADD2(a13.d,aa13.d,c1,zero.d,c2,cc2,t1,t2)
     MUL2(x2,xx2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
@@ -159,13 +162,13 @@ double tan(double x) {
       a2 = a*a;
       t2 = da+a*a2*(d3.d+a2*(d5.d+a2*(d7.d+a2*(d9.d+a2*d11.d))));
       if (n) {
-        /* First stage -cot */
-        EADD(a,t2,b,db)
-        DIV2(one.d,zero.d,b,db,c,dc,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
-        if ((y=c+(dc-u6.d*c))==c+(dc+u6.d*c))  return (-y); }
+	/* First stage -cot */
+	EADD(a,t2,b,db)
+	DIV2(one.d,zero.d,b,db,c,dc,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
+	if ((y=c+(dc-u6.d*c))==c+(dc+u6.d*c))  return (-y); }
       else {
-        /* First stage tan */
-        if ((y=a+(t2-u5.d*a))==a+(t2+u5.d*a))  return y; }
+	/* First stage tan */
+	if ((y=a+(t2-u5.d*a))==a+(t2+u5.d*a))  return y; }
       /* Second stage */
       /* Range reduction by algorithm ii */
       t = (x*hpinv.d + toint.d);
@@ -184,7 +187,7 @@ double tan(double x) {
       EADD(a,da,t1,t2)   a=t1;  da=t2;
       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+
-           x2*a27.d))))));
+	   x2*a27.d))))));
       ADD2(a13.d,aa13.d,c1,zero.d,c2,cc2,t1,t2)
       MUL2(x2,xx2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
       ADD2(a11.d,aa11.d,c1,cc1,c2,cc2,t1,t2)
@@ -201,12 +204,12 @@ double tan(double x) {
       ADD2(a  ,da  ,c2,cc2,c1,cc1,t1,t2)
 
       if (n) {
-        /* Second stage -cot */
-        DIV2(one.d,zero.d,c1,cc1,c2,cc2,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
-        if ((y=c2+(cc2-u8.d*c2)) == c2+(cc2+u8.d*c2))  return (-y); }
+	/* Second stage -cot */
+	DIV2(one.d,zero.d,c1,cc1,c2,cc2,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
+	if ((y=c2+(cc2-u8.d*c2)) == c2+(cc2+u8.d*c2))  return (-y); }
       else {
-        /* Second stage tan */
-        if ((y=c1+(cc1-u7.d*c1)) == c1+(cc1+u7.d*c1))  return y; }
+	/* Second stage tan */
+	if ((y=c1+(cc1-u7.d*c1)) == c1+(cc1+u7.d*c1))  return y; }
       return tanMp(x);
     }
 
@@ -287,18 +290,18 @@ double tan(double x) {
       a2 = a*a;
       t2 = da+a*a2*(d3.d+a2*(d5.d+a2*(d7.d+a2*(d9.d+a2*d11.d))));
       if (n) {
-        /* First stage -cot */
-        EADD(a,t2,b,db)
-        DIV2(one.d,zero.d,b,db,c,dc,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
-        if ((y=c+(dc-u14.d*c))==c+(dc+u14.d*c))  return (-y); }
+	/* First stage -cot */
+	EADD(a,t2,b,db)
+	DIV2(one.d,zero.d,b,db,c,dc,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
+	if ((y=c+(dc-u14.d*c))==c+(dc+u14.d*c))  return (-y); }
       else {
-        /* First stage tan */
-        if ((y=a+(t2-u13.d*a))==a+(t2+u13.d*a))  return y; }
+	/* First stage tan */
+	if ((y=a+(t2-u13.d*a))==a+(t2+u13.d*a))  return y; }
 
       /* Second stage */
       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+
-           x2*a27.d))))));
+	   x2*a27.d))))));
       ADD2(a13.d,aa13.d,c1,zero.d,c2,cc2,t1,t2)
       MUL2(x2,xx2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
       ADD2(a11.d,aa11.d,c1,cc1,c2,cc2,t1,t2)
@@ -315,12 +318,12 @@ double tan(double x) {
       ADD2(a  ,da  ,c2,cc2,c1,cc1,t1,t2)
 
       if (n) {
-        /* Second stage -cot */
-        DIV2(one.d,zero.d,c1,cc1,c2,cc2,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
-        if ((y=c2+(cc2-u16.d*c2)) == c2+(cc2+u16.d*c2))  return (-y); }
+	/* Second stage -cot */
+	DIV2(one.d,zero.d,c1,cc1,c2,cc2,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
+	if ((y=c2+(cc2-u16.d*c2)) == c2+(cc2+u16.d*c2))  return (-y); }
       else {
-        /* Second stage tan */
-        if ((y=c1+(cc1-u15.d*c1)) == c1+(cc1+u15.d*c1))  return (y); }
+	/* Second stage tan */
+	if ((y=c1+(cc1-u15.d*c1)) == c1+(cc1+u15.d*c1))  return (y); }
       return tanMp(x);
     }
 
@@ -404,7 +407,7 @@ double tan(double x) {
 
     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+
-         x2*a27.d))))));
+	 x2*a27.d))))));
     ADD2(a13.d,aa13.d,c1,zero.d,c2,cc2,t1,t2)
     MUL2(x2,xx2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
     ADD2(a11.d,aa11.d,c1,cc1,c2,cc2,t1,t2)
diff --git a/libc/sysdeps/ieee754/dbl-64/w_exp.c b/libc/sysdeps/ieee754/dbl-64/w_exp.c
index 121649209..f1becff94 100644
--- a/libc/sysdeps/ieee754/dbl-64/w_exp.c
+++ b/libc/sysdeps/ieee754/dbl-64/w_exp.c
@@ -1,55 +1,46 @@
-/* @(#)w_exp.c 5.1 93/09/24 */
-/*
- * ====================================================
- * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
- *
- * Developed at SunPro, a Sun Microsystems, Inc. business.
- * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
- * is preserved.
- * ====================================================
- */
-
-#if defined(LIBM_SCCS) && !defined(lint)
-static char rcsid[] = "$NetBSD: w_exp.c,v 1.6 1995/05/10 20:48:51 jtc Exp $";
-#endif
+/* Copyright (C) 2011 Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
+   Contributed by Ulrich Drepper <drepper@gmail.com>, 2011.
+
+   The GNU C Library 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.1 of the License, or (at your option) any later version.
 
-/*
- * wrapper exp(x)
- */
+   The GNU C Library 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 GNU
+   Lesser General Public License for more details.
 
-#include "math.h"
-#include "math_private.h"
+   You should have received a copy of the GNU Lesser General Public
+   License along with the GNU C Library; if not, write to the Free
+   Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
+   02111-1307 USA.  */
+
+#include <math.h>
+#include <math_private.h>
 
-#ifdef __STDC__
 static const double
-#else
-static double
-#endif
 o_threshold=  7.09782712893383973096e+02,  /* 0x40862E42, 0xFEFA39EF */
 u_threshold= -7.45133219101941108420e+02;  /* 0xc0874910, 0xD52D3051 */
 
-#ifdef __STDC__
-	double __exp(double x)		/* wrapper exp */
-#else
-	double __exp(x)			/* wrapper exp */
-	double x;
-#endif
+
+/* wrapper exp */
+double
+__exp (double x)
 {
-#ifdef _IEEE_LIBM
-	return __ieee754_exp(x);
-#else
-	double z;
-	z = __ieee754_exp(x);
-	if(_LIB_VERSION == _IEEE_) return z;
-	if(__finite(x)) {
-	    if(x>o_threshold)
-	        return __kernel_standard(x,x,6); /* exp overflow */
-	    else if(x<u_threshold)
-	        return __kernel_standard(x,x,7); /* exp underflow */
-	}
-	return z;
-#endif
+  if (__builtin_expect (x > o_threshold, 0))
+    {
+      if (_LIB_VERSION != _IEEE_)
+	return __kernel_standard_f (x, x, 6);
+    }
+  else if (__builtin_expect (x < u_threshold, 0))
+    {
+      if (_LIB_VERSION != _IEEE_)
+	return __kernel_standard_f (x, x, 7);
+    }
+
+  return __ieee754_exp (x);
 }
 hidden_def (__exp)
 weak_alias (__exp, exp)
diff --git a/libc/sysdeps/ieee754/dbl-64/wordsize-64/e_cosh.c b/libc/sysdeps/ieee754/dbl-64/wordsize-64/e_cosh.c
new file mode 100644
index 000000000..41dc42c0a
--- /dev/null
+++ b/libc/sysdeps/ieee754/dbl-64/wordsize-64/e_cosh.c
@@ -0,0 +1,82 @@
+/* Optimized by Ulrich Drepper <drepper@gmail.com>, 2011 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+/* __ieee754_cosh(x)
+ * Method :
+ * mathematically cosh(x) if defined to be (exp(x)+exp(-x))/2
+ *	1. Replace x by |x| (cosh(x) = cosh(-x)).
+ *	2.
+ *							[ exp(x) - 1 ]^2
+ *	    0        <= x <= ln2/2  :  cosh(x) := 1 + -------------------
+ *							   2*exp(x)
+ *
+ *						  exp(x) +  1/exp(x)
+ *	    ln2/2    <= x <= 22     :  cosh(x) := -------------------
+ *							  2
+ *	    22       <= x <= lnovft :  cosh(x) := exp(x)/2
+ *	    lnovft   <= x <= ln2ovft:  cosh(x) := exp(x/2)/2 * exp(x/2)
+ *	    ln2ovft  <  x	    :  cosh(x) := huge*huge (overflow)
+ *
+ * Special cases:
+ *	cosh(x) is |x| if x is +INF, -INF, or NaN.
+ *	only cosh(0)=1 is exact for finite x.
+ */
+
+#include "math.h"
+#include "math_private.h"
+
+static const double one = 1.0, half=0.5, huge = 1.0e300;
+
+double
+__ieee754_cosh (double x)
+{
+	double t,w;
+	int32_t ix;
+
+    /* High word of |x|. */
+	GET_HIGH_WORD(ix,x);
+	ix &= 0x7fffffff;
+
+    /* |x| in [0,22] */
+	if (ix < 0x40360000) {
+	    /* |x| in [0,0.5*ln2], return 1+expm1(|x|)^2/(2*exp(|x|)) */
+		if(ix<0x3fd62e43) {
+		    t = __expm1(fabs(x));
+		    w = one+t;
+		    if (ix<0x3c800000) return w;	/* cosh(tiny) = 1 */
+		    return one+(t*t)/(w+w);
+		}
+
+	    /* |x| in [0.5*ln2,22], return (exp(|x|)+1/exp(|x|)/2; */
+		t = __ieee754_exp(fabs(x));
+		return half*t+half/t;
+	}
+
+    /* |x| in [22, log(maxdouble)] return half*exp(|x|) */
+	if (ix < 0x40862e42)  return half*__ieee754_exp(fabs(x));
+
+    /* |x| in [log(maxdouble), overflowthresold] */
+	int64_t fix;
+	EXTRACT_WORDS64(fix, x);
+	if (fix <= UINT64_C(0x408633ce8fb9f87d)) {
+	    w = __ieee754_exp(half*fabs(x));
+	    t = half*w;
+	    return t*w;
+	}
+
+    /* x is INF or NaN */
+	if(ix>=0x7ff00000) return x*x;
+
+    /* |x| > overflowthresold, cosh(x) overflow */
+	return huge*huge;
+}
+strong_alias (__ieee754_cosh, __cosh_finite)
diff --git a/libc/sysdeps/ieee754/dbl-64/wordsize-64/s_ceil.c b/libc/sysdeps/ieee754/dbl-64/wordsize-64/s_ceil.c
index 9123fdc7b..e0e71558f 100644
--- a/libc/sysdeps/ieee754/dbl-64/wordsize-64/s_ceil.c
+++ b/libc/sysdeps/ieee754/dbl-64/wordsize-64/s_ceil.c
@@ -22,18 +22,10 @@
 #include "math.h"
 #include "math_private.h"
 
-#ifdef __STDC__
 static const double huge = 1.0e300;
-#else
-static double huge = 1.0e300;
-#endif
 
-#ifdef __STDC__
-	double __ceil(double x)
-#else
-	double __ceil(x)
-	double x;
-#endif
+double
+__ceil(double x)
 {
 	int64_t i0,i;
 	int32_t j0;
@@ -60,8 +52,10 @@ static double huge = 1.0e300;
 	INSERT_WORDS64(x,i0);
 	return x;
 }
+#ifndef __ceil
 weak_alias (__ceil, ceil)
-#ifdef NO_LONG_DOUBLE
+# ifdef NO_LONG_DOUBLE
 strong_alias (__ceil, __ceill)
 weak_alias (__ceil, ceill)
+# endif
 #endif
diff --git a/libc/sysdeps/ieee754/dbl-64/wordsize-64/s_finite.c b/libc/sysdeps/ieee754/dbl-64/wordsize-64/s_finite.c
new file mode 100644
index 000000000..93a39a683
--- /dev/null
+++ b/libc/sysdeps/ieee754/dbl-64/wordsize-64/s_finite.c
@@ -0,0 +1,33 @@
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+/*
+ * finite(x) returns 1 is x is finite, else 0;
+ * no branching!
+ */
+
+#include "math.h"
+#include "math_private.h"
+
+#undef __finite
+int
+__finite(double x)
+{
+  int64_t lx;
+  EXTRACT_WORDS64(lx,x);
+  return (int)((uint64_t)((lx&INT64_C(0x7fffffffffffffff))-INT64_C(0x7ff0000000000000))>>63);
+}
+hidden_def (__finite)
+weak_alias (__finite, finite)
+#ifdef NO_LONG_DOUBLE
+strong_alias (__finite, __finitel)
+weak_alias (__finite, finitel)
+#endif
diff --git a/libc/sysdeps/ieee754/dbl-64/wordsize-64/s_floor.c b/libc/sysdeps/ieee754/dbl-64/wordsize-64/s_floor.c
new file mode 100644
index 000000000..8b7300bb9
--- /dev/null
+++ b/libc/sysdeps/ieee754/dbl-64/wordsize-64/s_floor.c
@@ -0,0 +1,81 @@
+/* Round double to integer away from zero.
+   Copyright (C) 2011 Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
+   Contributed by Ulrich Drepper <drepper@cygnus.com>, 2011.
+
+   The GNU C Library 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.1 of the License, or (at your option) any later version.
+
+   The GNU C Library 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 GNU
+   Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public
+   License along with the GNU C Library; if not, write to the Free
+   Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
+   02111-1307 USA.  */
+
+/* Based on a version which carries the following copyright:  */
+
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+#include "math.h"
+#include "math_private.h"
+
+/*
+ * floor(x)
+ * Return x rounded toward -inf to integral value
+ * Method:
+ *	Bit twiddling.
+ * Exception:
+ *	Inexact flag raised if x not equal to floor(x).
+ */
+
+static const double huge = 1.0e300;
+
+
+double
+__floor (double x)
+{
+	int64_t i0;
+	EXTRACT_WORDS64(i0,x);
+	int32_t j0 = ((i0>>52)&0x7ff)-0x3ff;
+	if(__builtin_expect(j0<52, 1)) {
+	    if(j0<0) {	/* raise inexact if x != 0 */
+		if(huge+x>0.0) {/* return 0*sign(x) if |x|<1 */
+		    if(i0>=0) {i0=0;}
+		    else if((i0&0x7fffffffffffffffl)!=0)
+			{ i0=0xbff0000000000000l;}
+		}
+	    } else {
+		uint64_t i = (0x000fffffffffffffl)>>j0;
+		if((i0&i)==0) return x; /* x is integral */
+		if(huge+x>0.0) {	/* raise inexact flag */
+		    if(i0<0) i0 += (0x0010000000000000l)>>j0;
+		    i0 &= (~i);
+		}
+	    }
+	    INSERT_WORDS64(x,i0);
+	} else if (j0==0x400)
+	    return x+x;	/* inf or NaN */
+	return x;
+}
+#ifndef __floor
+weak_alias (__floor, floor)
+# ifdef NO_LONG_DOUBLE
+strong_alias (__floor, __floorl)
+weak_alias (__floor, floorl)
+# endif
+#endif
diff --git a/libc/sysdeps/ieee754/dbl-64/wordsize-64/s_frexp.c b/libc/sysdeps/ieee754/dbl-64/wordsize-64/s_frexp.c
new file mode 100644
index 000000000..76f761080
--- /dev/null
+++ b/libc/sysdeps/ieee754/dbl-64/wordsize-64/s_frexp.c
@@ -0,0 +1,67 @@
+/* Copyright (C) 2011 Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
+   Contributed by Ulrich Drepper <drepper@gmail.com>, 2011.
+
+   The GNU C Library 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.1 of the License, or (at your option) any later version.
+
+   The GNU C Library 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 GNU
+   Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public
+   License along with the GNU C Library; if not, write to the Free
+   Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
+   02111-1307 USA.  */
+
+#include <inttypes.h>
+#include "math.h"
+#include "math_private.h"
+
+/*
+ * for non-zero, finite x
+ *	x = frexp(arg,&exp);
+ * return a double fp quantity x such that 0.5 <= |x| <1.0
+ * and the corresponding binary exponent "exp". That is
+ *	arg = x*2^exp.
+ * If arg is inf, 0.0, or NaN, then frexp(arg,&exp) returns arg
+ * with *exp=0.
+ */
+
+
+double
+__frexp (double x, int *eptr)
+{
+  int64_t ix;
+  EXTRACT_WORDS64 (ix, x);
+  int32_t ex = 0x7ff & (ix >> 52);
+  int e = 0;
+
+  if (__builtin_expect (ex != 0x7ff && x != 0.0, 1))
+    {
+      /* Not zero and finite.  */
+      e = ex - 1022;
+      if (__builtin_expect (ex == 0, 0))
+	{
+	  /* Subnormal.  */
+	  x *= 0x1p54;
+	  EXTRACT_WORDS64 (ix, x);
+	  ex = 0x7ff & (ix >> 52);
+	  e = ex - 1022 - 54;
+	}
+
+      ix = (ix & INT64_C (0x800fffffffffffff)) | INT64_C (0x3fe0000000000000);
+      INSERT_WORDS64 (x, ix);
+    }
+
+  *eptr = e;
+  return x;
+}
+weak_alias (__frexp, frexp)
+#ifdef NO_LONG_DOUBLE
+strong_alias (__frexp, __frexpl)
+weak_alias (__frexp, frexpl)
+#endif
diff --git a/libc/sysdeps/ieee754/dbl-64/wordsize-64/s_isinf_ns.c b/libc/sysdeps/ieee754/dbl-64/wordsize-64/s_isinf_ns.c
new file mode 100644
index 000000000..09dcc9453
--- /dev/null
+++ b/libc/sysdeps/ieee754/dbl-64/wordsize-64/s_isinf_ns.c
@@ -0,0 +1,20 @@
+/*
+ * Written by Ulrich Drepper <drepper@gmail.com>.
+ */
+
+/*
+ * __isinf_ns(x) returns != 0 if x is ±inf, else 0;
+ * no branching!
+ */
+
+#include "math.h"
+#include "math_private.h"
+
+#undef __isinf_ns
+int
+__isinf_ns (double x)
+{
+	int64_t ix;
+	EXTRACT_WORDS64(ix,x);
+	return (ix & UINT64_C(0x7fffffffffffffff)) == UINT64_C(0x7ff0000000000000);
+}
diff --git a/libc/sysdeps/ieee754/dbl-64/wordsize-64/s_isnan.c b/libc/sysdeps/ieee754/dbl-64/wordsize-64/s_isnan.c
index 3b08c54dd..65dc8934f 100644
--- a/libc/sysdeps/ieee754/dbl-64/wordsize-64/s_isnan.c
+++ b/libc/sysdeps/ieee754/dbl-64/wordsize-64/s_isnan.c
@@ -18,6 +18,7 @@
 #include "math.h"
 #include "math_private.h"
 
+#undef __isnan
 #ifdef __STDC__
 	int __isnan(double x)
 #else
diff --git a/libc/sysdeps/ieee754/dbl-64/wordsize-64/s_logb.c b/libc/sysdeps/ieee754/dbl-64/wordsize-64/s_logb.c
new file mode 100644
index 000000000..741350436
--- /dev/null
+++ b/libc/sysdeps/ieee754/dbl-64/wordsize-64/s_logb.c
@@ -0,0 +1,44 @@
+/* Compute radix independent exponent.
+   Copyright (C) 2011 Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
+   Contributed by Ulrich Drepper <drepper@gmail.com>, 2011.
+
+   The GNU C Library 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.1 of the License, or (at your option) any later version.
+
+   The GNU C Library 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 GNU
+   Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public
+   License along with the GNU C Library; if not, write to the Free
+   Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
+   02111-1307 USA.  */
+
+#include <math.h>
+
+#include "math_private.h"
+
+
+double
+__logb (double x)
+{
+  int64_t ix;
+
+  EXTRACT_WORDS64 (ix, x);
+  ix &= UINT64_C(0x7fffffffffffffff);
+  if (ix == 0)
+    return -1.0 / fabs (x);
+  unsigned int ex = ix >> 52;
+  if (ex == 0x7ff)
+    return x * x;
+  return ex == 0 ? -1022.0 : (double) (ex - 1023);
+}
+weak_alias (__logb, logb)
+#ifdef NO_LONG_DOUBLE
+strong_alias (__logb, __logbl)
+weak_alias (__logb, logbl)
+#endif
diff --git a/libc/sysdeps/ieee754/dbl-64/wordsize-64/s_nearbyint.c b/libc/sysdeps/ieee754/dbl-64/wordsize-64/s_nearbyint.c
index cb49019dd..861da20b1 100644
--- a/libc/sysdeps/ieee754/dbl-64/wordsize-64/s_nearbyint.c
+++ b/libc/sysdeps/ieee754/dbl-64/wordsize-64/s_nearbyint.c
@@ -24,22 +24,14 @@
 #include "math.h"
 #include "math_private.h"
 
-#ifdef __STDC__
 static const double
-#else
-static double
-#endif
 TWO52[2]={
   4.50359962737049600000e+15, /* 0x43300000, 0x00000000 */
  -4.50359962737049600000e+15, /* 0xC3300000, 0x00000000 */
 };
 
-#ifdef __STDC__
-	double __nearbyint(double x)
-#else
-	double __nearbyint(x)
-	double x;
-#endif
+double
+__nearbyint(double x)
 {
 	fenv_t env;
 	int64_t i0,sx;
@@ -47,20 +39,19 @@ TWO52[2]={
 	EXTRACT_WORDS64(i0,x);
 	sx = (i0>>63)&1;
 	j0 = ((i0>>52)&0x7ff)-0x3ff;
-	if(j0<52) {
+	if(__builtin_expect(j0<52, 1)) {
 	    if(j0<0) {
 	      if((i0&UINT64_C(0x7fffffffffffffff))==0) return x;
 		uint64_t i = i0 & UINT64_C(0xfffffffffffff);
 		i0 &= UINT64_C(0xfffe000000000000);
 		i0 |= (((i|-i) >> 12) & UINT64_C(0x8000000000000));
 		INSERT_WORDS64(x,i0);
-		feholdexcept (&env);
+		libc_feholdexcept (&env);
 		double w = TWO52[sx]+x;
 		double t =  w-TWO52[sx];
-		fesetenv (&env);
-		EXTRACT_WORDS64(i0,t);
-		INSERT_WORDS64(t,(i0&UINT64_C(0x7fffffffffffffff))|(sx<<63));
-		return t;
+		math_opt_barrier(t);
+		libc_fesetenv (&env);
+		return copysign(t, x);
 	    } else {
 		uint64_t i = UINT64_C(0x000fffffffffffff)>>j0;
 		if((i0&i)==0) return x; /* x is integral */
@@ -73,10 +64,11 @@ TWO52[2]={
 	    else return x;		/* x is integral */
 	}
 	INSERT_WORDS64(x,i0);
-	feholdexcept (&env);
+	libc_feholdexcept (&env);
 	double w = TWO52[sx]+x;
 	double t = w-TWO52[sx];
-	fesetenv (&env);
+	math_opt_barrier (t);
+	libc_fesetenv (&env);
 	return t;
 }
 weak_alias (__nearbyint, nearbyint)
diff --git a/libc/sysdeps/ieee754/dbl-64/wordsize-64/s_rint.c b/libc/sysdeps/ieee754/dbl-64/wordsize-64/s_rint.c
index 4a60aa327..571b3811a 100644
--- a/libc/sysdeps/ieee754/dbl-64/wordsize-64/s_rint.c
+++ b/libc/sysdeps/ieee754/dbl-64/wordsize-64/s_rint.c
@@ -1,4 +1,3 @@
-/* @(#)s_rint.c 5.1 93/09/24 */
 /*
  * ====================================================
  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
@@ -23,22 +22,14 @@
 #include "math.h"
 #include "math_private.h"
 
-#ifdef __STDC__
 static const double
-#else
-static double
-#endif
 TWO52[2]={
   4.50359962737049600000e+15, /* 0x43300000, 0x00000000 */
  -4.50359962737049600000e+15, /* 0xC3300000, 0x00000000 */
 };
 
-#ifdef __STDC__
-	double __rint(double x)
-#else
-	double __rint(x)
-	double x;
-#endif
+double
+__rint(double x)
 {
 	int64_t i0,sx;
 	int32_t j0;
@@ -72,8 +63,10 @@ TWO52[2]={
 	double w = TWO52[sx]+x;
 	return w-TWO52[sx];
 }
+#ifndef __rint
 weak_alias (__rint, rint)
-#ifdef NO_LONG_DOUBLE
+# ifdef NO_LONG_DOUBLE
 strong_alias (__rint, __rintl)
 weak_alias (__rint, rintl)
+# endif
 #endif