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

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

21 files changed, 148 insertions, 155 deletions

diff --git a/libc/sysdeps/ieee754/ldbl-128ibm/e_acoshl.c b/libc/sysdeps/ieee754/ldbl-128ibm/e_acoshl.c
index 00576c76c..20d94eaa1 100644
--- a/libc/sysdeps/ieee754/ldbl-128ibm/e_acoshl.c
+++ b/libc/sysdeps/ieee754/ldbl-128ibm/e_acoshl.c
@@ -10,10 +10,6 @@
  * ====================================================
  */
 
-#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
@@ -31,20 +27,12 @@ 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 long double
-#else
-static long double
-#endif
 one	= 1.0L,
 ln2	= 6.93147180559945286227e-01L;  /* 0x3FE62E42, 0xFEFA39EF */
 
-#ifdef __STDC__
-	long double __ieee754_acoshl(long double x)
-#else
-	long double __ieee754_acoshl(x)
-	long double x;
-#endif
+long double
+__ieee754_acoshl(long double x)
 {
 	long double t;
 	int64_t hx;
@@ -54,7 +42,7 @@ ln2	= 6.93147180559945286227e-01L;  /* 0x3FE62E42, 0xFEFA39EF */
 	    return (x-x)/(x-x);
 	} else if(hx >=0x41b0000000000000LL) {	/* x > 2**28 */
 	    if(hx >=0x7ff0000000000000LL) {	/* x is inf of NaN */
-	        return x+x;
+		return x+x;
 	    } else
 		return __ieee754_logl(x)+ln2;	/* acosh(huge)=log(2x) */
 	} else if (((hx-0x3ff0000000000000LL)|(lx&0x7fffffffffffffffLL))==0) {
@@ -67,3 +55,4 @@ ln2	= 6.93147180559945286227e-01L;  /* 0x3FE62E42, 0xFEFA39EF */
 	    return __log1p(t+__sqrtl(2.0*t+t*t));
 	}
 }
+strong_alias (__ieee754_acoshl, __acoshl_finite)
diff --git a/libc/sysdeps/ieee754/ldbl-128ibm/e_acosl.c b/libc/sysdeps/ieee754/ldbl-128ibm/e_acosl.c
index 8823fd69b..1b37c9220 100644
--- a/libc/sysdeps/ieee754/ldbl-128ibm/e_acosl.c
+++ b/libc/sysdeps/ieee754/ldbl-128ibm/e_acosl.c
@@ -283,11 +283,11 @@ __ieee754_acosl (x)
       s = __ieee754_sqrtl (z);
       /* Compute an extended precision square root from
 	 the Newton iteration  s -> 0.5 * (s + z / s).
-         The change w from s to the improved value is
+	 The change w from s to the improved value is
 	    w = 0.5 * (s + z / s) - s  = (s^2 + z)/2s - s = (z - s^2)/2s.
-          Express s = f1 + f2 where f1 * f1 is exactly representable.
+	  Express s = f1 + f2 where f1 * f1 is exactly representable.
 	  w = (z - s^2)/2s = (z - f1^2 - 2 f1 f2 - f2^2)/2s .
-          s + w has extended precision.  */
+	  s + w has extended precision.  */
       u.value = s;
       u.parts32.w2 = 0;
       u.parts32.w3 = 0;
@@ -326,3 +326,4 @@ __ieee754_acosl (x)
       return 2.0 * w;
     }
 }
