Fix exp2l inaccuracy (bug 13824)

["Joseph S. Myers" <joseph at codesourcery dot com>]

I propose this patch to fix bug 13824, inaccuracy of the default exp2l
implementation (used for ldbl-128 and ldbl-128ibm).  This patch
handles overflow and underflow (and infinity and NaN arguments) much
as they are handled in the dbl-64 exp2, then for in-range arguments
splits the integer and fractional parts, computes expl of M_LN2L times
the fractional part and scales that exactly with scalbnl.  This (a)
avoids inaccurate results for integer arguments, as long as expl (0)
is exactly 1, and (b) avoids the large errors (hundreds of ulps) for
large arguments that arise from bits lost in the multiplication M_LN2l
* x being significant for the mantissa after exponentiation (this
problem arises for exponentials of large inexact numbers, not for
numbers with small absolute value).

Tested powerpc, x86 and x86_64 (the latter two in case of any ulps
updates needed for the new tests).  No ulps updates were needed for
any of those architectures.

2012-03-21  Joseph Myers  <joseph@codesourcery.com>

	[BZ #13824]
	* math/e_exp2l.c: Include <float.h>.
	(__ieee754_exp2l): Handle overflow and underflow cases
	separately.  Only pass fractional part of argument to
	__ieee754_expl.
	* math/libm-test.inc (exp2_test): Add more tests.

diff --git a/math/e_exp2l.c b/math/e_exp2l.c
index e7e4939..8904d3e 100644
--- a/math/e_exp2l.c
+++ b/math/e_exp2l.c
@@ -1,11 +1,49 @@
+/* Compute 2^x.
+   Copyright (C) 2012 Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
+
+   The GNU C Library is free software; you can redistribute it and/or
+   modify it under the terms of the GNU Lesser General Public
+   License as published by the Free Software Foundation; either
+   version 2.1 of the License, or (at your option) any later version.
+
+   The GNU C Library is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+   Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public
+   License along with the GNU C Library; if not, see
+   <http://www.gnu.org/licenses/>.  */
+
 #include <math.h>
 #include <math_private.h>
+#include <float.h>
 
 long double
 __ieee754_exp2l (long double x)
 {
-  /* This is a very stupid and inprecise implementation.  It'll get
-     replaced sometime (soon?).  */
-  return __ieee754_expl (M_LN2l * x);
+  if (__builtin_expect (isless (x, (long double) LDBL_MAX_EXP), 1))
+    {
+      if (__builtin_expect (isgreaterequal (x, (long double) (LDBL_MIN_EXP
+							      - LDBL_MANT_DIG
+							      - 1)), 1))
+	{
+	  int intx = (int) x;
+	  long double fractx = x - intx;
+	  return __scalbnl (__ieee754_expl (M_LN2l * fractx), intx);
+	}
+      else
+	{
+	  /* Underflow or exact zero.  */
+	  if (__isinfl (x))
+	    return 0;
+	  else
+	    return LDBL_MIN * LDBL_MIN;
+	}
+    }
+  else
+    /* Infinity, NaN or overflow.  */
+    return LDBL_MAX * x;
 }
 strong_alias (__ieee754_exp2l, __exp2l_finite)
diff --git a/math/libm-test.inc b/math/libm-test.inc
index 05a000e..fad767d 100644
--- a/math/libm-test.inc
+++ b/math/libm-test.inc
@@ -3127,6 +3127,21 @@ exp2_test (void)
   TEST_f_f (exp2, -1e6, 0);
   TEST_f_f (exp2, 0.75L, 1.68179283050742908606225095246642979L);
 
+  TEST_f_f (exp2, 100.5, 1.792728671193156477399422023278661496394e+30L);
+  TEST_f_f (exp2, 127, 0x1p127);
+  TEST_f_f (exp2, -149, 0x1p-149);
+
+#ifndef TEST_FLOAT
+  TEST_f_f (exp2, 1000.25, 1.274245659452564874772384918171765416737e+301L);
+  TEST_f_f (exp2, 1023, 0x1p1023);
+  TEST_f_f (exp2, -1074, 0x1p-1074);
+#endif
+
+#if defined TEST_LDOUBLE && LDBL_MAX_EXP >= 16384
+  TEST_f_f (exp2, 16383, 0x1p16383L);
+  TEST_f_f (exp2, -16400, 0x1p-16400L);
+#endif
+
   END (exp2);
 }
 

-- 
Joseph S. Myers
joseph@codesourcery.com