This is the mail archive of the libc-alpha@sourceware.org mailing list for the glibc project.


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]
Other format: [Raw text]

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


This patch fixes bug 16271, incorrect rounding by dbl-64/e_sqrt.c in
non-default rounding modes.  The calculation is carried out in
round-to-nearest (since that's what any error analysis would have been
based on) and then the round-to-nearest result is adjusted for other
rounding modes based on the result of a division that is carried out
anyway to get "inexact" exceptions right.

In principle the call to fegetround (and consequent libm_hidden_proto
/ libm_hidden_def changes) shouldn't be needed as the rounding mode
information could be extracted from the fenv_t from
libc_feholdexcept_setround.  However, it looks like no architecture
with optimized versions of libc_fe* actually uses this file, so it
doesn't seem worthwhile implementing any optimization infrastructure
at this point for architectures to describe how to extract a rounding
mode without a separate fesetround call; if someone adds libc_fe*
optimizations for such an architecture, they can also optimize this
case if desired.

This patch does not depend on any of my other recent patches relating
to sqrt and sqrt tests, except that those patches provide the
testcases that show the bug being fixed here.

Tested for ARM.

2013-11-27  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.

ports/ChangeLog.aarch64:
2013-11-27  Joseph Myers  <joseph@codesourcery.com>

	* sysdeps/aarch64/fpu/fegetround.c (fegetround): Use
	libm_hidden_def.

ports/ChangeLog.alpha:
2013-11-27  Joseph Myers  <joseph@codesourcery.com>

	* sysdeps/alpha/fpu/fegetround.c (fegetround): Use
	libm_hidden_def.

ports/ChangeLog.am33:
2013-11-27  Joseph Myers  <joseph@codesourcery.com>

	* sysdeps/am33/fpu/fegetround.c (fegetround): Use libm_hidden_def.

ports/ChangeLog.arm:
2013-11-27  Joseph Myers  <joseph@codesourcery.com>

	* sysdeps/arm/fegetround.c (fegetround): Use libm_hidden_def.

ports/ChangeLog.hppa:
2013-11-27  Joseph Myers  <joseph@codesourcery.com>

	* sysdeps/hppa/fpu/fegetround.c (fegetround): Use libm_hidden_def.

ports/ChangeLog.ia64:
2013-11-27  Joseph Myers  <joseph@codesourcery.com>

	* sysdeps/ia64/fpu/fegetround.c (fegetround): Use libm_hidden_def.

ports/ChangeLog.m68k:
2013-11-27  Joseph Myers  <joseph@codesourcery.com>

	* sysdeps/m68k/fpu/fegetround.c (fegetround): Use libm_hidden_def.

ports/ChangeLog.microblaze:
2013-11-27  Joseph Myers  <joseph@codesourcery.com>

	* sysdeps/microblaze/fegetround.c (fegetround): Use
	libm_hidden_def.

ports/ChangeLog.mips:
2013-11-27  Joseph Myers  <joseph@codesourcery.com>

	* sysdeps/mips/fpu/fegetround.c (fegetround): Use libm_hidden_def.

diff --git a/include/fenv.h b/include/fenv.h
index 925d4b5..bd2c99d 100644
--- a/include/fenv.h
+++ b/include/fenv.h
@@ -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)
diff --git a/math/fegetround.c b/math/fegetround.c
index 24bbd16..140e698 100644
--- a/math/fegetround.c
+++ b/math/fegetround.c
@@ -28,4 +28,5 @@ fegetround (void)
   return 0;
 #endif
 }
+libm_hidden_def (fegetround)
 stub_warning (fegetround)
diff --git a/ports/sysdeps/aarch64/fpu/fegetround.c b/ports/sysdeps/aarch64/fpu/fegetround.c
index 3b5b306..370caa1 100644
--- a/ports/sysdeps/aarch64/fpu/fegetround.c
+++ b/ports/sysdeps/aarch64/fpu/fegetround.c
@@ -26,3 +26,4 @@ fegetround (void)
   _FPU_GETCW (fpcr);
   return fpcr & FE_TOWARDZERO;
 }