+strong_alias (__ieee754_acosl, __acosl_finite)
diff --git a/libc/sysdeps/ieee754/ldbl-128ibm/e_asinl.c b/libc/sysdeps/ieee754/ldbl-128ibm/e_asinl.c
index 3696694f7..6c61232c0 100644
--- a/libc/sysdeps/ieee754/ldbl-128ibm/e_asinl.c
+++ b/libc/sysdeps/ieee754/ldbl-128ibm/e_asinl.c
@@ -12,9 +12,9 @@
 /*
   Long double expansions are
   Copyright (C) 2001 Stephen L. Moshier <moshier@na-net.ornl.gov>
-  and are incorporated herein by permission of the author.  The author 
-  reserves the right to distribute this material elsewhere under different 
-  copying permissions.  These modifications are distributed here under the 
+  and are incorporated herein by permission of the author.  The author
+  reserves the right to distribute this material elsewhere under different
+  copying permissions.  These modifications are distributed here under the
   following terms:
 
     This library is free software; you can redistribute it and/or
@@ -263,3 +263,4 @@ __ieee754_asinl (x)
   else
     return -t;
 }
+strong_alias (__ieee754_asinl, __asinl_finite)
diff --git a/libc/sysdeps/ieee754/ldbl-128ibm/e_atan2l.c b/libc/sysdeps/ieee754/ldbl-128ibm/e_atan2l.c
index a4bb53df0..a5b662100 100644
--- a/libc/sysdeps/ieee754/ldbl-128ibm/e_atan2l.c
+++ b/libc/sysdeps/ieee754/ldbl-128ibm/e_atan2l.c
@@ -17,7 +17,7 @@
  * Method :
  *	1. Reduce y to positive by atan2l(y,x)=-atan2l(-y,x).
  *	2. Reduce x to positive by (if x and y are unexceptional):
- *		ARG (x+iy) = arctan(y/x)   	   ... if x > 0,
+ *		ARG (x+iy) = arctan(y/x)	   ... if x > 0,
  *		ARG (x+iy) = pi - arctan[y/(-x)]   ... if x < 0,
  *
  * Special cases:
@@ -43,11 +43,7 @@
 #include "math.h"
 #include "math_private.h"
 
-#ifdef __STDC__
 static const long double
-#else
-static long double
-#endif
 tiny  = 1.0e-300L,
 zero  = 0.0,
 pi_o_4  = 7.85398163397448309615660845819875699e-01L, /* 3ffe921fb54442d18469898cc51701b8 */
@@ -55,12 +51,8 @@ pi_o_2  = 1.57079632679489661923132169163975140e+00L, /* 3fff921fb54442d18469898
 pi      = 3.14159265358979323846264338327950280e+00L, /* 4000921fb54442d18469898cc51701b8 */
 pi_lo   = 8.67181013012378102479704402604335225e-35L; /* 3f8dcd129024e088a67cc74020bbea64 */
 
-#ifdef __STDC__
-	long double __ieee754_atan2l(long double y, long double x)
-#else
-	long double __ieee754_atan2l(y,x)
-	long double  y,x;
-#endif
+long double
+__ieee754_atan2l(long double y, long double x)
 {
 	long double z;
 	int64_t k,m,hx,hy,ix,iy;
@@ -80,7 +72,7 @@ pi_lo   = 8.67181013012378102479704402604335225e-35L; /* 3f8dcd129024e088a67cc74
 	if((iy|(ly&0x7fffffffffffffffLL))==0) {
 	    switch(m) {
 		case 0:
-		case 1: return y; 	/* atan(+-0,+anything)=+-0 */
+		case 1: return y;	/* atan(+-0,+anything)=+-0 */
 		case 2: return  pi+tiny;/* atan(+0,-anything) = pi */
 		case 3: return -pi-tiny;/* atan(-0,-anything) =-pi */
 	    }
@@ -111,14 +103,15 @@ pi_lo   = 8.67181013012378102479704402604335225e-35L; /* 3f8dcd129024e088a67cc74
 
     /* compute y/x */
 	k = (iy-ix)>>52;
-	if(k > 120) z=pi_o_2+0.5L*pi_lo; 	/* |y/x| >  2**120 */
-	else if(hx<0&&k<-120) z=0.0L; 		/* |y|/x < -2**120 */
+	if(k > 120) z=pi_o_2+0.5L*pi_lo;	/* |y/x| >  2**120 */
+	else if(hx<0&&k<-120) z=0.0L;		/* |y|/x < -2**120 */
 	else z=__atanl(fabsl(y/x));		/* safe to do y/x */
 	switch (m) {
 	    case 0: return       z  ;	/* atan(+,+) */
 	    case 1: return      -z  ;	/* atan(-,+) */
 	    case 2: return  pi-(z-pi_lo);/* atan(+,-) */
 	    default: /* case 3 */
-	    	    return  (z-pi_lo)-pi;/* atan(-,-) */
+		    return  (z-pi_lo)-pi;/* atan(-,-) */
 	}
 }
+strong_alias (__ieee754_atan2l, __atan2l_finite)
diff --git a/libc/sysdeps/ieee754/ldbl-128ibm/e_atanhl.c b/libc/sysdeps/ieee754/ldbl-128ibm/e_atanhl.c
index a801bd6f7..c879e4518 100644
--- a/libc/sysdeps/ieee754/ldbl-128ibm/e_atanhl.c
+++ b/libc/sysdeps/ieee754/ldbl-128ibm/e_atanhl.c
@@ -10,10 +10,6 @@
  * ====================================================
  */
 
-#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)
@@ -22,7 +18,7 @@ static char rcsid[] = "$NetBSD: e_atanh.c,v 1.8 1995/05/10 20:44:55 jtc Exp $";
  *	atanh(x) = --- * log(1 + -------) = 0.5 * log1p(2 * --------)
  *                  2             1 - x                      1 - x
  *
