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

[Szabolcs Nagy <Szabolcs dot Nagy at arm dot com>]

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.

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