Fix spurious "inexact" exceptions from dbl-64 sqrt (bug 15631)

["Joseph S. Myers" <joseph at codesourcery dot com>]

Bug 15631 is spurious "inexact" exceptions from the dbl-64 sqrt
implementation, which show up as testsuite failures for architectures
not using hardware square root instructions.  This patch fixes these
by saving and restoring exception state, which then requires an
explicit check for whether the computed result is the exact square
root so "inexact" can be raised when required.  Tested MIPS64 (n32
ABI).  (MIPS is in fact an architecture that *should* be using
hardware square root instructions, but as a result of sysdeps
directory ordering it isn't at present; I've filed bug 15632 for
that.)

(I suspect there may be directed rounding issues in this code - that
is, it may not always give correctly rounded results outside
round-to-nearest mode - but haven't searched for them yet.)

2013-06-14  Joseph Myers  <joseph@codesourcery.com>

	[BZ #15631]
	* sysdeps/ieee754/dbl-64/e_sqrt.c (__ieee754_sqrt): Save and
	restore exception state around main square root computation, then
	check for inexactness explicitly.

diff --git a/sysdeps/ieee754/dbl-64/e_sqrt.c b/sysdeps/ieee754/dbl-64/e_sqrt.c
index f90ea0c..54610ee 100644
--- a/sysdeps/ieee754/dbl-64/e_sqrt.c
+++ b/sysdeps/ieee754/dbl-64/e_sqrt.c
@@ -63,6 +63,9 @@ double __ieee754_sqrt(double x) {
   s=a.x;
   /*----------------- 2^-1022  <= | x |< 2^1024  -----------------*/
   if (k>0x000fffff && k<0x7ff00000) {
+    fenv_t env;
+    libc_feholdexcept (&env);
+    double ret;
     y=1.0-t*(t*s);
     t=t*(rt0+y*(rt1+y*(rt2+y*rt3)));
     c.i[HIGH_HALF]=0x20000000+((k&0x7fe00000)>>1);
@@ -70,12 +73,22 @@ double __ieee754_sqrt(double x) {
     hy=(y+big)-big;
     del=0.5*t*((s-hy*hy)-(y-hy)*(y+hy));
     res=y+del;
-    if (res == (res+1.002*((y-res)+del))) return res*c.x;
+    if (res == (res+1.002*((y-res)+del))) ret = res*c.x;
     else {
       res1=res+1.5*((y-res)+del);
       EMULV(res,res1,z,zz,p,hx,tx,hy,ty);  /* (z+zz)=res*res1 */
-      return ((((z-s)+zz)<0)?max(res,res1):min(res,res1))*c.x;
+      ret = ((((z-s)+zz)<0)?max(res,res1):min(res,res1))*c.x;
     }
+    math_force_eval (ret);
+    libc_fesetenv (&env);
+    if (x / ret != ret)
+      {
+	double force_inexact = 1.0 / 3.0;
+	math_force_eval (force_inexact);
+      }
+    /* Otherwise (x / ret == ret), either the square root was exact or
+       the division was inexact.  */
+    return ret;
   }
   else {
     if ((k & 0x7ff00000) == 0x7ff00000)

-- 
Joseph S. Myers
joseph@codesourcery.com