- * 	For x<0.5
+ *	For x<0.5
  *	atanh(x) = 0.5*log1p(2x+2x*x/(1-x))
  *
  * Special cases:
@@ -35,24 +31,12 @@ static char rcsid[] = "$NetBSD: e_atanh.c,v 1.8 1995/05/10 20:44:55 jtc Exp $";
 #include "math.h"
 #include "math_private.h"
 
-#ifdef __STDC__
 static const long double one = 1.0L, huge = 1e300L;
-#else
-static long double one = 1.0L, huge = 1e300L;
-#endif
 
-#ifdef __STDC__
 static const long double zero = 0.0L;
-#else
-static long double zero = 0.0L;
-#endif
 
-#ifdef __STDC__
-	long double __ieee754_atanhl(long double x)
-#else
-	long double __ieee754_atanhl(x)
-	long double x;
-#endif
+long double
+__ieee754_atanhl(long double x)
 {
 	long double t;
 	int64_t hx,ix;
@@ -61,7 +45,7 @@ static long double zero = 0.0L;
 	ix = hx&0x7fffffffffffffffLL;
 	if (ix >= 0x3ff0000000000000LL) { /* |x|>=1 */
 	    if (ix > 0x3ff0000000000000LL)
-	    	return (x-x)/(x-x);
+		return (x-x)/(x-x);
 	    t = fabsl (x);
 	    if (t > one)
 		return (x-x)/(x-x);
@@ -77,3 +61,4 @@ static long double zero = 0.0L;
 	    t = 0.5*__log1pl((x+x)/(one-x));
 	if(hx>=0) return t; else return -t;
 }
+strong_alias (__ieee754_atanhl, __atanhl_finite)
diff --git a/libc/sysdeps/ieee754/ldbl-128ibm/e_coshl.c b/libc/sysdeps/ieee754/ldbl-128ibm/e_coshl.c
index 73cb47892..ebc943639 100644
--- a/libc/sysdeps/ieee754/ldbl-128ibm/e_coshl.c
+++ b/libc/sysdeps/ieee754/ldbl-128ibm/e_coshl.c
@@ -10,22 +10,18 @@
  * ====================================================
  */
 
-#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 :
  * 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
+ *							[ 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
+ *							  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,18 +34,10 @@ 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 long double one = 1.0L, half=0.5L, huge = 1.0e300L;
-#else
-static long double one = 1.0L, half=0.5L, huge = 1.0e300L;
-#endif
 
-#ifdef __STDC__
-	long double __ieee754_coshl(long double x)
-#else
-	long double __ieee754_coshl(x)
-	long double x;
-#endif
+long double
+__ieee754_coshl (long double x)
 {
 	long double t,w;
 	int64_t ix;
@@ -79,7 +67,7 @@ static long double one = 1.0L, half=0.5L, huge = 1.0e300L;
 	if (ix < 0x40862e42fefa39efLL)  return half*__ieee754_expl(fabsl(x));
 
     /* |x| in [log(maxdouble), overflowthresold] */
-        if (ix < 0x408633ce8fb9f87dLL) {
+	if (ix < 0x408633ce8fb9f87dLL) {
 	    w = __ieee754_expl(half*fabsl(x));
 	    t = half*w;
 	    return t*w;
@@ -88,3 +76,4 @@ static long double one = 1.0L, half=0.5L, huge = 1.0e300L;
     /* |x| > overflowthresold, cosh(x) overflow */
 	return huge*huge;
 }
+strong_alias (__ieee754_coshl, __coshl_finite)
diff --git a/libc/sysdeps/ieee754/ldbl-128ibm/e_expl.c b/libc/sysdeps/ieee754/ldbl-128ibm/e_expl.c
index daf2cba32..9e03eae23 100644
--- a/libc/sysdeps/ieee754/ldbl-128ibm/e_expl.c
+++ b/libc/sysdeps/ieee754/ldbl-128ibm/e_expl.c
@@ -1,5 +1,5 @@
 /* Quad-precision floating point e^x.
-   Copyright (C) 1999,2004,2006, 2008 Free Software Foundation, Inc.
+   Copyright (C) 1999,2004,2006, 2008, 2011 Free Software Foundation, Inc.
    This file is part of the GNU C Library.
    Contributed by Jakub Jelinek <jj@ultra.linux.cz>
    Partly based on double-precision code
@@ -255,3 +255,4 @@ __ieee754_expl (long double x)
     /* Return x, if x is a NaN or Inf; or overflow, otherwise.  */
     return TWO1023*x;
 }
+strong_alias (__ieee754_expl, __expl_finite)
diff --git a/libc/sysdeps/ieee754/ldbl-128ibm/e_fmodl.c b/libc/sysdeps/ieee754/ldbl-128ibm/e_fmodl.c
index e99b0bac3..4ad59a091 100644
--- a/libc/sysdeps/ieee754/ldbl-128ibm/e_fmodl.c
+++ b/libc/sysdeps/ieee754/ldbl-128ibm/e_fmodl.c
@@ -22,18 +22,10 @@
 #include "math_private.h"
 #include <ieee754.h>
 
-#ifdef __STDC__
 static const long double one = 1.0, Zero[] = {0.0, -0.0,};
-#else
-static long double one = 1.0, Zero[] = {0.0, -0.0,};
-#endif
 
-#ifdef __STDC__
-	long double __ieee754_fmodl(long double x, long double y)
-#else
-	long double __ieee754_fmodl(x,y)
-	long double x,y;
-#endif
+long double
+__ieee754_fmodl (long double x, long double y)
 {
 	int64_t n,hx,hy,hz,ix,iy,sx,i;
 	u_int64_t lx,ly,lz;
@@ -76,8 +68,8 @@ static long double one = 1.0, Zero[] = {0.0, -0.0,};
     /* Make the IBM extended format 105 bit mantissa look like the ieee854 112
        bit mantissa so the following operatations will give the correct
        result.  */
-        ldbl_extract_mantissa(&hx, &lx, &temp, x);
-        ldbl_extract_mantissa(&hy, &ly, &temp, y);
+	ldbl_extract_mantissa(&hx, &lx, &temp, x);
+	ldbl_extract_mantissa(&hy, &ly, &temp, y);
 
     /* set up {hx,lx}, {hy,ly} and align y to x */
 	if(ix >= -1022)
@@ -85,8 +77,8 @@ static long double one = 1.0, Zero[] = {0.0, -0.0,};
 	else {		/* subnormal x, shift x to normal */
 	    n = -1022-ix;
 	    if(n<=63) {
-	        hx = (hx<<n)|(lx>>(64-n));
-	        lx <<= n;
+		hx = (hx<<n)|(lx>>(64-n));
+		lx <<= n;
 	    } else {
 		hx = lx<<(n-64);
 		lx = 0;
@@ -97,8 +89,8 @@ static long double one = 1.0, Zero[] = {0.0, -0.0,};
 	else {		/* subnormal y, shift y to normal */
 	    n = -1022-iy;
 	    if(n<=63) {
-	        hy = (hy<<n)|(ly>>(64-n));
-	        ly <<= n;
+		hy = (hy<<n)|(ly>>(64-n));
+		ly <<= n;
 	    } else {
 		hy = ly<<(n-64);
 		ly = 0;
@@ -111,16 +103,16 @@ static long 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>>63); lx = lx+lx;}
 	    else {
-	    	if((hz|(lz&0x7fffffffffffffff))==0) 		/* return sign(x)*0 */
+		if((hz|(lz&0x7fffffffffffffff))==0)		/* return sign(x)*0 */
 		    return Zero[(u_int64_t)sx>>63];
-	    	hx = hz+hz+(lz>>63); lx = lz+lz;
+		hx = hz+hz+(lz>>63); 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&0x7fffffffffffffff))==0) 			/* return sign(x)*0 */
+	if((hx|(lx&0x7fffffffffffffff))==0)			/* return sign(x)*0 */
 	    return Zero[(u_int64_t)sx>>63];
 	while(hx<0x0001000000000000LL) {	/* normalize x */
 	    hx = hx+hx+(lx>>63); lx = lx+lx;
@@ -143,3 +135,4 @@ static long double one = 1.0, Zero[] = {0.0, -0.0,};
 	}
 	return x;		/* exact output */
 }
+strong_alias (__ieee754_fmodl, __fmodl_finite)
diff --git a/libc/sysdeps/ieee754/ldbl-128ibm/e_gammal_r.c b/libc/sysdeps/ieee754/ldbl-128ibm/e_gammal_r.c
index 03bcb2129..f20ea9e05 100644
--- a/libc/sysdeps/ieee754/ldbl-128ibm/e_gammal_r.c
+++ b/libc/sysdeps/ieee754/ldbl-128ibm/e_gammal_r.c
@@ -1,8 +1,8 @@
 /* Implementation of gamma function according to ISO C.
-   Copyright (C) 1997,1999,2002,2004,2006 Free Software Foundation, Inc.
+   Copyright (C) 1997,1999,2002,2004,2006,2011 Free Software Foundation, Inc.
    This file is part of the GNU C Library.
    Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997 and
-   		  Jakub Jelinek <jj@ultra.linux.cz, 1999.
+		  Jakub Jelinek <jj@ultra.linux.cz, 1999.
 
    The GNU C Library is free software; you can redistribute it and/or
    modify it under the terms of the GNU Lesser General Public
@@ -56,3 +56,4 @@ __ieee754_gammal_r (long double x, int *signgamp)
   /* XXX FIXME.  */
   return __ieee754_expl (__ieee754_lgammal_r (x, signgamp));
 }
+strong_alias (__ieee754_gammal_r, __gammal_r_finite)
diff --git a/libc/sysdeps/ieee754/ldbl-128ibm/e_hypotl.c b/libc/sysdeps/ieee754/ldbl-128ibm/e_hypotl.c
index 4330f28b9..4ef076741 100644
--- a/libc/sysdeps/ieee754/ldbl-128ibm/e_hypotl.c
+++ b/libc/sysdeps/ieee754/ldbl-128ibm/e_hypotl.c
@@ -10,10 +10,6 @@
  * ====================================================
  */
 
-#if defined(LIBM_SCCS) && !defined(lint)
-static char rcsid[] = "$NetBSD: e_hypotl.c,v 1.9 1995/05/12 04:57:27 jtc Exp $";
-#endif
-
 /* __ieee754_hypotl(x,y)
  *
  * Method :
@@ -42,8 +38,8 @@ static char rcsid[] = "$NetBSD: e_hypotl.c,v 1.9 1995/05/12 04:57:27 jtc Exp $";
  *	hypotl(x,y) is NAN if x or y is NAN.
  *
  * Accuracy:
- * 	hypotl(x,y) returns sqrtl(x^2+y^2) with error less
- * 	than 1 ulps (units in the last place)
+ *	hypotl(x,y) returns sqrtl(x^2+y^2) with error less
+ *	than 1 ulps (units in the last place)
  */
 
 #include "math.h"
@@ -52,12 +48,8 @@ static char rcsid[] = "$NetBSD: e_hypotl.c,v 1.9 1995/05/12 04:57:27 jtc Exp $";
 static const long double two600 = 0x1.0p+600L;
 static const long double two1022 = 0x1.0p+1022L;
 
-#ifdef __STDC__
-	long double __ieee754_hypotl(long double x, long double y)
-#else
-	long double __ieee754_hypotl(x,y)
-	long double x, y;
-#endif
+long double
+__ieee754_hypotl(long double x, long double y)
 {
 	long double a,b,t1,t2,y1,y2,w,kld;
 	int64_t j,k,ha,hb;
@@ -93,7 +85,7 @@ static const long double two1022 = 0x1.0p+1022L;
 	}
 	if(hb < 0x20b0000000000000LL) {	/* b < 2**-500 */
 	    if(hb <= 0x000fffffffffffffLL) {	/* subnormal b or 0 */
-	        u_int64_t low;
+		u_int64_t low;
 		GET_LDOUBLE_LSW64(low,b);
 		if((hb|(low&0x7fffffffffffffffLL))==0) return a;
 		t1=two1022;	/* t1=2^1022 */
@@ -102,7 +94,7 @@ static const long double two1022 = 0x1.0p+1022L;
 		k -= 1022;
 		kld = kld / two1022;
 	    } else {		/* scale a and b by 2^600 */
-	        ha += 0x2580000000000000LL; 	/* a *= 2^600 */
+		ha += 0x2580000000000000LL;	/* a *= 2^600 */
 		hb += 0x2580000000000000LL;	/* b *= 2^600 */
 		k -= 600;
 		a *= two600;
@@ -129,3 +121,4 @@ static const long double two1022 = 0x1.0p+1022L;
 	else
 	    return w;
 }
+strong_alias (__ieee754_hypotl, __hypotl_finite)
diff --git a/libc/sysdeps/ieee754/ldbl-128ibm/e_jnl.c b/libc/sysdeps/ieee754/ldbl-128ibm/e_jnl.c
index 0eea74567..2114753f8 100644
--- a/libc/sysdeps/ieee754/ldbl-128ibm/e_jnl.c
+++ b/libc/sysdeps/ieee754/ldbl-128ibm/e_jnl.c
@@ -286,7 +286,16 @@ __ieee754_jnl (n, x)
 		    }
 		}
 	    }
