Merge changes between r22552 and r22663 from /fsf/trunk.

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

11 files changed, 132 insertions, 74 deletions

diff --git a/libc/sysdeps/ieee754/bits/nan.h b/libc/sysdeps/ieee754/bits/nan.h
index 935271a7c..41f47ba09 100644
--- a/libc/sysdeps/ieee754/bits/nan.h
+++ b/libc/sysdeps/ieee754/bits/nan.h
@@ -39,14 +39,14 @@
 # include <endian.h>
 
 # if __BYTE_ORDER == __BIG_ENDIAN
-#  define __nan_bytes		{ 0x7f, 0xc0, 0, 0 }
+#  define __qnan_bytes		{ 0x7f, 0xc0, 0, 0 }
 # endif
 # if __BYTE_ORDER == __LITTLE_ENDIAN
-#  define __nan_bytes		{ 0, 0, 0xc0, 0x7f }
+#  define __qnan_bytes		{ 0, 0, 0xc0, 0x7f }
 # endif
 
-static union { unsigned char __c[4]; float __d; } __nan_union
-  __attribute__ ((__unused__)) = { __nan_bytes };
-# define NAN	(__nan_union.__d)
+static union { unsigned char __c[4]; float __d; } __qnan_union
+  __attribute__ ((__unused__)) = { __qnan_bytes };
+# define NAN	(__qnan_union.__d)
 
 #endif	/* GCC.  */
diff --git a/libc/sysdeps/ieee754/dbl-64/e_j0.c b/libc/sysdeps/ieee754/dbl-64/e_j0.c
index f393a762b..d641a0914 100644
--- a/libc/sysdeps/ieee754/dbl-64/e_j0.c
+++ b/libc/sysdeps/ieee754/dbl-64/e_j0.c
@@ -293,7 +293,8 @@ pzero(double x)
 	int32_t ix;
 	GET_HIGH_WORD(ix,x);
 	ix &= 0x7fffffff;
-	if(ix>=0x40200000)     {p = pR8; q= pS8;}
+	if (ix>=0x41b00000)    {return one;}
+	else 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;}
@@ -400,7 +401,8 @@ qzero(double x)
 	int32_t ix;
 	GET_HIGH_WORD(ix,x);
 	ix &= 0x7fffffff;
-	if(ix>=0x40200000)     {p = qR8; q= qS8;}
+	if (ix>=0x41b00000)    {return -.125/x;}
+	else if(ix>=0x40200000){p = qR8; q= qS8;}
 	else if(ix>=0x40122E8B){p = qR5; q= qS5;}
 	else if(ix>=0x4006DB6D){p = qR3; q= qS3;}
 	else if(ix>=0x40000000){p = qR2; q= qS2;}
diff --git a/libc/sysdeps/ieee754/dbl-64/e_j1.c b/libc/sysdeps/ieee754/dbl-64/e_j1.c
index cba4d46b1..cca5f20b4 100644
--- a/libc/sysdeps/ieee754/dbl-64/e_j1.c
+++ b/libc/sysdeps/ieee754/dbl-64/e_j1.c
@@ -291,7 +291,8 @@ pone(double x)
 	int32_t ix;
 	GET_HIGH_WORD(ix,x);
 	ix &= 0x7fffffff;
-	if(ix>=0x40200000)     {p = pr8; q= ps8;}
+	if (ix>=0x41b00000)    {return one;}
+	else 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;}
@@ -399,7 +400,8 @@ qone(double x)
 	int32_t ix;
 	GET_HIGH_WORD(ix,x);
 	ix &= 0x7fffffff;
-	if(ix>=0x40200000)     {p = qr8; q= qs8;}
+	if (ix>=0x41b00000)    {return .375/x;}
+	else if(ix>=0x40200000){p = qr8; q= qs8;}
 	else if(ix>=0x40122E8B){p = qr5; q= qs5;}
 	else if(ix>=0x4006DB6D){p = qr3; q= qs3;}
 	else if(ix>=0x40000000){p = qr2; q= qs2;}
