aboutsummaryrefslogtreecommitdiff
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
glibc-3c1c46a64ad1037d616ec39514c4e55133997c9f.tar.gz
glibc-3c1c46a64ad1037d616ec39514c4e55133997c9f.tar.bz2
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. */