-	  b = (t * __ieee754_j0l (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_j0l (x);
+	  w = __ieee754_j1l (x);
+	  if (fabsl (z) >= fabsl (w))
+	    b = (t * z / b);
+	  else
+	    b = (t * w / a);
 	}
     }
   if (sgn == 1)
@@ -294,6 +303,7 @@ __ieee754_jnl (n, x)
   else
     return b;
 }
+strong_alias (__ieee754_jnl, __jnl_finite)
 
 #ifdef __STDC__
 long double
@@ -400,3 +410,4 @@ __ieee754_ynl (n, x)
   else
     return -b;
 }
+strong_alias (__ieee754_ynl, __ynl_finite)
diff --git a/libc/sysdeps/ieee754/ldbl-128ibm/e_log10l.c b/libc/sysdeps/ieee754/ldbl-128ibm/e_log10l.c
index 27e2c71b8..de57be39d 100644
--- a/libc/sysdeps/ieee754/ldbl-128ibm/e_log10l.c
+++ b/libc/sysdeps/ieee754/ldbl-128ibm/e_log10l.c
@@ -179,8 +179,7 @@ deval (long double x, const long double *p, int n)
 
 
 long double
-__ieee754_log10l (x)
-     long double x;
+__ieee754_log10l (long double x)
 {
   long double z;
   long double y;
@@ -256,3 +255,4 @@ done:
   z += e * L102A;
   return (z);
 }
