about summary refs log tree commit diff
path: root/sysdeps/ieee754
diff options
context:
space:
mode:
authorJoseph Myers <joseph@codesourcery.com>2013-11-28 16:50:38 +0000
committerJoseph Myers <joseph@codesourcery.com>2013-11-28 16:50:38 +0000
commit3c1c46a64ad1037d616ec39514c4e55133997c9f (patch)
tree717fc72e73a01fb74a4cdc97e6332086c2c3ee38 /sysdeps/ieee754
parent5a4c6d53f50b264d60cf6453576ca2810c7890b7 (diff)
downloadglibc-3c1c46a64ad1037d616ec39514c4e55133997c9f.tar.gz
glibc-3c1c46a64ad1037d616ec39514c4e55133997c9f.tar.xz
glibc-3c1c46a64ad1037d616ec39514c4e55133997c9f.zip
Fix dbl-64 e_sqrt.c for non-default rounding modes (bug 16271).
Diffstat (limited to 'sysdeps/ieee754')
-rw-r--r--sysdeps/ieee754/dbl-64/e_sqrt.c38
1 files changed, 34 insertions, 4 deletions
diff --git a/sysdeps/ieee754/dbl-64/e_sqrt.c b/sysdeps/ieee754/dbl-64/e_sqrt.c
index 854ae38c41..88809daa76 100644
--- a/sysdeps/ieee754/dbl-64/e_sqrt.c
+++ b/sysdeps/ieee754/dbl-64/e_sqrt.c
@@ -66,8 +66,9 @@ __ieee754_sqrt (double x)
   /*----------------- 2^-1022  <= | x |< 2^1024  -----------------*/
   if (k > 0x000fffff && k < 0x7ff00000)
     {
+      int rm = fegetround ();
       fenv_t env;
-      libc_feholdexcept (&env);
+      libc_feholdexcept_setround (&env, FE_TONEAREST);
       double ret;
       y = 1.0 - t * (t * s);
       t = t * (rt0 + y * (rt1 + y * (rt2 + y * rt3)));
@@ -82,15 +83,44 @@ __ieee754_sqrt (double x)
 	{
 	  res1 = res + 1.5 * ((y - res) + del);
 	  EMULV (res, res1, z, zz, p, hx, tx, hy, ty); /* (z+zz)=res*res1 */
-	  ret = ((((z - s) + zz) < 0) ? max (res, res1) :
-					min (res, res1)) * c.x;
+	  res = ((((z - s) + zz) < 0) ? max (res, res1) :
+					min (res, res1));
+	  ret = res * c.x;
 	}
       math_force_eval (ret);
       libc_fesetenv (&env);
-      if (x / ret != ret)
+      double dret = x / ret;
+      if (dret != ret)
 	{
 	  double force_inexact = 1.0 / 3.0;
 	  math_force_eval (force_inexact);
+	  /* The square root is inexact, ret is the round-to-nearest
+	     value which may need adjusting for other rounding
+	     modes.  */
+	  switch (rm)
+	    {
+#ifdef FE_UPWARD
+	    case FE_UPWARD:
+	      if (dret > ret)
+		ret = (res + 0x1p-1022) * c.x;
+	      break;
+#endif
+
+#ifdef FE_DOWNWARD
+	    case FE_DOWNWARD:
+#endif
+#ifdef FE_TOWARDZERO
+	    case FE_TOWARDZERO:
+#endif
+#if defined FE_DOWNWARD || defined FE_TOWARDZERO
+	      if (dret < ret)
+		ret = (res - 0x1p-1022) * c.x;
+	      break;
+#endif
+
+	    default:
+	      break;
+	    }
 	}
       /* Otherwise (x / ret == ret), either the square root was exact or
          the division was inexact.  */