Fix dbl-64 e_sqrt.c for non-default rounding modes (bug 16271).

9 files changed, 42 insertions, 4 deletions

diff --git a/sysdeps/i386/fpu/fegetround.c b/sysdeps/i386/fpu/fegetround.c
index d0170d3c86..cd96ae99d3 100644
--- a/sysdeps/i386/fpu/fegetround.c
+++ b/sysdeps/i386/fpu/fegetround.c
@@ -28,3 +28,4 @@ fegetround (void)
 
   return cw & 0xc00;
 }
+libm_hidden_def (fegetround)
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.  */
diff --git a/sysdeps/powerpc/fpu/fegetround.c b/sysdeps/powerpc/fpu/fegetround.c
index bcb6caab9d..078911f4a3 100644
--- a/sysdeps/powerpc/fpu/fegetround.c
+++ b/sysdeps/powerpc/fpu/fegetround.c
@@ -24,3 +24,4 @@ fegetround (void)
 {
   return __fegetround();
 }
+libm_hidden_def (fegetround)
diff --git a/sysdeps/powerpc/nofpu/fegetround.c b/sysdeps/powerpc/nofpu/fegetround.c
index 016602fac6..2c7bdbe5f6 100644
--- a/sysdeps/powerpc/nofpu/fegetround.c
+++ b/sysdeps/powerpc/nofpu/fegetround.c
@@ -26,3 +26,4 @@ fegetround (void)
 {
   return __sim_round_mode_thread;
 }
+libm_hidden_def (fegetround)
diff --git a/sysdeps/powerpc/powerpc32/e500/nofpu/fegetround.c b/sysdeps/powerpc/powerpc32/e500/nofpu/fegetround.c
index f69e9a5bdb..1e894e7523 100644
--- a/sysdeps/powerpc/powerpc32/e500/nofpu/fegetround.c
+++ b/sysdeps/powerpc/powerpc32/e500/nofpu/fegetround.c
@@ -27,3 +27,4 @@ fegetround (void)
   fpescr = fegetenv_register ();
   return fpescr & 3;
 }
+libm_hidden_def (fegetround)
diff --git a/sysdeps/s390/fpu/fegetround.c b/sysdeps/s390/fpu/fegetround.c
index 4843a56d26..94482f6318 100644
--- a/sysdeps/s390/fpu/fegetround.c
+++ b/sysdeps/s390/fpu/fegetround.c
@@ -29,3 +29,4 @@ fegetround (void)
 
   return cw & FPC_RM_MASK;
 }
+libm_hidden_def (fegetround)
diff --git a/sysdeps/sh/sh4/fpu/fegetround.c b/sysdeps/sh/sh4/fpu/fegetround.c
index be4833f017..0523321b2d 100644
--- a/sysdeps/sh/sh4/fpu/fegetround.c
+++ b/sysdeps/sh/sh4/fpu/fegetround.c
@@ -30,3 +30,4 @@ fegetround (void)
 
   return cw & 0x1;
 }
+libm_hidden_def (fegetround)
diff --git a/sysdeps/sparc/fpu/fegetround.c b/sysdeps/sparc/fpu/fegetround.c
index c4987e8b3e..c2d5f5af03 100644
--- a/sysdeps/sparc/fpu/fegetround.c
+++ b/sysdeps/sparc/fpu/fegetround.c
@@ -27,3 +27,4 @@ fegetround (void)
 
   return tmp & __FE_ROUND_MASK;
 }
+libm_hidden_def (fegetround)
diff --git a/sysdeps/x86_64/fpu/fegetround.c b/sysdeps/x86_64/fpu/fegetround.c
index 1a52b7ea67..c7cd046f39 100644
--- a/sysdeps/x86_64/fpu/fegetround.c
+++ b/sysdeps/x86_64/fpu/fegetround.c
@@ -30,3 +30,4 @@ fegetround (void)
 
   return cw & 0xc00;
 }
+libm_hidden_def (fegetround)