+strong_alias (__ieee754_log10l, __log10l_finite)
diff --git a/libc/sysdeps/ieee754/ldbl-128ibm/e_log2l.c b/libc/sysdeps/ieee754/ldbl-128ibm/e_log2l.c
index fe8a8e1d6..9737e13f0 100644
--- a/libc/sysdeps/ieee754/ldbl-128ibm/e_log2l.c
+++ b/libc/sysdeps/ieee754/ldbl-128ibm/e_log2l.c
@@ -248,3 +248,4 @@ done:
   z += e;
   return (z);
 }
+strong_alias (__ieee754_log2l, __log2l_finite)
diff --git a/libc/sysdeps/ieee754/ldbl-128ibm/e_logl.c b/libc/sysdeps/ieee754/ldbl-128ibm/e_logl.c
index aa5fc3740..683260806 100644
--- a/libc/sysdeps/ieee754/ldbl-128ibm/e_logl.c
+++ b/libc/sysdeps/ieee754/ldbl-128ibm/e_logl.c
@@ -253,7 +253,7 @@ __ieee754_logl(long double x)
       /* log(u) = log( t u/t ) = log(t) + log(u/t)
 	 log(t) is tabulated in the lookup table.
 	 Express log(u/t) = log(1+z),  where z = u/t - 1 = (u-t)/t.
-         cf. Cody & Waite. */
+	 cf. Cody & Waite. */
       z = (u.value - t.value) / t.value;
     }
   /* Series expansion of log(1+z).  */
