Fix dbl-64 e_sqrt.c for non-default rounding modes (bug 16271).
authorJoseph Myers <joseph@codesourcery.com>
Thu, 28 Nov 2013 16:50:38 +0000 (16:50 +0000)
committerJoseph Myers <joseph@codesourcery.com>
Thu, 28 Nov 2013 16:50:38 +0000 (16:50 +0000)
31 files changed:
ChangeLog
NEWS
include/fenv.h
math/fegetround.c
ports/ChangeLog.aarch64
ports/ChangeLog.alpha
ports/ChangeLog.am33
ports/ChangeLog.arm
ports/ChangeLog.hppa
ports/ChangeLog.ia64
ports/ChangeLog.m68k
ports/ChangeLog.microblaze
ports/ChangeLog.mips
ports/sysdeps/aarch64/fpu/fegetround.c
ports/sysdeps/alpha/fpu/fegetround.c
ports/sysdeps/am33/fpu/fegetround.c
ports/sysdeps/arm/fegetround.c
ports/sysdeps/hppa/fpu/fegetround.c
ports/sysdeps/ia64/fpu/fegetround.c
ports/sysdeps/m68k/fpu/fegetround.c
ports/sysdeps/microblaze/fegetround.c
ports/sysdeps/mips/fpu/fegetround.c
sysdeps/i386/fpu/fegetround.c
sysdeps/ieee754/dbl-64/e_sqrt.c
sysdeps/powerpc/fpu/fegetround.c
sysdeps/powerpc/nofpu/fegetround.c
sysdeps/powerpc/powerpc32/e500/nofpu/fegetround.c
sysdeps/s390/fpu/fegetround.c
sysdeps/sh/sh4/fpu/fegetround.c
sysdeps/sparc/fpu/fegetround.c
sysdeps/x86_64/fpu/fegetround.c

