Merge branch 'master' into gcngcn

1 files changed, 66 insertions, 37 deletions

diff --git a/libquadmath/math/j1q.c b/libquadmath/math/j1q.c
index eb599c949a9..5eb705084e2 100644
--- a/libquadmath/math/j1q.c
+++ b/libquadmath/math/j1q.c
@@ -95,6 +95,7 @@
     License along with this library; if not, write to the Free Software
     Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307  USA */
 
+#include <errno.h>
 #include "quadmath-imp.h"
 
 /* 1 / sqrt(pi) */
@@ -687,13 +688,21 @@ j1q (__float128 x)
   if (! finiteq (x))
     {
       if (x != x)
-	return x;
+	return x + x;
       else
 	return 0.0Q;
     }
   if (x == 0.0Q)
     return x;
   xx = fabsq (x);
+  if (xx <= 0x1p-58Q)
+    {
+      __float128 ret = x * 0.5Q;
+      math_check_force_underflow (ret);
+      if (ret == 0)
+	errno = ERANGE;
+      return ret;
+    }
   if (xx <= 2.0Q)
     {
       /* 0 <= x <= 2 */
@@ -705,6 +714,32 @@ j1q (__float128 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.  */
+  sincosq (xx, &s, &c);
+  ss = -s - c;
+  cc = s - c;
+  if (xx <= FLT128_MAX / 2.0Q)
+    {
+      z = cosq (xx + xx);
+      if ((s * c) > 0)
+	cc = z / ss;
+      else
+	ss = z / cc;
+    }
+
+  if (xx > 0x1p256Q)
+    {
+      z = ONEOSQPI * cc / sqrtq (xx);
+      if (x < 0)
+	z = -z;
+      return z;
+    }
+
   xinv = 1.0Q / xx;
   z = xinv * xinv;
   if (xinv <= 0.25)
@@ -766,20 +801,6 @@ j1q (__float128 x)
   p = 1.0Q + z * p;
   q = z * q;
   q = q * xinv + 0.375Q * 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.  */
-  sincosq (xx, &s, &c);
-  ss = -s - c;
-  cc = s - c;
-  z = cosq (xx + xx);
-  if ((s * c) > 0)
-    cc = z / ss;
-  else
-    ss = z / cc;
   z = ONEOSQPI * (p * cc - q * ss) / sqrtq (xx);
   if (x < 0)
     z = -z;
@@ -823,24 +844,25 @@ y1q (__float128 x)
   __float128 xx, xinv, z, p, q, c, s, cc, ss;
 
   if (! finiteq (x))
-    {
-      if (x != x)
-	return x;
-      else
-	return 0.0Q;
-    }
+    return 1 / (x + x * x);
   if (x <= 0.0Q)
     {
       if (x < 0.0Q)
 	return (zero / (zero * x));
-      return -HUGE_VALQ + x;
+      return -1 / zero; /* -inf and divide by zero exception.  */
     }
   xx = fabsq (x);
   if (xx <= 0x1p-114)
-    return -TWOOPI / x;
+    {
+      z = -TWOOPI / x;
+      if (isinfq (z))
+	errno = ERANGE;
+      return z;
+    }
   if (xx <= 2.0Q)
    {
       /* 0 <= x <= 2 */
+      /* FIXME:  SET_RESTORE_ROUNDL (FE_TONEAREST);  */
       z = xx * xx;
       p = xx * neval (z, Y0_2N, NY0_2N) / deval (z, Y0_2D, NY0_2D);
       p = -TWOOPI / xx + p;
@@ -848,6 +870,27 @@ y1q (__float128 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.  */
+  sincosq (xx, &s, &c);
+  ss = -s - c;
+  cc = s - c;
+  if (xx <= FLT128_MAX / 2.0Q)
+    {
+      z = cosq (xx + xx);
+      if ((s * c) > 0)
+	cc = z / ss;
+      else
+	ss = z / cc;
+    }
+
+  if (xx > 0x1p256Q)
+    return ONEOSQPI * ss / sqrtq (xx);
+
   xinv = 1.0Q / xx;
   z = xinv * xinv;
   if (xinv <= 0.25)
@@ -909,20 +952,6 @@ y1q (__float128 x)
   p = 1.0Q + z * p;
   q = z * q;
   q = q * xinv + 0.375Q * 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.  */
-  sincosq (xx, &s, &c);
-  ss = -s - c;
-  cc = s - c;
-  z = cosq (xx + xx);
-  if ((s * c) > 0)
-    cc = z / ss;
-  else
-    ss = z / cc;
   z = ONEOSQPI * (p * ss + q * cc) / sqrtq (xx);
   return z;
 }