Fix exp10 spurious overflows (bug 13924).

4 files changed, 26 insertions, 16 deletions

diff --git a/math/e_exp10.c b/math/e_exp10.c
index 7518583..ce319fd 100644
--- a/math/e_exp10.c
+++ b/math/e_exp10.c
@@ -1,4 +1,4 @@
-/* Copyright (C) 1998, 2011 Free Software Foundation, Inc.
+/* Copyright (C) 1998-2012 Free Software Foundation, Inc.
    This file is part of the GNU C Library.
    Contributed by Ulrich Drepper <drepper@cygnus.com>, 1998.
 
@@ -18,13 +18,16 @@
 
 #include <math.h>
 #include <math_private.h>
-
+#include <float.h>
 
 double
 __ieee754_exp10 (double arg)
 {
-  /* This is a very stupid and inprecise implementation.  It'll get
-     replaced sometime (soon?).  */
-  return __ieee754_exp (M_LN10 * arg);
+  if (__finite (arg) && arg < DBL_MIN_10_EXP - DBL_DIG - 10)
+    return DBL_MIN * DBL_MIN;
+  else
+    /* This is a very stupid and inprecise implementation.  It'll get
+       replaced sometime (soon?).  */
+    return __ieee754_exp (M_LN10 * arg);
 }
 strong_alias (__ieee754_exp10, __exp10_finite)
diff --git a/math/e_exp10f.c b/math/e_exp10f.c
index 7b94005..582824f 100644
--- a/math/e_exp10f.c
+++ b/math/e_exp10f.c
@@ -1,4 +1,4 @@
-/* Copyright (C) 1998, 2011 Free Software Foundation, Inc.
+/* Copyright (C) 1998-2012 Free Software Foundation, Inc.
    This file is part of the GNU C Library.
    Contributed by Ulrich Drepper <drepper@cygnus.com>, 1998.
 
@@ -23,8 +23,10 @@
 float
 __ieee754_exp10f (float arg)
 {
-  /* This is a very stupid and inprecise implementation.  It'll get
-     replaced sometime (soon?).  */
-  return __ieee754_expf (M_LN10 * arg);
+  /* The argument to exp needs to be calculated in enough precision
+     that the fractional part has as much precision as float, in
+     addition to the bits in the integer part; using double ensures
+     this.  */
+  return __ieee754_exp (M_LN10 * arg);
 }
 strong_alias (__ieee754_exp10f, __exp10f_finite)
diff --git a/math/e_exp10l.c b/math/e_exp10l.c
index e3dad0a..d57c8cf 100644
--- a/math/e_exp10l.c
+++ b/math/e_exp10l.c
@@ -1,4 +1,4 @@
-/* Copyright (C) 1998, 2011 Free Software Foundation, Inc.
+/* Copyright (C) 1998-2012 Free Software Foundation, Inc.
    This file is part of the GNU C Library.
    Contributed by Ulrich Drepper <drepper@cygnus.com>, 1998.
 
@@ -18,13 +18,16 @@
 
 #include <math.h>
 #include <math_private.h>
-
+#include <float.h>
 
 long double
 __ieee754_exp10l (long double arg)
 {
-  /* This is a very stupid and inprecise implementation.  It'll get
-     replaced sometime (soon?).  */
-  return __ieee754_expl (M_LN10l * arg);
+  if (__finitel (arg) && arg < LDBL_MIN_10_EXP - LDBL_DIG - 10)
+    return LDBL_MIN * LDBL_MIN;
+  else
+    /* This is a very stupid and inprecise implementation.  It'll get
+       replaced sometime (soon?).  */
+    return __ieee754_expl (M_LN10l * arg);
 }
 strong_alias (__ieee754_exp10l, __exp10l_finite)
diff --git a/math/libm-test.inc b/math/libm-test.inc
index 1b5d1c7..cd627cd 100644
--- a/math/libm-test.inc
+++ b/math/libm-test.inc
@@ -3196,12 +3196,14 @@ exp10_test (void)
   TEST_f_f (exp10, nan_value, nan_value);
   TEST_f_f (exp10, 3, 1000);
   TEST_f_f (exp10, -1, 0.1L);
+#ifdef TEST_FLOAT /* Bug 13884: inaccurate results except for float.  */
+  TEST_f_f (exp10, 36, 1.0e36L);
+#endif
   TEST_f_f (exp10, 1e6, plus_infty, OVERFLOW_EXCEPTION);
   TEST_f_f (exp10, -1e6, 0);
 #ifndef TEST_LDOUBLE /* Bug 13914: spurious exceptions.  */
   TEST_f_f (exp10, max_value, plus_infty, OVERFLOW_EXCEPTION);
-  /* Bug 13924: spurious OVERFLOW exception may be present.  */
-  TEST_f_f (exp10, -max_value, 0, OVERFLOW_EXCEPTION_OK);
+  TEST_f_f (exp10, -max_value, 0);
 #endif
   TEST_f_f (exp10, 0.75L, 5.62341325190349080394951039776481231L);