@@ -279,3 +279,4 @@ __ieee754_logl(long double x)
   y += e * ln2a;
   return y;
 }
+strong_alias (__ieee754_logl, __logl_finite)
diff --git a/libc/sysdeps/ieee754/ldbl-128ibm/e_powl.c b/libc/sysdeps/ieee754/ldbl-128ibm/e_powl.c
index feeaa8ce2..9b1f2be1d 100644
--- a/libc/sysdeps/ieee754/ldbl-128ibm/e_powl.c
+++ b/libc/sysdeps/ieee754/ldbl-128ibm/e_powl.c
@@ -144,14 +144,8 @@ static const long double
   cp_h = 9.6179669392597555432899980587535537779331E-1L,
   cp_l = 5.0577616648125906047157785230014751039424E-17L;
 
-#ifdef __STDC__
 long double
 __ieee754_powl (long double x, long double y)
-#else
-long double
-__ieee754_powl (x, y)
-     long double x, y;
-#endif
 {
   long double z, ax, z_h, z_l, p_h, p_l;
   long double y1, t1, t2, r, s, t, u, v, w;
@@ -390,7 +384,7 @@ __ieee754_powl (x, y)
     {
       /* if z > 16384 */
       if (((j - 0x40d00000) | o.parts32.w1
-        | (o.parts32.w2 & 0x7fffffff) | o.parts32.w3) != 0)
+	| (o.parts32.w2 & 0x7fffffff) | o.parts32.w3) != 0)
 	return s * huge * huge;	/* overflow */
       else
 	{
@@ -402,7 +396,7 @@ __ieee754_powl (x, y)
     {
       /* z < -16495 */
       if (((j - 0xc0d01bc0) | o.parts32.w1
-         | (o.parts32.w2 & 0x7fffffff) | o.parts32.w3) != 0)
+	 | (o.parts32.w2 & 0x7fffffff) | o.parts32.w3) != 0)
 	return s * tiny * tiny;	/* underflow */
       else
 	{
@@ -439,3 +433,4 @@ __ieee754_powl (x, y)
   z = __scalbnl (z, n);
   return s * z;
 }
