diff mbox

Fix powf overflow handling in non-nearest rounding mode [BZ #23961]

Message ID b5887546-4923-8286-9cb8-3bb077acef08@arm.com
State New, archived
Headers

Commit Message

Szabolcs Nagy Dec. 10, 2018, 3:08 p.m. UTC 
  The threshold value at which powf overflows depends on the rounding mode
and the current check did not take this into account. So when the result
was rounded away from zero it could become infinity without setting
errno to ERANGE.

Example: pow(0x1.7ac7cp+5, 23) is 0x1.fffffep+127 + 0.1633ulp

If the result goes above 0x1.fffffep+127 + 0.5ulp then errno is set,
which is fine in nearest rounding mode, but

  powf(0x1.7ac7cp+5, 23) is inf in upward rounding mode
  powf(-0x1.7ac7cp+5, 23) is -inf in downward rounding mode

and the previous implementation did not set errno in these cases.

The fix tries to avoid affecting the common code path or calling a
function that may introduce a stack frame, so float arithmetics is used
to check the rounding mode and the threshold is selected accordingly.

2018-12-10  Szabolcs Nagy  <szabolcs.nagy@arm.com>

	[BZ #23961]
	* math/auto-libm-test-in: Add new test case.
	* math/auto-libm-test-out-pow: Regenerated.
	* sysdeps/ieee754/flt-32/e_powf.c (__powf): Fix overflow check.

Comments

Joseph Myers Dec. 10, 2018, 10:15 p.m. UTC | #1 
On Mon, 10 Dec 2018, Szabolcs Nagy wrote:

> 2018-12-10  Szabolcs Nagy  <szabolcs.nagy@arm.com>
> 
> 	[BZ #23961]
> 	* math/auto-libm-test-in: Add new test case.
> 	* math/auto-libm-test-out-pow: Regenerated.
> 	* sysdeps/ieee754/flt-32/e_powf.c (__powf): Fix overflow check.

OK.

(It's possible some architecture with its own powf implementation could 
need to XFAIL the new tests if that architecture doesn't get overflow 
exactly right for those cases: the test infrastructure doesn't currently 
support marking a test as very close to the overflow threshold and then 
allowing appropriately for computed results on the other side of that 
threshold only being slightly inaccurate, so OK, as long as they handle 
errno / exceptions consistently with the returned value.)
diff mbox

Patch

diff --git a/math/auto-libm-test-in b/math/auto-libm-test-in
index 785328562f..c488e7225c 100644
--- a/math/auto-libm-test-in
+++ b/math/auto-libm-test-in
@@ -7154,6 +7154,8 @@  pow 0x1.430d4cp+0 0x5.0e462p+4
 pow 0x9.8b82ap-4 -0x1.99907ap+12
 pow 0xd.73035p-4 -0x1.47bb8p+8
 pow 0x1.059c76p+0 0x1.ff80bep+11
+pow 0x1.7ac7cp+5 23
+pow -0x1.7ac7cp+5 23
 
 sin 0
 sin -0
diff --git a/math/auto-libm-test-out-pow b/math/auto-libm-test-out-pow
index 4a9f1fc66a..09ec53e49e 100644
--- a/math/auto-libm-test-out-pow
+++ b/math/auto-libm-test-out-pow
@@ -44171,3 +44171,53 @@  pow 0x1.059c76p+0 0x1.ff80bep+11
 = pow tonearest ibm128 0x1.059c76p+0 0xf.fc05fp+8 : 0xf.ffe5535a38f9be648255c105d4p+124 : inexact-ok
 = pow towardzero ibm128 0x1.059c76p+0 0xf.fc05fp+8 : 0xf.ffe5535a38f9be648255c105d4p+124 : inexact-ok
 = pow upward ibm128 0x1.059c76p+0 0xf.fc05fp+8 : 0xf.ffe5535a38f9be648255c105d8p+124 : inexact-ok
+pow 0x1.7ac7cp+5 23
+= pow downward binary32 0x2.f58f8p+4 0x1.7p+4 : 0xf.fffffp+124 : inexact-ok
+= pow tonearest binary32 0x2.f58f8p+4 0x1.7p+4 : 0xf.fffffp+124 : inexact-ok
+= pow towardzero binary32 0x2.f58f8p+4 0x1.7p+4 : 0xf.fffffp+124 : inexact-ok
+= pow upward binary32 0x2.f58f8p+4 0x1.7p+4 : plus_infty : inexact-ok overflow errno-erange
+= pow downward binary64 0x2.f58f8p+4 0x1.7p+4 : 0xf.fffff29cf02e8p+124 : inexact-ok
+= pow tonearest binary64 0x2.f58f8p+4 0x1.7p+4 : 0xf.fffff29cf02fp+124 : inexact-ok
+= pow towardzero binary64 0x2.f58f8p+4 0x1.7p+4 : 0xf.fffff29cf02e8p+124 : inexact-ok
+= pow upward binary64 0x2.f58f8p+4 0x1.7p+4 : 0xf.fffff29cf02fp+124 : inexact-ok
+= pow downward intel96 0x2.f58f8p+4 0x1.7p+4 : 0xf.fffff29cf02eeecp+124 : inexact-ok
+= pow tonearest intel96 0x2.f58f8p+4 0x1.7p+4 : 0xf.fffff29cf02eeecp+124 : inexact-ok
+= pow towardzero intel96 0x2.f58f8p+4 0x1.7p+4 : 0xf.fffff29cf02eeecp+124 : inexact-ok
+= pow upward intel96 0x2.f58f8p+4 0x1.7p+4 : 0xf.fffff29cf02eeedp+124 : inexact-ok
+= pow downward m68k96 0x2.f58f8p+4 0x1.7p+4 : 0xf.fffff29cf02eeecp+124 : inexact-ok
+= pow tonearest m68k96 0x2.f58f8p+4 0x1.7p+4 : 0xf.fffff29cf02eeecp+124 : inexact-ok
+= pow towardzero m68k96 0x2.f58f8p+4 0x1.7p+4 : 0xf.fffff29cf02eeecp+124 : inexact-ok
+= pow upward m68k96 0x2.f58f8p+4 0x1.7p+4 : 0xf.fffff29cf02eeedp+124 : inexact-ok
+= pow downward binary128 0x2.f58f8p+4 0x1.7p+4 : 0xf.fffff29cf02eeec4a7cde7b5a2ep+124 : inexact-ok
+= pow tonearest binary128 0x2.f58f8p+4 0x1.7p+4 : 0xf.fffff29cf02eeec4a7cde7b5a2ep+124 : inexact-ok
+= pow towardzero binary128 0x2.f58f8p+4 0x1.7p+4 : 0xf.fffff29cf02eeec4a7cde7b5a2ep+124 : inexact-ok
+= pow upward binary128 0x2.f58f8p+4 0x1.7p+4 : 0xf.fffff29cf02eeec4a7cde7b5a2e8p+124 : inexact-ok
+= pow downward ibm128 0x2.f58f8p+4 0x1.7p+4 : 0xf.fffff29cf02eeec4a7cde7b5ap+124 : inexact-ok
+= pow tonearest ibm128 0x2.f58f8p+4 0x1.7p+4 : 0xf.fffff29cf02eeec4a7cde7b5a4p+124 : inexact-ok
+= pow towardzero ibm128 0x2.f58f8p+4 0x1.7p+4 : 0xf.fffff29cf02eeec4a7cde7b5ap+124 : inexact-ok
+= pow upward ibm128 0x2.f58f8p+4 0x1.7p+4 : 0xf.fffff29cf02eeec4a7cde7b5a4p+124 : inexact-ok
+pow -0x1.7ac7cp+5 23
+= pow downward binary32 -0x2.f58f8p+4 0x1.7p+4 : minus_infty : inexact-ok overflow errno-erange
+= pow tonearest binary32 -0x2.f58f8p+4 0x1.7p+4 : -0xf.fffffp+124 : inexact-ok
+= pow towardzero binary32 -0x2.f58f8p+4 0x1.7p+4 : -0xf.fffffp+124 : inexact-ok
+= pow upward binary32 -0x2.f58f8p+4 0x1.7p+4 : -0xf.fffffp+124 : inexact-ok
+= pow downward binary64 -0x2.f58f8p+4 0x1.7p+4 : -0xf.fffff29cf02fp+124 : inexact-ok
+= pow tonearest binary64 -0x2.f58f8p+4 0x1.7p+4 : -0xf.fffff29cf02fp+124 : inexact-ok
+= pow towardzero binary64 -0x2.f58f8p+4 0x1.7p+4 : -0xf.fffff29cf02e8p+124 : inexact-ok
+= pow upward binary64 -0x2.f58f8p+4 0x1.7p+4 : -0xf.fffff29cf02e8p+124 : inexact-ok
+= pow downward intel96 -0x2.f58f8p+4 0x1.7p+4 : -0xf.fffff29cf02eeedp+124 : inexact-ok
+= pow tonearest intel96 -0x2.f58f8p+4 0x1.7p+4 : -0xf.fffff29cf02eeecp+124 : inexact-ok
+= pow towardzero intel96 -0x2.f58f8p+4 0x1.7p+4 : -0xf.fffff29cf02eeecp+124 : inexact-ok
+= pow upward intel96 -0x2.f58f8p+4 0x1.7p+4 : -0xf.fffff29cf02eeecp+124 : inexact-ok
+= pow downward m68k96 -0x2.f58f8p+4 0x1.7p+4 : -0xf.fffff29cf02eeedp+124 : inexact-ok
+= pow tonearest m68k96 -0x2.f58f8p+4 0x1.7p+4 : -0xf.fffff29cf02eeecp+124 : inexact-ok
+= pow towardzero m68k96 -0x2.f58f8p+4 0x1.7p+4 : -0xf.fffff29cf02eeecp+124 : inexact-ok
+= pow upward m68k96 -0x2.f58f8p+4 0x1.7p+4 : -0xf.fffff29cf02eeecp+124 : inexact-ok
+= pow downward binary128 -0x2.f58f8p+4 0x1.7p+4 : -0xf.fffff29cf02eeec4a7cde7b5a2e8p+124 : inexact-ok
+= pow tonearest binary128 -0x2.f58f8p+4 0x1.7p+4 : -0xf.fffff29cf02eeec4a7cde7b5a2ep+124 : inexact-ok
+= pow towardzero binary128 -0x2.f58f8p+4 0x1.7p+4 : -0xf.fffff29cf02eeec4a7cde7b5a2ep+124 : inexact-ok
+= pow upward binary128 -0x2.f58f8p+4 0x1.7p+4 : -0xf.fffff29cf02eeec4a7cde7b5a2ep+124 : inexact-ok
+= pow downward ibm128 -0x2.f58f8p+4 0x1.7p+4 : -0xf.fffff29cf02eeec4a7cde7b5a4p+124 : inexact-ok
+= pow tonearest ibm128 -0x2.f58f8p+4 0x1.7p+4 : -0xf.fffff29cf02eeec4a7cde7b5a4p+124 : inexact-ok
+= pow towardzero ibm128 -0x2.f58f8p+4 0x1.7p+4 : -0xf.fffff29cf02eeec4a7cde7b5ap+124 : inexact-ok
+= pow upward ibm128 -0x2.f58f8p+4 0x1.7p+4 : -0xf.fffff29cf02eeec4a7cde7b5ap+124 : inexact-ok
diff --git a/sysdeps/ieee754/flt-32/e_powf.c b/sysdeps/ieee754/flt-32/e_powf.c
index d0c912188f..965f1fc2e7 100644
--- a/sysdeps/ieee754/flt-32/e_powf.c
+++ b/sysdeps/ieee754/flt-32/e_powf.c
@@ -17,6 +17,8 @@ 
    <http://www.gnu.org/licenses/>.  */
 
 #include <math.h>
+#include <math-barriers.h>
+#include <math-narrow-eval.h>
 #include <stdint.h>
 #include <shlib-compat.h>
 #include <libm-alias-float.h>
@@ -207,7 +209,17 @@  __powf (float x, float y)
     {
       /* |y*log(x)| >= 126.  */
       if (ylogx > 0x1.fffffffd1d571p+6 * POWF_SCALE)
+	/* |x^y| > 0x1.ffffffp127.  */
 	return __math_oflowf (sign_bias);
+      if (WANT_ROUNDING && WANT_ERRNO
+	  && ylogx > 0x1.fffffffa3aae2p+6 * POWF_SCALE)
+	/* |x^y| > 0x1.fffffep127, check if we round away from 0.  */
+	if ((!sign_bias
+	     && math_narrow_eval (1.0f + math_opt_barrier (0x1p-25f)) != 1.0f)
+	    || (sign_bias
+		&& math_narrow_eval (-1.0f - math_opt_barrier (0x1p-25f))
+		     != -1.0f))
+	  return __math_oflowf (sign_bias);
       if (ylogx <= -150.0 * POWF_SCALE)
 	return __math_uflowf (sign_bias);
 #if WANT_ERRNO_UFLOW