This is the mail archive of the libc-alpha@sourceware.org mailing list for the glibc project.


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]
Other format: [Raw text]

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


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.

---

2012-05-31  Adhemerval Zanella  <azanella@linux.vnet.ibm.com>

	[BZ #13889]
	* sysdeps/ieee754/ldbl-128ibm/e_expl.c: 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: Fix threshold constants to
	check for underflow and overflow.


diff --git a/sysdeps/ieee754/ldbl-128ibm/e_expl.c b/sysdeps/ieee754/ldbl-128ibm/e_expl.c
index 5618eb3..2c13735 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..c00b311 100644
--- a/sysdeps/ieee754/ldbl-128ibm/w_expl.c
+++ b/sysdeps/ieee754/ldbl-128ibm/w_expl.c
@@ -1,6 +1,44 @@
-/* Looks like we can use ieee854 w_expl.c as is for IBM extended format. */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+
+#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.6.0.2

-- 
Adhemerval Zanella Netto
  Software Engineer
  Linux Technology Center Brazil
  Toolchain / GLIBC on Power Architecture
  azanella@linux.vnet.ibm.com / azanella@br.ibm.com
  +55 61 8642-9890


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]