atanq.c (atanq): Update from GLIBC.

2012-11-22  David S. Miller  <davem@davemloft.net>
            Tobias Burnus  <burnus@net-b.de>
            Joseph Myers  <joseph@codesourcery.com>

        * math/atanq.c (atanq): Update from GLIBC. Handle tiny and
        very large arguments properly.
        * math/j0q.c (y0q): Update from GLIBC. Avoid arithmetic
        underflow when 'x' is very small.
        * math/j1q.c (y1q): Ditto.
        * math/log1pq.c (log1pq): Update from GLIBC. Saturate
        nonzero exponents with absolute value below 0x1p-128 to
        +/- 0x1p-128.
        * math/powq.c (powq): Update from GLIBC. If xm1 is
        smaller than LDBL_EPSILON/2.0L, just return xm1.


Co-Authored-By: Joseph Myers <joseph@codesourcery.com>
Co-Authored-By: Tobias Burnus <burnus@net-b.de>

From-SVN: r193716

libquadmath/ChangeLog

@@ -1,3 +1,18 @@
+2012-11-22  David S. Miller  <davem@davemloft.net>
+	    Tobias Burnus  <burnus@net-b.de>
+	    Joseph Myers  <joseph@codesourcery.com>
+
+	* math/atanq.c (atanq): Update from GLIBC. Handle tiny and
+	very large arguments properly.
+	* math/j0q.c (y0q): Update from GLIBC. Avoid arithmetic
+	underflow when 'x' is very small.
+	* math/j1q.c (y1q): Ditto.
+	* math/log1pq.c (log1pq): Update from GLIBC. Saturate
+	nonzero exponents with absolute value below 0x1p-128 to
+	+/- 0x1p-128.
+	* math/powq.c (powq): Update from GLIBC. If xm1 is
+	smaller than LDBL_EPSILON/2.0L, just return xm1.
+
 2012-11-21  Tobias Burnus  <burnus@net-b.de>
 
 	PR libquadmath/55225

libquadmath/math/atanq.c

@@ -167,6 +167,7 @@ static const __float128
   q4 = 2.173623741810414221251136181221172551416E1Q;
   /* q5 = 1.000000000000000000000000000000000000000E0 */
 
+static const long double huge = 1.0e4930Q;
 
 __float128
 atanq (__float128 x)
@@ -197,6 +198,22 @@ atanq (__float128 x)
 	return atantbl[83];
     }
 
+  if (k <= 0x3fc50000) /* |x| < 2**-58 */
+    {
+      /* Raise inexact. */
+      if (huge + x > 0.0)
+	return x;
+    }
+
+  if (k >= 0x40720000) /* |x| > 2**115 */
+    {
+      /* Saturate result to {-,+}pi/2 */
+      if (sign)
+	return -atantbl[83];
+      else
+	return atantbl[83];
+    }
+
   if (sign)
       x = -x;
 

libquadmath/math/j0q.c

@@ -807,6 +807,7 @@ static __float128 Y0_2D[NY0_2D + 1] = {
  /* 1.000000000000000000000000000000000000000E0 */
 };
 
+static const long double U0 = -7.3804295108687225274343927948483016310862e-02Q;
 
 /* Bessel function of the second kind, order zero.  */
 
@@ -829,6 +830,8 @@ y0q (__float128 x)
       return -HUGE_VALQ + x;
     }
   xx = fabsq (x);
+  if (xx <= 0x1p-57)
+    return U0 + TWOOPI * logq (x);
   if (xx <= 2.0Q)
     {
       /* 0 <= x <= 2 */

libquadmath/math/j1q.c

@@ -836,8 +836,10 @@ y1q (__float128 x)
       return -HUGE_VALQ + x;
     }
   xx = fabsq (x);
+  if (xx <= 0x1p-114)
+    return -TWOOPI / x;
   if (xx <= 2.0Q)
-    {
+   {
       /* 0 <= x <= 2 */
       z = xx * xx;
       p = xx * neval (z, Y0_2N, NY0_2N) / deval (z, Y0_2D, NY0_2D);

libquadmath/math/log1pq.c

@@ -136,6 +136,12 @@ log1pq (__float128 xm1)
       && (u.words32.w1 | u.words32.w2 | u.words32.w3) == 0)
     return xm1;
 
+  if ((hx & 0x7fffffff) < 0x3f8e0000)
+    {
+      if ((int) xm1 == 0)
+       return xm1;
+    }
+
   x = xm1 + 1.0Q;
 
   /* log1p(-1) = -inf */

libquadmath/math/powq.c

@@ -148,7 +148,7 @@ powq (__float128 x, __float128 y)
 {
   __float128 z, ax, z_h, z_l, p_h, p_l;
   __float128 y1, t1, t2, r, s, t, u, v, w;
-  __float128 s2, s_h, s_l, t_h, t_l;
+  __float128 s2, s_h, s_l, t_h, t_l, ay;
   int32_t i, j, k, yisint, n;
   uint32_t ix, iy;
   int32_t hx, hy;
@@ -281,6 +281,10 @@ powq (__float128 x, __float128 y)
 	return (hy > 0) ? huge * huge : tiny * tiny;
     }
 
+  ay = y > 0 ? y : -y;
+  if (ay < 0x1p-128)
+    y = y < 0 ? -0x1p-128 : 0x1p-128;
+
   n = 0;
   /* take care subnormal number */
   if (ix < 0x00010000)