From 48e44791e4d4d755bf7a7dd083d87584dc4779e4 Mon Sep 17 00:00:00 2001 From: Joseph Myers Date: Thu, 22 Mar 2012 12:55:19 +0000 Subject: Fix exp2l inaccuracy (bug 13824). --- math/e_exp2l.c | 44 +++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 41 insertions(+), 3 deletions(-) (limited to 'math/e_exp2l.c') 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 + . */ + #include #include +#include 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) -- cgit v1.1