Fix clog / clog10 sign of zero result in round-downward mode (bug 16789).

This patch fixes bug 16789, incorrect sign of (real part) zero result
from clog and clog10 in round-downward mode, arising from that real
part being computed as 0 - 0.  To ensure that an underflow exception
occurred, the code used an underflowing value (the next term in the
series for log1p) in arithmetic computing the real part of the result,
yielding the problematic 0 - 0 computation in some cases even when the
mathematical result would be small but positive.  The patch changes
this code to use the math_force_eval approach to ensuring that an
underflowing computation actually occurs.  Tests of clog and clog10
are enabled in all rounding modes.

Tested x86_64 and x86 and ulps updated accordingly.

	[BZ #16789]
	* math/s_clog.c (__clog): Use math_force_eval to ensure underflow
	instead of using underflowing value in computing result.
	* math/s_clog10.c (__clog10): Likewise.
	* math/s_clog10f.c (__clog10f): Likewise.
	* math/s_clog10l.c (__clog10l): Likewise.
	* math/s_clogf.c (__clogf): Likewise.
	* math/s_clogl.c (__clogl): Likewise.
	* math/libm-test.inc (clog_test): Use ALL_RM_TEST.
	(clog10_test): Likewise.
	* sysdeps/i386/fpu/libm-test-ulps: Update.
	* sysdeps/x86_64/fpu/libm-test-ulps: Likewise.

ChangeLog

@@ -1,3 +1,18 @@
+2014-04-02  Joseph Myers  <joseph@codesourcery.com>
+
+	[BZ #16789]
+	* math/s_clog.c (__clog): Use math_force_eval to ensure underflow
+	instead of using underflowing value in computing result.
+	* math/s_clog10.c (__clog10): Likewise.
+	* math/s_clog10f.c (__clog10f): Likewise.
+	* math/s_clog10l.c (__clog10l): Likewise.
+	* math/s_clogf.c (__clogf): Likewise.
+	* math/s_clogl.c (__clogl): Likewise.
+	* math/libm-test.inc (clog_test): Use ALL_RM_TEST.
+	(clog10_test): Likewise.
+	* sysdeps/i386/fpu/libm-test-ulps: Update.
+	* sysdeps/x86_64/fpu/libm-test-ulps: Likewise.
+
 2014-04-02  Alan Modra  <amodra@gmail.com>
 
 	[BZ #16739]

NEWS

@@ -13,7 +13,8 @@ Version 2.20
   16357, 16362, 16447, 16532, 16545, 16574, 16599, 16600, 16609, 16610,
   16611, 16613, 16623, 16632, 16634, 16639, 16642, 16648, 16649, 16670,
   16674, 16677, 16680, 16683, 16689, 16695, 16701, 16706, 16707, 16712,
-  16713, 16714, 16731, 16739, 16743, 16758, 16759, 16760, 16770, 16786.
+  16713, 16714, 16731, 16739, 16743, 16758, 16759, 16760, 16770, 16786,
+  16789.
 
 * Running the testsuite no longer terminates as soon as a test fails.
   Instead, a file tests.sum (xtests.sum from "make xcheck") is generated,

math/libm-test.inc

@@ -6122,9 +6122,7 @@ static const struct test_c_c_data clog_test_data[] =
 static void
 clog_test (void)
 {
-  START (clog, 0);
+  ALL_RM_TEST (clog, 0, clog_test_data, RUN_TEST_LOOP_c_c, END_COMPLEX);
-  RUN_TEST_LOOP_c_c (clog, clog_test_data, );
-  END_COMPLEX;
 }
 
 
@@ -6184,9 +6182,7 @@ static const struct test_c_c_data clog10_test_data[] =
 static void
 clog10_test (void)
 {
-  START (clog10, 0);
+  ALL_RM_TEST (clog10, 0, clog10_test_data, RUN_TEST_LOOP_c_c, END_COMPLEX);
-  RUN_TEST_LOOP_c_c (clog10, clog10_test_data, );
-  END_COMPLEX;
 }
 
 

math/s_clog.c

@@ -68,12 +68,9 @@ __clog (__complex__ double x)
 	  double absy2 = absy * absy;
 	  if (absy2 <= DBL_MIN * 2.0)
 	    {
-#if __FLT_EVAL_METHOD__ == 0
+	      double force_underflow = absy2 * absy2;
-	      __real__ result = absy2 / 2.0 - absy2 * absy2 / 4.0;
+	      __real__ result = absy2 / 2.0;
-#else
+	      math_force_eval (force_underflow);
-	      volatile double force_underflow = absy2 * absy2 / 4.0;
-	      __real__ result = absy2 / 2.0 - force_underflow;
-#endif
 	    }
 	  else
 	    __real__ result = __log1p (absy2) / 2.0;

math/s_clog10.c

@@ -74,12 +74,9 @@ __clog10 (__complex__ double x)
 	  double absy2 = absy * absy;
 	  if (absy2 <= DBL_MIN * 2.0 * M_LN10)
 	    {
-#if __FLT_EVAL_METHOD__ == 0
+	      double force_underflow = absy2 * absy2;
-	      __real__ result = (absy2 / 2.0 - absy2 * absy2 / 4.0) * M_LOG10E;
+	      __real__ result = absy2 * (M_LOG10E / 2.0);
-#else
+	      math_force_eval (force_underflow);
-	      volatile double force_underflow = absy2 * absy2 / 4.0;
-	      __real__ result = (absy2 / 2.0 - force_underflow) * M_LOG10E;
-#endif
 	    }
 	  else
 	    __real__ result = __log1p (absy2) * (M_LOG10E / 2.0);

math/s_clog10f.c

@@ -74,14 +74,9 @@ __clog10f (__complex__ float x)
 	  float absy2 = absy * absy;
 	  if (absy2 <= FLT_MIN * 2.0f * (float) M_LN10)
 	    {
-#if __FLT_EVAL_METHOD__ == 0
+	      float force_underflow = absy2 * absy2;
-	      __real__ result
+	      __real__ result = absy2 * ((float) M_LOG10E / 2.0f);
-		= (absy2 / 2.0f - absy2 * absy2 / 4.0f) * (float) M_LOG10E;
+	      math_force_eval (force_underflow);
-#else
-	      volatile float force_underflow = absy2 * absy2 / 4.0f;
-	      __real__ result
-		= (absy2 / 2.0f - force_underflow) * (float) M_LOG10E;
-#endif
 	    }
 	  else
 	    __real__ result = __log1pf (absy2) * ((float) M_LOG10E / 2.0f);

math/s_clog10l.c

@@ -80,8 +80,11 @@ __clog10l (__complex__ long double x)
 	{
 	  long double absy2 = absy * absy;
 	  if (absy2 <= LDBL_MIN * 2.0L * M_LN10l)
-	    __real__ result
+	    {
-	      = (absy2 / 2.0L - absy2 * absy2 / 4.0L) * M_LOG10El;
+	      long double force_underflow = absy2 * absy2;
+	      __real__ result = absy2 * (M_LOG10El / 2.0);
+	      math_force_eval (force_underflow);
+	    }
 	  else
 	    __real__ result = __log1pl (absy2) * (M_LOG10El / 2.0L);
 	}

math/s_clogf.c

@@ -68,12 +68,9 @@ __clogf (__complex__ float x)
 	  float absy2 = absy * absy;
 	  if (absy2 <= FLT_MIN * 2.0f)
 	    {
-#if __FLT_EVAL_METHOD__ == 0
+	      float force_underflow = absy2 * absy2;
-	      __real__ result = absy2 / 2.0f - absy2 * absy2 / 4.0f;
+	      __real__ result = absy2 / 2.0f;
-#else
+	      math_force_eval (force_underflow);
-	      volatile float force_underflow = absy2 * absy2 / 4.0f;
-	      __real__ result = absy2 / 2.0f - force_underflow;
-#endif
 	    }
 	  else
 	    __real__ result = __log1pf (absy2) / 2.0f;

math/s_clogl.c

@@ -74,7 +74,11 @@ __clogl (__complex__ long double x)
 	{
 	  long double absy2 = absy * absy;
 	  if (absy2 <= LDBL_MIN * 2.0L)
-	    __real__ result = absy2 / 2.0L - absy2 * absy2 / 4.0L;
+	    {
+	      long double force_underflow = absy2 * absy2;
+	      __real__ result = absy2 / 2.0L;
+	      math_force_eval (force_underflow);
+	    }
 	  else
 	    __real__ result = __log1pl (absy2) / 2.0L;
 	}

sysdeps/i386/fpu/libm-test-ulps

@@ -593,6 +593,102 @@ ifloat: 1
 ildouble: 2
 ldouble: 2
 
+Function: Real part of "clog10_downward":
+double: 2
+float: 2
+idouble: 2
+ifloat: 2
+ildouble: 4
+ldouble: 4
+
+Function: Imaginary part of "clog10_downward":
+double: 2
+float: 2
+idouble: 2
+ifloat: 2
+ildouble: 2
+ldouble: 2
+
+Function: Real part of "clog10_towardzero":
+double: 2
+float: 1
+idouble: 2
+ifloat: 1
+ildouble: 4
+ldouble: 4
+
+Function: Imaginary part of "clog10_towardzero":
+double: 2
+float: 2
+idouble: 2
+ifloat: 2
+ildouble: 2
+ldouble: 2
+
+Function: Real part of "clog10_upward":
+double: 2
+float: 2
+idouble: 2
+ifloat: 2
+ildouble: 4
+ldouble: 4
+
+Function: Imaginary part of "clog10_upward":
+double: 2
+float: 2
+idouble: 2
+ifloat: 2
+ildouble: 2
+ldouble: 2
+
+Function: Real part of "clog_downward":
+double: 2
+float: 1
+idouble: 2
+ifloat: 1
+ildouble: 2
+ldouble: 2
+
+Function: Imaginary part of "clog_downward":
+double: 1
+float: 1
+idouble: 1
+ifloat: 1
+ildouble: 1
+ldouble: 1
+
+Function: Real part of "clog_towardzero":
+double: 2
+float: 1
+idouble: 2
+ifloat: 1
+ildouble: 2
+ldouble: 2
+
+Function: Imaginary part of "clog_towardzero":
+double: 1
+float: 1
+idouble: 1
+ifloat: 1
+ildouble: 1
+ldouble: 1
+
+Function: Real part of "clog_upward":
+double: 1
+float: 1
+idouble: 1
+ifloat: 1
+ildouble: 2
+ldouble: 2
+
+Function: Imaginary part of "clog_upward":
+double: 1
+float: 1
+idouble: 1
+ifloat: 1
+ildouble: 1
+ldouble: 1
+
 Function: "cos":
 ildouble: 1
 ldouble: 1

sysdeps/x86_64/fpu/libm-test-ulps

@@ -624,6 +624,102 @@ ifloat: 1
 ildouble: 2
 ldouble: 2
 
+Function: Real part of "clog10_downward":
+double: 3
+float: 3
+idouble: 3
+ifloat: 3
+ildouble: 4
+ldouble: 4
+
+Function: Imaginary part of "clog10_downward":
+double: 3
+float: 2
+idouble: 3
+ifloat: 2
+ildouble: 2
+ldouble: 2
+
+Function: Real part of "clog10_towardzero":
+double: 3
+float: 2
+idouble: 3
+ifloat: 2
+ildouble: 4
+ldouble: 4
+
+Function: Imaginary part of "clog10_towardzero":
+double: 3
+float: 2
+idouble: 3
+ifloat: 2
+ildouble: 2
+ldouble: 2
+
+Function: Real part of "clog10_upward":
+double: 4
+float: 3
+idouble: 4
+ifloat: 3
+ildouble: 4
+ldouble: 4
+
+Function: Imaginary part of "clog10_upward":
+double: 2
+float: 2
+idouble: 2
+ifloat: 2
+ildouble: 2
+ldouble: 2
+
+Function: Real part of "clog_downward":
+double: 2
+float: 2
+idouble: 2
+ifloat: 2
+ildouble: 2
+ldouble: 2
+
+Function: Imaginary part of "clog_downward":
+double: 1
+float: 2
+idouble: 1
+ifloat: 2
+ildouble: 1
+ldouble: 1
+
+Function: Real part of "clog_towardzero":
+double: 2
+float: 2
+idouble: 2
+ifloat: 2
+ildouble: 2
+ldouble: 2
+
+Function: Imaginary part of "clog_towardzero":
+double: 1
+float: 2
+idouble: 1
+ifloat: 2
+ildouble: 1
+ldouble: 1
+
+Function: Real part of "clog_upward":
+double: 2
+float: 1
+idouble: 2
+ifloat: 1
+ildouble: 2
+ldouble: 2
+
+Function: Imaginary part of "clog_upward":
+double: 2
+float: 2
+idouble: 2
+ifloat: 2
+ildouble: 1
+ldouble: 1
+
 Function: "cos":
 ildouble: 1
 ldouble: 1