From 35eb879e3b7849e86f06658bdb888f7858b30759 Mon Sep 17 00:00:00 2001
From: Andreas Schwab <schwab@linux-m68k.org>
Date: Wed, 6 Jun 2012 18:00:52 +0200
Subject: [PATCH] sysdeps/ieee754/ldbl-128ibm/e_expl.c

---
 sysdeps/ieee754/ldbl-128ibm/e_expl.c | 34 ++++++++++++++++----------
 sysdeps/ieee754/ldbl-128ibm/w_expl.c | 46 ++++++++++++++++++++++++++++++++----
 2 files changed, 64 insertions(+), 16 deletions(-)

diff --git a/sysdeps/ieee754/ldbl-128ibm/e_expl.c b/sysdeps/ieee754/ldbl-128ibm/e_expl.c
index 8236390..60229d5 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 underflow 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);
-- 
2.7.4