+libm_hidden_def (fegetround)
diff --git a/ports/sysdeps/alpha/fpu/fegetround.c b/ports/sysdeps/alpha/fpu/fegetround.c
index aba657a..03a55ee 100644
--- a/ports/sysdeps/alpha/fpu/fegetround.c
+++ b/ports/sysdeps/alpha/fpu/fegetround.c
@@ -28,3 +28,4 @@ fegetround (void)
 
   return (fpcr >> FPCR_ROUND_SHIFT) & 3;
 }
+libm_hidden_def (fegetround)
diff --git a/ports/sysdeps/am33/fpu/fegetround.c b/ports/sysdeps/am33/fpu/fegetround.c
index b309c92..49cae00 100644
--- a/ports/sysdeps/am33/fpu/fegetround.c
+++ b/ports/sysdeps/am33/fpu/fegetround.c
@@ -32,3 +32,4 @@ fegetround (void)
 
   return (cw & ROUND_MASK);
 }
+libm_hidden_def (fegetround)
diff --git a/ports/sysdeps/arm/fegetround.c b/ports/sysdeps/arm/fegetround.c
index 78a3795..149a989 100644
--- a/ports/sysdeps/arm/fegetround.c
+++ b/ports/sysdeps/arm/fegetround.c
@@ -37,3 +37,4 @@ fegetround (void)
   /* The current soft-float implementation only handles TONEAREST.  */
   return FE_TONEAREST;
 }
+libm_hidden_def (fegetround)
diff --git a/ports/sysdeps/hppa/fpu/fegetround.c b/ports/sysdeps/hppa/fpu/fegetround.c
index 67dd7c4..3815fbd 100644
--- a/ports/sysdeps/hppa/fpu/fegetround.c
+++ b/ports/sysdeps/hppa/fpu/fegetround.c
@@ -24,3 +24,4 @@ fegetround (void)
 {
   return get_rounding_mode ();
 }
+libm_hidden_def (fegetround)
diff --git a/ports/sysdeps/ia64/fpu/fegetround.c b/ports/sysdeps/ia64/fpu/fegetround.c
index 5c9b343..f6dfea7 100644
--- a/ports/sysdeps/ia64/fpu/fegetround.c
+++ b/ports/sysdeps/ia64/fpu/fegetround.c
@@ -24,3 +24,4 @@ fegetround (void)
 {
   return get_rounding_mode ();
 }
+libm_hidden_def (fegetround)
diff --git a/ports/sysdeps/m68k/fpu/fegetround.c b/ports/sysdeps/m68k/fpu/fegetround.c
index f1227fe..54fa7df 100644
--- a/ports/sysdeps/m68k/fpu/fegetround.c
+++ b/ports/sysdeps/m68k/fpu/fegetround.c
@@ -28,3 +28,4 @@ fegetround (void)
 
   return fpcr & FE_UPWARD;
 }
+libm_hidden_def (fegetround)
diff --git a/ports/sysdeps/microblaze/fegetround.c b/ports/sysdeps/microblaze/fegetround.c
index 4f47dd1..b1039e8 100644
--- a/ports/sysdeps/microblaze/fegetround.c
+++ b/ports/sysdeps/microblaze/fegetround.c
@@ -22,3 +22,4 @@ fegetround (void)
 {
   return FE_TONEAREST;
 }
+libm_hidden_def (fegetround)
diff --git a/ports/sysdeps/mips/fpu/fegetround.c b/ports/sysdeps/mips/fpu/fegetround.c
index 17cd3e9..011d27f 100644
--- a/ports/sysdeps/mips/fpu/fegetround.c
+++ b/ports/sysdeps/mips/fpu/fegetround.c
@@ -30,3 +30,4 @@ fegetround (void)
 
   return cw & _FPU_RC_MASK;
 }
+libm_hidden_def (fegetround)
diff --git a/sysdeps/i386/fpu/fegetround.c b/sysdeps/i386/fpu/fegetround.c
index d0170d3..cd96ae9 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 854ae38..88809da 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 bcb6caa..078911f 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 016602f..2c7bdbe 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 f69e9a5..1e894e7 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 4843a56..94482f6 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 be4833f..0523321 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 c4987e8..c2d5f5a 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 1a52b7e..c7cd046 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)

-- 
Joseph S. Myers
joseph@codesourcery.com


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]