about summary refs log tree commit diff
path: root/math/e_exp2l.c
diff options
context:
space:
mode:
authorJoseph Myers <joseph@codesourcery.com>2012-03-22 12:55:19 +0000
committerJoseph Myers <joseph@codesourcery.com>2012-03-22 12:55:19 +0000
commit48e44791e4d4d755bf7a7dd083d87584dc4779e4 (patch)
treeb35542729a07abdd56d4f502b73b5e561559665e /math/e_exp2l.c
parentc0df8e693f34b535bd6ee1b691bc4ca6bc3b4579 (diff)
downloadglibc-48e44791e4d4d755bf7a7dd083d87584dc4779e4.tar.gz
glibc-48e44791e4d4d755bf7a7dd083d87584dc4779e4.tar.xz
glibc-48e44791e4d4d755bf7a7dd083d87584dc4779e4.zip
Fix exp2l inaccuracy (bug 13824).
Diffstat (limited to 'math/e_exp2l.c')
-rw-r--r--math/e_exp2l.c44
1 files changed, 41 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)