index 7e46a66..16aa48a 100644 (file)
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,20 @@
+2013-11-28  Joseph Myers  <joseph@codesourcery.com>
+
+       [BZ #16271]
+       * sysdeps/ieee754/dbl-64/e_sqrt.c (__ieee754_sqrt): Set
+       round-to-nearest then adjust result for other rounding modes.
+       * include/fenv.h (fegetround): Use libm_hidden_proto.
+       * math/fegetround.c (fegetround): Use libm_hidden_def.
+       * sysdeps/i386/fpu/fegetround.c (fegetround): Likewise.
+       * sysdeps/powerpc/fpu/fegetround.c (fegetround): Likewise.
+       * sysdeps/powerpc/nofpu/fegetround.c (fegetround): Likewise.
+       * sysdeps/powerpc/powerpc32/e500/nofpu/fegetround.c (fegetround):
+       Likewise.
+       * sysdeps/s390/fpu/fegetround.c (fegetround): Likewise.
+       * sysdeps/sh/sh4/fpu/fegetround.c (fegetround): Likewise.
+       * sysdeps/sparc/fpu/fegetround.c (fegetround): Likewise.
+       * sysdeps/x86_64/fpu/fegetround.c (fegetround): Likewise.
+
 2013-11-28  Siddhesh Poyarekar  <siddhesh@redhat.com>
 
        [BZ #16077]
diff --git a/NEWS b/NEWS
index b45fa53..4a8e3c0 100644 (file)
--- a/NEWS
+++ b/NEWS
@@ -19,7 +19,7 @@ Version 2.19
   15897, 15905, 15909, 15917, 15919, 15921, 15923, 15939, 15948, 15963,
   15966, 15985, 15988, 15997, 16032, 16034, 16036, 16037, 16041, 16055,
   16071, 16072, 16074, 16077, 16078, 16103, 16112, 16143, 16144, 16146,
-  16150, 16151, 16153, 16167, 16172, 16245.
+  16150, 16151, 16153, 16167, 16172, 16245, 16271.
 
 * The public headers no longer use __unused nor __block.  This change is to
   support compiling programs that are derived from BSD sources and use
index 925d4b5..bd2c99d 100644 (file)
@@ -16,6 +16,7 @@ extern int __feupdateenv (const fenv_t *__envp);
 
 libm_hidden_proto (feraiseexcept)
 libm_hidden_proto (fegetenv)
+libm_hidden_proto (fegetround)
 libm_hidden_proto (fesetenv)
 libm_hidden_proto (fesetround)
 libm_hidden_proto (feholdexcept)
index 24bbd16..140e698 100644 (file)
@@ -28,4 +28,5 @@ fegetround (void)
   return 0;
 #endif
 }
+libm_hidden_def (fegetround)
 stub_warning (fegetround)
index c6331a4..4cbc11f 100644 (file)
@@ -1,3 +1,8 @@
+2013-11-28  Joseph Myers  <joseph@codesourcery.com>
+
+       * sysdeps/aarch64/fpu/fegetround.c (fegetround): Use
+       libm_hidden_def.
+
 2013-11-26  Will Newton  <will.newton@linaro.org>
 
        * sysdeps/aarch64/dl-irel.h: Include ldsodefs.h.
index 05c6ba4..95b575e 100644 (file)
@@ -1,3 +1,8 @@
+2013-11-28  Joseph Myers  <joseph@codesourcery.com>
+
+       * sysdeps/alpha/fpu/fegetround.c (fegetround): Use
+       libm_hidden_def.
+
 2013-11-26  Ondřej Bílka  <neleai@seznam.cz>
        * sysdeps/unix/sysv/linux/alpha/bits/ipc.h: Use __glibc_reserved instead __unused.
        * sysdeps/unix/sysv/linux/alpha/bits/msq.h: Likewise.
index 15b7a1b..317dd8e 100644 (file)
@@ -1,3 +1,7 @@
+2013-11-28  Joseph Myers  <joseph@codesourcery.com>
+
+       * sysdeps/am33/fpu/fegetround.c (fegetround): Use libm_hidden_def.
+
 2013-10-30  Mike Frysinger  <vapier@gentoo.org>
 
        * sysdeps/unix/sysv/linux/am33/configure.in: Moved to ...
index 4124e47..4a4d319 100644 (file)
@@ -1,3 +1,7 @@
+2013-11-28  Joseph Myers  <joseph@codesourcery.com>
+
+       * sysdeps/arm/fegetround.c (fegetround): Use libm_hidden_def.
+
 2013-11-26  Ondřej Bílka  <neleai@seznam.cz>
        * sysdeps/unix/sysv/linux/arm/bits/shm.h: Use __glibc_reserved instead __unused.
 
index 6a28bfe..98d36a3 100644 (file)
@@ -1,3 +1,7 @@
+2013-11-28  Joseph Myers  <joseph@codesourcery.com>
+
+       * sysdeps/hppa/fpu/fegetround.c (fegetround): Use libm_hidden_def.
+
 2013-11-26  Ondřej Bílka  <neleai@seznam.cz>
        * sysdeps/unix/sysv/linux/hppa/bits/ipc.h: Use __glibc_reserved instead __unused.
        * sysdeps/unix/sysv/linux/hppa/bits/msq.h: Likewise.
index d2825f1..1aded58 100644 (file)
@@ -1,3 +1,7 @@
+2013-11-28  Joseph Myers  <joseph@codesourcery.com>
+
+       * sysdeps/ia64/fpu/fegetround.c (fegetround): Use libm_hidden_def.
+
 2013-11-26  Ondřej Bílka  <neleai@seznam.cz>
        * sysdeps/unix/sysv/linux/ia64/bits/ipc.h: Use __glibc_reserved instead __unused.
        * sysdeps/unix/sysv/linux/ia64/bits/msq.h: Likewise.
index f78f384..92ff973 100644 (file)
@@ -1,3 +1,7 @@
+2013-11-28  Joseph Myers  <joseph@codesourcery.com>
+
+       * sysdeps/m68k/fpu/fegetround.c (fegetround): Use libm_hidden_def.
+
 2013-11-26  Ondřej Bílka  <neleai@seznam.cz>
        * sysdeps/unix/sysv/linux/m68k/bits/stat.h: Use __glibc_reserved instead __unused.
 
index 0762dd3..7d47e04 100644 (file)
@@ -1,3 +1,8 @@
+2013-11-28  Joseph Myers  <joseph@codesourcery.com>
+
+       * sysdeps/microblaze/fegetround.c (fegetround): Use
+       libm_hidden_def.
+
 2013-11-26  Ondřej Bílka  <neleai@seznam.cz>
        * sysdeps/unix/sysv/linux/microblaze/bits/stat.h: Use __glibc_reserved instead __unused.
        * sysdeps/unix/sysv/linux/microblaze/kernel_stat.h: Likewise.
index 8456799..21308e7 100644 (file)
@@ -1,3 +1,7 @@
+2013-11-28  Joseph Myers  <joseph@codesourcery.com>
+
+       * sysdeps/mips/fpu/fegetround.c (fegetround): Use libm_hidden_def.
+
 2013-11-27  Aurelien Jarno <aurelien@aurel32.net>
 
        * sysdeps/unix/sysv/linux/mips/bits/resource.h (RLIM64_INFINITY): Fix
index 3b5b306..370caa1 100644 (file)
@@ -26,3 +26,4 @@ fegetround (void)
   _FPU_GETCW (fpcr);
   return fpcr & FE_TOWARDZERO;
 }
