Fix exp2l inaccuracy (bug 13824).

2 files changed, 56 insertions, 3 deletions

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);
 }