diff --git a/libc/sysdeps/ieee754/dbl-64/mpa.c b/libc/sysdeps/ieee754/dbl-64/mpa.c
index 8fc2626f7..076647654 100644
--- a/libc/sysdeps/ieee754/dbl-64/mpa.c
+++ b/libc/sysdeps/ieee754/dbl-64/mpa.c
@@ -611,6 +611,7 @@ __sub (const mp_no *x, const mp_no *y, mp_no *z, int p)
     }
 }
 
+#ifndef NO__MUL
 /* Multiply *X and *Y and store result in *Z.  X and Y may overlap but not X
    and Z or Y and Z.  For P in [1, 2, 3], the exact result is truncated to P
    digits.  In case P > 3 the error is bounded by 1.001 ULP.  */
@@ -761,7 +762,9 @@ __mul (const mp_no *x, const mp_no *y, mp_no *z, int p)
   EZ = e;
   Z[0] = X[0] * Y[0];
 }
+#endif
 
+#ifndef NO__SQR
 /* Square *X and store result in *Y.  X and Y may not overlap.  For P in
    [1, 2, 3], the exact result is truncated to P digits.  In case P > 3 the
    error is bounded by 1.001 ULP.  This is a faster special case of
@@ -862,6 +865,7 @@ __sqr (const mp_no *x, mp_no *y, int p)
 
   EY = e;
 }
+#endif
 
 /* Invert *X and store in *Y.  Relative error bound:
    - For P = 2: 1.001 * R ^ (1 - P)
diff --git a/libc/sysdeps/ieee754/dbl-64/slowexp.c b/libc/sysdeps/ieee754/dbl-64/slowexp.c
index c423fc311..8f353f634 100644
--- a/libc/sysdeps/ieee754/dbl-64/slowexp.c
+++ b/libc/sysdeps/ieee754/dbl-64/slowexp.c
@@ -27,20 +27,23 @@
 /*Converting from double precision to Multi-precision and calculating     */
 /* e^x                                                                    */
 /**************************************************************************/
-#include "mpa.h"
 #include <math_private.h>
 
+#ifndef USE_LONG_DOUBLE_FOR_MP
+# include "mpa.h"
+void __mpexp (mp_no *x, mp_no *y, int p);
+#endif
+
 #ifndef SECTION
 # define SECTION
 #endif
 
-void __mpexp (mp_no *x, mp_no *y, int p);
-
 /*Converting from double precision to Multi-precision and calculating  e^x */
 double
 SECTION
 __slowexp (double x)
 {
+#ifndef USE_LONG_DOUBLE_FOR_MP
   double w, z, res, eps = 3.0e-26;
   int p;
   mp_no mpx, mpy, mpz, mpw, mpeps, mpcor;
@@ -66,4 +69,7 @@ __slowexp (double x)
       __mp_dbl (&mpy, &res, p);
       return res;
     }
+#else
+  return (double) __ieee754_expl((long double)x);
+#endif
 }
diff --git a/libc/sysdeps/ieee754/dbl-64/slowpow.c b/libc/sysdeps/ieee754/dbl-64/slowpow.c
index cccc7e32c..a379728b1 100644
--- a/libc/sysdeps/ieee754/dbl-64/slowpow.c
+++ b/libc/sysdeps/ieee754/dbl-64/slowpow.c
@@ -59,6 +59,23 @@ __slowpow (double x, double y, double z)
   if (res >= 0)
     return res;
 
+  /* Compute pow as long double.  This is currently only used by powerpc, where
+     one may get 106 bits of accuracy.  */
+#ifdef USE_LONG_DOUBLE_FOR_MP
+  long double ldw, ldz, ldpp;
+  static const long double ldeps = 0x4.0p-96;
+
+  ldz = __ieee754_logl ((long double) x);
+  ldw = (long double) y *ldz;
+  ldpp = __ieee754_expl (ldw);
+  res = (double) (ldpp + ldeps);
+  res1 = (double) (ldpp - ldeps);
+
+  /* Return the result if it is accurate enough.  */
+  if (res == res1)
+    return res;
+#endif
+
   /* Or else, calculate using multiple precision.  P = 10 implies accuracy of
      240 bits accuracy, since MP_NO has a radix of 2^24.  */
   p = 10;
