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

diff --git a/libquadmath/ChangeLog b/libquadmath/ChangeLog
index b97a45868b3b1d0523d2d29678d2787c3f0e8108..ffd5927b72e57c973d7341394fe15039508df57e 100644 (file)
--- a/libquadmath/ChangeLog
+++ b/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

diff --git a/libquadmath/math/atanq.c b/libquadmath/math/atanq.c
index cb38a340a205272e1b556210ef08f1a6a8c1dae5..8eccdc3317dab7bc5b389251e4e1c5e13c1af444 100644 (file)
--- a/libquadmath/math/atanq.c
+++ b/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;
 

diff --git a/libquadmath/math/j0q.c b/libquadmath/math/j0q.c
index 3ef9201a8a42407aeafdd360c925a42010827673..8c6f811125e702c230576079ad2725e25d2c0cba 100644 (file)
--- a/libquadmath/math/j0q.c
+++ b/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 */

diff --git a/libquadmath/math/j1q.c b/libquadmath/math/j1q.c
index 8a18e2779b2b79c77ed1098523e66002a88156a1..eb599c949a90e3aa97c13bb8b28c921084a689db 100644 (file)
--- a/libquadmath/math/j1q.c
+++ b/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);

diff --git a/libquadmath/math/log1pq.c b/libquadmath/math/log1pq.c
index 2de13fe2adcc94067db062505c0d7ce29492496d..d8bff405dff32c59a901ef3fd71f8ba1a44d252b 100644 (file)
--- a/libquadmath/math/log1pq.c
+++ b/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 */

diff --git a/libquadmath/math/powq.c b/libquadmath/math/powq.c
index 12b87d536d7250660f8b7b12f58094b9ed1918fb..dd44b7c175a011f1a1e661105d76a2ce5df600b6 100644 (file)
--- a/libquadmath/math/powq.c
+++ b/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)