[PATCH] BZ#13889: expl (709.75) wrongly overflows for ldbl-128ibm

[Adhemerval Zanella <azanella at linux dot vnet dot ibm dot com>]

Re-sending a patch on my backlog that I forgot to follow:

--

The patch increase the high value to check if expl overflows. Current high
mark value is not really correct, the algorithm accepts high values.

It also corrects spurious underflow exceptions generated with subnormal number
multiplication. I followed the dbl-64 implementation that does pretty much the
same with SET_RESTORE_ROUND macro.

And finally I also corrected the wrapper function that was using the ldbl-128
threshold values to check for overflow and underflow.

Tested on ppc32 and ppc64.

--

2013-03-04  Adhemerval Zanella  <azanella@linux.vnet.ibm.com>

	[BZ #13889]
	* sysdeps/ieee754/ldbl-128ibm/e_expl.c (__ieee754_expl): Increase the high
	value to check if expl overflow. Also corrects spurious underflow
	exceptions generated with subnormal number multiplication.
	* sysdeps/ieee754/ldbl-128ibm/w_expl.c (__expl): Fix threshold constants
	to check for underflow and overflow.
	* math/libm-test.inc: Add exp test.

--

diff --git a/math/libm-test.inc b/math/libm-test.inc
index 6ac3cd2..76dd686 100644
--- a/math/libm-test.inc
+++ b/math/libm-test.inc
@@ -4308,6 +4308,7 @@ exp_test (void)
   TEST_f_f (exp, 0.75L, 2.11700001661267466854536981983709561L);
   TEST_f_f (exp, 50.0L, 5184705528587072464087.45332293348538L);
   TEST_f_f (exp, 88.72269439697265625L, 3.40233126623160774937554134772290447915e38L);
+  TEST_f_f (exp, 709.75L, 1.73983687326416055769825271167383e+308L);
 #if defined TEST_LDOUBLE && __LDBL_MAX_EXP__ > 1024
   /* The result can only be represented in sane long double.  */
   TEST_f_f (exp, 1000.0L, 0.197007111401704699388887935224332313e435L);
diff --git a/sysdeps/ieee754/ldbl-128ibm/e_expl.c b/sysdeps/ieee754/ldbl-128ibm/e_expl.c
index 8236390..fa9408c 100644
--- a/sysdeps/ieee754/ldbl-128ibm/e_expl.c
+++ b/sysdeps/ieee754/ldbl-128ibm/e_expl.c
@@ -70,11 +70,11 @@
 static const long double C[] = {
 /* Smallest integer x for which e^x overflows.  */
 #define himark C[0]
- 709.08956571282405153382846025171462914L,
+ 709.78271289338399678773454114191496482L,
 
 /* Largest integer x for which e^x underflows.  */
 #define lomark C[1]
--744.44007192138121808966388925909996033L,
+-744.44007192138126231410729844608163411L,
 
 /* 3x2^96 */
 #define THREEp96 C[2]
@@ -138,14 +138,12 @@ __ieee754_expl (long double x)
   if (isless (x, himark) && isgreater (x, lomark))
     {
       int tval1, tval2, unsafe, n_i, exponent2;
-      long double x22, n, result, xl;
-      union ibm_extended_long_double ex2_u, scale_u;
+      long double x22, n, xl;
+      union ibm_extended_long_double ex2_u, scale_u, result;
       fenv_t oldenv;
 
       feholdexcept (&oldenv);
-#ifdef FE_TONEAREST
       fesetround (FE_TONEAREST);
-#endif
 
       n = __roundl (x*M_1_LN2);
       x = x-n*M_LN2_0;
@@ -166,7 +164,7 @@ __ieee754_expl (long double x)
 		* __expl_table[T_EXPL_RES2 + tval2];
       n_i = (int)n;
       /* 'unsafe' is 1 iff n_1 != 0.  */
-      unsafe = fabsl(n_i) >= -LDBL_MIN_EXP - 1;
+      unsafe = abs (n_i) >= -LDBL_MIN_EXP - 1;
       ex2_u.ieee.exponent += n_i >> unsafe;
       /* Fortunately, there are no subnormal lowpart doubles in
 	 __expl_table, only normal values and zeros.
@@ -204,7 +202,7 @@ __ieee754_expl (long double x)
       /* Return result.  */
       fesetenv (&oldenv);
 
-      result = x22 * ex2_u.d + ex2_u.d;
+      result.d = x22 * ex2_u.d + ex2_u.d;
 
       /* Now we can test whether the result is ultimate or if we are unsure.
 	 In the later case we should probably call a mpn based routine to give
@@ -235,10 +233,22 @@ __ieee754_expl (long double x)
 	    return __ieee754_expl_proc2 (origx);
 	  }
        */
-      if (!unsafe)
-	return result;
-      else
-	return result * scale_u.d;
+      if (unsafe)
+	{
+	  /* The scale with denormal number might generated undeflow for
+	     first high double multiplication. The exception holding is
+	     to avoid it.  */
+	  if (result.ieee.exponent + scale_u.ieee.exponent
+		< IBM_EXTENDED_LONG_DOUBLE_BIAS)
+	    {
+	      feholdexcept (&oldenv);
+	      result.d *= scale_u.d;
+	      fesetenv (&oldenv);
+	    }
+	  else
+	    result.d *= scale_u.d;
+	}
+      return result.d;
     }
   /* Exceptional cases:  */
   else if (isless (x, himark))
diff --git a/sysdeps/ieee754/ldbl-128ibm/w_expl.c b/sysdeps/ieee754/ldbl-128ibm/w_expl.c
index a5e72b2..df8ef37 100644
--- a/sysdeps/ieee754/ldbl-128ibm/w_expl.c
+++ b/sysdeps/ieee754/ldbl-128ibm/w_expl.c
@@ -1,6 +1,32 @@
-/* Looks like we can use ieee854 w_expl.c as is for IBM extended format. */
+#include <math.h>
+#include <math_private.h>
 #include <math_ldbl_opt.h>
-#undef weak_alias
-#define weak_alias(n,a)
-#include <sysdeps/ieee754/ldbl-128/w_expl.c>
+
+/*
+ * wrapper expl(x)
+ */
+static const long double
+o_threshold = 709.78271289338399678773454114191496482L,
+u_threshold = -744.44007192138126231410729844608163411L;
+
+long double __expl(long double x)	/* wrapper exp */
+{
+#ifdef _IEEE_LIBM
+  return __ieee754_expl(x);
+#else
+  long double z;
+  z = __ieee754_expl(x);
+  if (_LIB_VERSION == _IEEE_)
+    return z;
+  if (__finitel(x))
+    {
+      if (x >= o_threshold)
+	return __kernel_standard_l(x,x,206); /* exp overflow */
+      else if (x <= u_threshold)
+	return __kernel_standard_l(x,x,207); /* exp underflow */
+    }
+  return z;
+#endif
+}
+hidden_def (__expl)
 long_double_symbol (libm, __expl, expl);
-- 
1.7.1