diff --git a/libc/sysdeps/ieee754/ldbl-128/e_j0l.c b/libc/sysdeps/ieee754/ldbl-128/e_j0l.c
index 1b1828958..9e7880c49 100644
--- a/libc/sysdeps/ieee754/ldbl-128/e_j0l.c
+++ b/libc/sysdeps/ieee754/ldbl-128/e_j0l.c
@@ -700,6 +700,25 @@ __ieee754_j0l (long double x)
       return p;
     }
 
+  /* X = x - pi/4
+     cos(X) = cos(x) cos(pi/4) + sin(x) sin(pi/4)
+     = 1/sqrt(2) * (cos(x) + sin(x))
+     sin(X) = sin(x) cos(pi/4) - cos(x) sin(pi/4)
+     = 1/sqrt(2) * (sin(x) - cos(x))
+     sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
+     cf. Fdlibm.  */
+  __sincosl (xx, &s, &c);
+  ss = s - c;
+  cc = s + c;
+  z = -__cosl (xx + xx);
+  if ((s * c) < 0)
+    cc = z / ss;
+  else
+    ss = z / cc;
+
+  if (xx > 0x1p256L)
+    return ONEOSQPI * cc / __ieee754_sqrtl (xx);
+
   xinv = 1.0L / xx;
   z = xinv * xinv;
   if (xinv <= 0.25)
@@ -761,21 +780,6 @@ __ieee754_j0l (long double x)
   p = 1.0L + z * p;
   q = z * xinv * q;
   q = q - 0.125L * xinv;
-  /* X = x - pi/4
-     cos(X) = cos(x) cos(pi/4) + sin(x) sin(pi/4)
-     = 1/sqrt(2) * (cos(x) + sin(x))
-     sin(X) = sin(x) cos(pi/4) - cos(x) sin(pi/4)
-     = 1/sqrt(2) * (sin(x) - cos(x))
-     sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
-     cf. Fdlibm.  */
-  __sincosl (xx, &s, &c);
-  ss = s - c;
-  cc = s + c;
-  z = -__cosl (xx + xx);
-  if ((s * c) < 0)
-    cc = z / ss;
-  else
-    ss = z / cc;
   z = ONEOSQPI * (p * cc - q * ss) / __ieee754_sqrtl (xx);
   return z;
 }
@@ -843,6 +847,25 @@ long double
       return p;
     }
 
+  /* X = x - pi/4
+     cos(X) = cos(x) cos(pi/4) + sin(x) sin(pi/4)
+     = 1/sqrt(2) * (cos(x) + sin(x))
+     sin(X) = sin(x) cos(pi/4) - cos(x) sin(pi/4)
+     = 1/sqrt(2) * (sin(x) - cos(x))
+     sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
+     cf. Fdlibm.  */
+  __sincosl (x, &s, &c);
+  ss = s - c;
+  cc = s + c;
+  z = -__cosl (x + x);
+  if ((s * c) < 0)
+    cc = z / ss;
+  else
+    ss = z / cc;
+
+  if (xx > 0x1p256L)
+    return ONEOSQPI * ss / __ieee754_sqrtl (x);
+
   xinv = 1.0L / xx;
   z = xinv * xinv;
   if (xinv <= 0.25)
@@ -904,21 +927,6 @@ long double
   p = 1.0L + z * p;
   q = z * xinv * q;
   q = q - 0.125L * xinv;
-  /* X = x - pi/4
-     cos(X) = cos(x) cos(pi/4) + sin(x) sin(pi/4)
-     = 1/sqrt(2) * (cos(x) + sin(x))
-     sin(X) = sin(x) cos(pi/4) - cos(x) sin(pi/4)
-     = 1/sqrt(2) * (sin(x) - cos(x))
-     sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
-     cf. Fdlibm.  */
-  __sincosl (x, &s, &c);
-  ss = s - c;
-  cc = s + c;
-  z = -__cosl (x + x);
-  if ((s * c) < 0)
-    cc = z / ss;
-  else
-    ss = z / cc;
   z = ONEOSQPI * (p * ss + q * cc) / __ieee754_sqrtl (x);
   return z;
 }