+strong_alias (__ieee754_powl, __powl_finite)
diff --git a/libc/sysdeps/ieee754/ldbl-128ibm/e_remainderl.c b/libc/sysdeps/ieee754/ldbl-128ibm/e_remainderl.c
index b7fa68f32..d4a847dbe 100644
--- a/libc/sysdeps/ieee754/ldbl-128ibm/e_remainderl.c
+++ b/libc/sysdeps/ieee754/ldbl-128ibm/e_remainderl.c
@@ -14,8 +14,8 @@
 
 /* __ieee754_remainderl(x,p)
  * Return :
- * 	returns  x REM p  =  x - [x/p]*p as if in infinite
- * 	precise arithmetic, where [x/p] is the (infinite bit)
+ *	returns  x REM p  =  x - [x/p]*p as if in infinite
+ *	precise arithmetic, where [x/p] is the (infinite bit)
  *	integer nearest x/p (in half way case choose the even one).
  * Method :
  *	Based on fmodl() return x-[x/p]chopped*p exactlp.
@@ -24,19 +24,11 @@
 #include "math.h"
 #include "math_private.h"
 
-#ifdef __STDC__
 static const long double zero = 0.0L;
-#else
-static long double zero = 0.0L;
-#endif
 
 
-#ifdef __STDC__
-	long double __ieee754_remainderl(long double x, long double p)
-#else
-	long double __ieee754_remainderl(x,p)
-	long double x,p;
-#endif
+long double
+__ieee754_remainderl(long double x, long double p)
 {
 	int64_t hx,hp;
 	u_int64_t sx,lx,lp;
@@ -49,7 +41,7 @@ static long double zero = 0.0L;
 	hx &= 0x7fffffffffffffffLL;
 
     /* purge off exception values */
-	if((hp|(lp&0x7fffffffffffffff))==0) return (x*p)/(x*p); 	/* p = 0 */
+	if((hp|(lp&0x7fffffffffffffff))==0) return (x*p)/(x*p);	/* p = 0 */
 	if((hx>=0x7ff0000000000000LL)||			/* x not finite */
 	  ((hp>=0x7ff0000000000000LL)&&			/* p is NaN */
 	  (((hp-0x7ff0000000000000LL)|lp)!=0)))
@@ -76,3 +68,4 @@ static long double zero = 0.0L;
 	SET_LDOUBLE_MSW64(x,hx^sx);
 	return x;
 }
+strong_alias (__ieee754_remainderl, __remainderl_finite)
diff --git a/libc/sysdeps/ieee754/ldbl-128ibm/e_sinhl.c b/libc/sysdeps/ieee754/ldbl-128ibm/e_sinhl.c
index 38ae71d4b..b8e86b1a6 100644
--- a/libc/sysdeps/ieee754/ldbl-128ibm/e_sinhl.c
+++ b/libc/sysdeps/ieee754/ldbl-128ibm/e_sinhl.c
@@ -10,18 +10,14 @@
  */
 
 
-#if defined(LIBM_SCCS) && !defined(lint)
-static char rcsid[] = "$NetBSD: e_sinh.c,v 1.7 1995/05/10 20:46:13 jtc Exp $";
-#endif
-
 /* __ieee754_sinh(x)
  * 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)
+ *						    E + E/(E+1)
  *	    0        <= x <= 22     :  sinh(x) := --------------, E=expm1(x)
- *			       			        2
+ *							2
  *
  *	    22       <= x <= lnovft :  sinh(x) := exp(x)/2
  *	    lnovft   <= x <= ln2ovft:  sinh(x) := exp(x/2)/2 * exp(x/2)
@@ -35,18 +31,10 @@ 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 long double one = 1.0, shuge = 1.0e307;
-#else
-static long double one = 1.0, shuge = 1.0e307;
-#endif
 
-#ifdef __STDC__
-	long double __ieee754_sinhl(long double x)
-#else
-	long double __ieee754_sinhl(x)
-	long double x;
-#endif
+long double
+__ieee754_sinhl(long double x)
 {
 	long double t,w,h;
 	int64_t ix,jx;
@@ -62,7 +50,7 @@ static long double one = 1.0, shuge = 1.0e307;
 	if (jx<0) h = -h;
     /* |x| in [0,22], return sign(x)*0.5*(E+E/(E+1))) */
 	if (ix < 0x4036000000000000LL) {	/* |x|<22 */
-	    if (ix<0x3e20000000000000LL) 	/* |x|<2**-29 */
+	    if (ix<0x3e20000000000000LL)	/* |x|<2**-29 */
 		if(shuge+x>one) return x;/* sinhl(tiny) = tiny with inexact */
 	    t = __expm1l(fabsl(x));
 	    if(ix<0x3ff0000000000000LL) return h*(2.0*t-t*t/(t+one));
@@ -82,3 +70,4 @@ static long double one = 1.0, shuge = 1.0e307;
     /* |x| > overflowthresold, sinh(x) overflow */
 	return x*shuge;
 }
