Fix exp2l inaccuracy (bug 13824).

4 files changed, 69 insertions, 9 deletions

diff --git a/ChangeLog b/ChangeLog
index a993315e9e..599c340ce4 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,5 +1,12 @@
 2012-03-22  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.
+
 	* sysdeps/ieee754/ldbl-128/k_cosl.c (__kernel_cosl): Negate y if
 	negating x to take absolute value.
 	* sysdeps/ieee754/ldbl-128/k_sincosl.c (__kernel_sincosl):
diff --git a/NEWS b/NEWS
index 78fac7b732..5f9648e553 100644
--- a/NEWS
+++ b/NEWS
@@ -11,12 +11,12 @@ Version 2.16
 
   174, 350, 411, 2541, 2547, 2548, 2551, 2552, 2553, 2554, 2562, 2563, 2565,
   2566, 2576, 3335, 3976, 3992, 4026, 4108, 4596, 4822, 5077, 5461, 5805,
-  5993, 6471, 6730, 6884, 6907, 6911, 9739, 9902, 10110, 10135, 10140, 10210,
-  10545, 10716, 11174, 11322, 11365, 11451, 11494, 12047, 13058, 13525,
-  13526, 13527, 13528, 13529, 13530, 13531, 13532, 13533, 13547, 13551,
-  13552, 13553, 13555, 13559, 13566, 13583, 13618, 13637, 13656, 13658,
-  13673, 13695, 13704, 13706, 13726, 13738, 13786, 13792, 13806, 13840,
-  13841, 13844, 13846, 13851, 13852, 13854, 13871, 13883
+  5993, 6471, 6730, 6884, 6907, 6911, 9739, 9902, 10110, 10135, 10140,
+  10210, 10545, 10716, 11174, 11322, 11365, 11451, 11494, 12047, 13058,
+  13525, 13526, 13527, 13528, 13529, 13530, 13531, 13532, 13533, 13547,
+  13551, 13552, 13553, 13555, 13559, 13566, 13583, 13618, 13637, 13656,
+  13658, 13673, 13695, 13704, 13706, 13726, 13738, 13786, 13792, 13806,
+  13824, 13840, 13841, 13844, 13846, 13851, 13852, 13854, 13871, 13883
 
 * ISO C11 support:
 
diff --git a/math/e_exp2l.c b/math/e_exp2l.c
index e7e493933a..8904d3e4a0 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 05a000e0c1..fad767dd02 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);
 }