diff --git a/libc/sysdeps/ieee754/ldbl-128/e_j1l.c b/libc/sysdeps/ieee754/ldbl-128/e_j1l.c
index f16343b26..95e01a39c 100644
--- a/libc/sysdeps/ieee754/ldbl-128/e_j1l.c
+++ b/libc/sysdeps/ieee754/ldbl-128/e_j1l.c
@@ -706,6 +706,29 @@ __ieee754_j1l (long double x)
       return p;
     }
 
+  /* X = x - 3 pi/4
+     cos(X) = cos(x) cos(3 pi/4) + sin(x) sin(3 pi/4)
+     = 1/sqrt(2) * (-cos(x) + sin(x))
+     sin(X) = sin(x) cos(3 pi/4) - cos(x) sin(3 pi/4)
+     = -1/sqrt(2) * (sin(x) + cos(x))
+     cf. Fdlibm.  */
+  __sincosl (xx, &s, &c);
+  ss = -s - c;
+  cc = s - c;
+  z = __cosl (xx + xx);
+  if ((s * c) > 0)
+    cc = z / ss;
+  else
+    ss = z / cc;
+
+  if (xx > 0x1p256L)
+    {
+      z = ONEOSQPI * cc / __ieee754_sqrtl (xx);
+      if (x < 0)
+	z = -z;
+      return z;
+    }
+
   xinv = 1.0L / xx;
   z = xinv * xinv;
   if (xinv <= 0.25)
@@ -767,20 +790,6 @@ __ieee754_j1l (long double x)
   p = 1.0L + z * p;
   q = z * q;
   q = q * xinv + 0.375L * xinv;
-  /* X = x - 3 pi/4
-     cos(X) = cos(x) cos(3 pi/4) + sin(x) sin(3 pi/4)
-     = 1/sqrt(2) * (-cos(x) + sin(x))
-     sin(X) = sin(x) cos(3 pi/4) - cos(x) sin(3 pi/4)
-     = -1/sqrt(2) * (sin(x) + cos(x))
-     cf. Fdlibm.  */
-  __sincosl (xx, &s, &c);
-  ss = -s - c;
-  cc = s - c;
-  z = __cosl (xx + xx);
-  if ((s * c) > 0)
-    cc = z / ss;
-  else
-    ss = z / cc;
   z = ONEOSQPI * (p * cc - q * ss) / __ieee754_sqrtl (xx);
   if (x < 0)
     z = -z;
@@ -850,6 +859,24 @@ __ieee754_y1l (long double x)
       return p;
     }
 
+  /* X = x - 3 pi/4
+     cos(X) = cos(x) cos(3 pi/4) + sin(x) sin(3 pi/4)
+     = 1/sqrt(2) * (-cos(x) + sin(x))
+     sin(X) = sin(x) cos(3 pi/4) - cos(x) sin(3 pi/4)
+     = -1/sqrt(2) * (sin(x) + cos(x))
+     cf. Fdlibm.  */
+  __sincosl (xx, &s, &c);
+  ss = -s - c;
+  cc = s - c;
+  z = __cosl (xx + xx);
+  if ((s * c) > 0)
+    cc = z / ss;
+  else
+    ss = z / cc;
+
+  if (xx > 0x1p256L)
+    return ONEOSQPI * ss / __ieee754_sqrtl (xx);
+
   xinv = 1.0L / xx;
   z = xinv * xinv;
   if (xinv <= 0.25)
@@ -911,20 +938,6 @@ __ieee754_y1l (long double x)
   p = 1.0L + z * p;
   q = z * q;
   q = q * xinv + 0.375L * xinv;
