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

diff --git a/ChangeLog b/ChangeLog
index bf4e137..af0d243 100644 (file)
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,5 +1,10 @@
 2013-06-15  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.
+
        * math/libm-test.inc (fma_test_data): Add another test.
 
 2013-06-15  Siddhesh Poyarekar  <siddhesh@redhat.com>

diff --git a/NEWS b/NEWS
index ea9c6a0..63f2892 100644 (file)
--- a/NEWS
+++ b/NEWS
@@ -20,7 +20,7 @@ Version 2.18
   15366, 15380, 15381, 15394, 15395, 15405, 15406, 15409, 15416, 15418,
   15419, 15423, 15424, 15426, 15429, 15431, 15432, 15441, 15442, 15448,
   15465, 15480, 15485, 15488, 15490, 15492, 15493, 15497, 15506, 15529,
-  15536, 15553, 15577, 15583, 15618, 15627.
+  15536, 15553, 15577, 15583, 15618, 15627, 15631.
 
 * CVE-2013-0242 Buffer overrun in regexp matcher has been fixed (Bugzilla
   #15078).

diff --git a/sysdeps/ieee754/dbl-64/e_sqrt.c b/sysdeps/ieee754/dbl-64/e_sqrt.c
index f90ea0c..54610ee 100644 (file)
--- 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)