+libm_hidden_def (fegetround)
index aba657a..03a55ee 100644 (file)
@@ -28,3 +28,4 @@ fegetround (void)
 
   return (fpcr >> FPCR_ROUND_SHIFT) & 3;
 }
+libm_hidden_def (fegetround)
index b309c92..49cae00 100644 (file)
@@ -32,3 +32,4 @@ fegetround (void)
 
   return (cw & ROUND_MASK);
 }
+libm_hidden_def (fegetround)
index 78a3795..149a989 100644 (file)
@@ -37,3 +37,4 @@ fegetround (void)
   /* The current soft-float implementation only handles TONEAREST.  */
   return FE_TONEAREST;
 }
+libm_hidden_def (fegetround)
index 67dd7c4..3815fbd 100644 (file)
@@ -24,3 +24,4 @@ fegetround (void)
 {
   return get_rounding_mode ();
 }
+libm_hidden_def (fegetround)
index 5c9b343..f6dfea7 100644 (file)
@@ -24,3 +24,4 @@ fegetround (void)
 {
   return get_rounding_mode ();
 }
+libm_hidden_def (fegetround)
index f1227fe..54fa7df 100644 (file)
@@ -28,3 +28,4 @@ fegetround (void)
 
   return fpcr & FE_UPWARD;
 }
+libm_hidden_def (fegetround)
index 4f47dd1..b1039e8 100644 (file)
@@ -22,3 +22,4 @@ fegetround (void)
 {
   return FE_TONEAREST;
 }
+libm_hidden_def (fegetround)
index 17cd3e9..011d27f 100644 (file)
@@ -30,3 +30,4 @@ fegetround (void)
 
   return cw & _FPU_RC_MASK;
 }
+libm_hidden_def (fegetround)
index d0170d3..cd96ae9 100644 (file)
@@ -28,3 +28,4 @@ fegetround (void)
 
   return cw & 0xc00;
 }
+libm_hidden_def (fegetround)
index 854ae38..88809da 100644 (file)
@@ -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.  */
index bcb6caa..078911f 100644 (file)
@@ -24,3 +24,4 @@ fegetround (void)
 {
   return __fegetround();
 }
+libm_hidden_def (fegetround)
index 016602f..2c7bdbe 100644 (file)
@@ -26,3 +26,4 @@ fegetround (void)
 {
   return __sim_round_mode_thread;
 }
+libm_hidden_def (fegetround)
index f69e9a5..1e894e7 100644 (file)
@@ -27,3 +27,4 @@ fegetround (void)
   fpescr = fegetenv_register ();
   return fpescr & 3;
 }
+libm_hidden_def (fegetround)
index 4843a56..94482f6 100644 (file)
@@ -29,3 +29,4 @@ fegetround (void)
 
   return cw & FPC_RM_MASK;
 }
+libm_hidden_def (fegetround)
index be4833f..0523321 100644 (file)
@@ -30,3 +30,4 @@ fegetround (void)
 
   return cw & 0x1;
 }
+libm_hidden_def (fegetround)
index c4987e8..c2d5f5a 100644 (file)
@@ -27,3 +27,4 @@ fegetround (void)
 
   return tmp & __FE_ROUND_MASK;
 }
+libm_hidden_def (fegetround)
index 1a52b7e..c7cd046 100644 (file)
@@ -30,3 +30,4 @@ fegetround (void)
 
   return cw & 0xc00;
 }
+libm_hidden_def (fegetround)