-  /* X = x - 3 pi/4
-     cos(X) = cos(x) cos(3 pi/4) + sin(x) sin(3 pi/4)
-     = 1/sqrt(2) * (-cos(x) + sin(x))
-     sin(X) = sin(x) cos(3 pi/4) - cos(x) sin(3 pi/4)
-     = -1/sqrt(2) * (sin(x) + cos(x))
-     cf. Fdlibm.  */
-  __sincosl (xx, &s, &c);
-  ss = -s - c;
-  cc = s - c;
-  z = __cosl (xx + xx);
-  if ((s * c) > 0)
-    cc = z / ss;
-  else
-    ss = z / cc;
   z = ONEOSQPI * (p * ss + q * cc) / __ieee754_sqrtl (xx);
   return z;
 }
diff --git a/libc/sysdeps/ieee754/ldbl-128ibm/e_acoshl.c b/libc/sysdeps/ieee754/ldbl-128ibm/e_acoshl.c
index 117bd0f05..abc78a35b 100644
--- a/libc/sysdeps/ieee754/ldbl-128ibm/e_acoshl.c
+++ b/libc/sysdeps/ieee754/ldbl-128ibm/e_acoshl.c
@@ -52,7 +52,7 @@ __ieee754_acoshl(long double x)
 	    return __ieee754_logl(2.0*x-one/(x+__ieee754_sqrtl(t-one)));
 	} else {			/* 1<x<2 */
 	    t = x-one;
-	    return __log1p(t+__sqrtl(2.0*t+t*t));
+	    return __log1p(t+__ieee754_sqrtl(2.0*t+t*t));
 	}
 }
 strong_alias (__ieee754_acoshl, __acoshl_finite)
diff --git a/libc/sysdeps/ieee754/ldbl-128ibm/math_ldbl.h b/libc/sysdeps/ieee754/ldbl-128ibm/math_ldbl.h
index be9ac71cb..1cce1fc4d 100644
--- a/libc/sysdeps/ieee754/ldbl-128ibm/math_ldbl.h
+++ b/libc/sysdeps/ieee754/ldbl-128ibm/math_ldbl.h
@@ -125,7 +125,7 @@ ldbl_insert_mantissa (int sign, int exp, int64_t hi64, u_int64_t lo64)
 /* Handy utility functions to pack/unpack/cononicalize and find the nearbyint
    of long double implemented as double double.  */
 static inline long double
-ldbl_pack (double a, double aa)
+default_ldbl_pack (double a, double aa)
 {
   union ibm_extended_long_double u;
   u.dd[0] = a;
@@ -134,7 +134,7 @@ ldbl_pack (double a, double aa)
 }
 
 static inline void
-ldbl_unpack (long double l, double *a, double *aa)
+default_ldbl_unpack (long double l, double *a, double *aa)
 {
   union ibm_extended_long_double u;
   u.d = l;
@@ -142,6 +142,12 @@ ldbl_unpack (long double l, double *a, double *aa)
   *aa = u.dd[1];
 }
 
+#ifndef ldbl_pack
+# define ldbl_pack   default_ldbl_pack
+#endif
+#ifndef ldbl_unpack
+# define ldbl_unpack default_ldbl_unpack
+#endif
 
 /* Convert a finite long double to canonical form.
    Does not handle +/-Inf properly.  */
diff --git a/libc/sysdeps/ieee754/ldbl-96/e_j1l.c b/libc/sysdeps/ieee754/ldbl-96/e_j1l.c
index 785c0b067..4c13018ae 100644
--- a/libc/sysdeps/ieee754/ldbl-96/e_j1l.c
+++ b/libc/sysdeps/ieee754/ldbl-96/e_j1l.c
@@ -203,7 +203,7 @@ __ieee754_y1l (long double x)
       __sincosl (x, &s, &c);
       ss = -s - c;
       cc = s - c;
-      if (ix < 0x7fe00000)
+      if (ix < 0x7ffe)
 	{			/* make sure x+x not overflow */
 	  z = __cosl (x + x);
 	  if ((s * c) > zero)