+strong_alias (__ieee754_sinhl, __sinhl_finite)
diff --git a/libc/sysdeps/ieee754/ldbl-128ibm/e_sqrtl.c b/libc/sysdeps/ieee754/ldbl-128ibm/e_sqrtl.c
index fe6bb55b0..68aa18f55 100644
--- a/libc/sysdeps/ieee754/ldbl-128ibm/e_sqrtl.c
+++ b/libc/sysdeps/ieee754/ldbl-128ibm/e_sqrtl.c
@@ -1,7 +1,7 @@
 /*
  * IBM Accurate Mathematical Library
  * written by International Business Machines Corp.
- * Copyright (C) 2001, 2004, 2006 Free Software Foundation
+ * Copyright (C) 2001, 2004, 2006, 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
@@ -50,7 +50,7 @@ twom54 = 5.55111512312578270212e-17; /* 0x3C90000000000000 */
 /* it computes the correctly rounded (to nearest) value of square    */
 /* root of x.                                                        */
 /*********************************************************************/
-long double __ieee754_sqrtl(long double x) 
+long double __ieee754_sqrtl(long double x)
 {
   static const long double big = 134217728.0, big1 = 134217729.0;
   long double t,s,i;
@@ -107,3 +107,4 @@ long double __ieee754_sqrtl(long double x)
     return tm256.x*__ieee754_sqrtl(x*t512.x);
   }
 }
+strong_alias (__ieee754_sqrtl, __sqrtl_finite)
diff --git a/libc/sysdeps/ieee754/ldbl-128ibm/s_atanl.c b/libc/sysdeps/ieee754/ldbl-128ibm/s_atanl.c
index b6195f10b..db31e4f90 100644
--- a/libc/sysdeps/ieee754/ldbl-128ibm/s_atanl.c
+++ b/libc/sysdeps/ieee754/ldbl-128ibm/s_atanl.c
@@ -59,6 +59,7 @@
     Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307  USA */
 
 
+#include <math.h>
 #include "math_private.h"
 #include <math_ldbl_opt.h>
 
diff --git a/libc/sysdeps/ieee754/ldbl-128ibm/s_fmal.c b/libc/sysdeps/ieee754/ldbl-128ibm/s_fmal.c
new file mode 100644
index 000000000..e59d2c349
--- /dev/null
+++ b/libc/sysdeps/ieee754/ldbl-128ibm/s_fmal.c
@@ -0,0 +1,39 @@
+/* Compute x * y + z as ternary operation.
+   Copyright (C) 2011 Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
+   Contributed by David Flaherty <flaherty@linux.vnet.ibm.com>.
+
+   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_ldbl_opt.h>
+
+long double
+__fmal (long double x, long double y, long double z)
+{
+	/* An IBM long double 128 is really just 2 IEEE64 doubles, and in
+	 * the case of inf/nan only the first double counts. So we use the
+	 * (double) cast to avoid any data movement.   */
+       if ((finite ((double)x) && finite ((double)y)) && isinf ((double)z))
+               return (z);
+
+       return (x * y) + z;
+}
+#ifdef IS_IN_libm
+long_double_symbol (libm, __fmal, fmal);
+#else
+long_double_symbol (libc, __fmal, fmal);
+#endif
diff --git a/libc/sysdeps/ieee754/ldbl-128ibm/s_isinf_nsl.c b/libc/sysdeps/ieee754/ldbl-128ibm/s_isinf_nsl.c
new file mode 100644
index 000000000..edeaba5f7
--- /dev/null
+++ b/libc/sysdeps/ieee754/ldbl-128ibm/s_isinf_nsl.c
@@ -0,0 +1,16 @@
+/*
+ * __isinf_nsl(x) returns != 0 if x is ±inf, else 0;
+ * no branching!
+ */
+
+#include "math.h"
+#include "math_private.h"
+
+int
+__isinf_nsl (long double x)
+{
+	int64_t hx,lx;
+	GET_LDOUBLE_WORDS64(hx,lx,x);
+	return !((lx & 0x7fffffffffffffffLL)
+		 | ((hx & 0x7fffffffffffffffLL) ^ 0x7ff0000000000000LL));
+}