Fix exp10 inaccuracy and exceptions (bugs 13884, 13914).

11 files changed, 255 insertions, 62 deletions

diff --git a/ChangeLog b/ChangeLog
index 4109a7587b..cd40e73a46 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,5 +1,35 @@
 2012-05-06  Joseph Myers  <joseph@codesourcery.com>
 
+	[BZ #13884]
+	[BZ #13914]
+	* sysdeps/i386/fpu/e_expl.S (IEEE754_EXPL): Define conditional on
+	USE_AS_EXP10L.
+	(EXPL_FINITE): Likewise.
+	(FLDLOG): Likewise.
+	(c0): Likewise.
+	(c1): Likewise.
+	(__ieee754_expl): Change to IEEE754_EXPL.  Use FLDLOG macro.
+	Adjust comments for base varying.
+	(__expl_finite): Change alias to EXPL_FINITE.
+	* sysdeps/i386/fpu/e_exp10l.S: Define USE_AS_EXP10L and include
+	e_expl.S.
+	* sysdeps/ieee754/dbl-64/e_exp10.c: New file.
+	* sysdeps/ieee754/ldbl-128/e_exp10l.c: Likewise.
+	* sysdeps/ieee754/ldbl-128ibm/e_exp10l.c: Likewise.
+	* sysdeps/x86_64/fpu/e_exp10l.S: Likewise.
+	* sysdeps/x86_64/fpu/e_expl.S (IEEE754_EXPL): Define conditional on
+	USE_AS_EXP10L.
+	(EXPL_FINITE): Likewise.
+	(FLDLOG): Likewise.
+	(c0): Likewise.
+	(c1): Likewise.
+	(__ieee754_expl): Change to IEEE754_EXPL.  Use FLDLOG macro.
+	Adjust comments for base varying.
+	(__expl_finite): Change alias to EXPL_FINITE.
+	* math/libm-test.inc (exp10_test): Add more tests.  Do not disable
+	tests for bugs.
+	* sysdeps/x86_64/fpu/libm-test-ulps: Update.
+
 	[BZ #14064]
 	* math/libm-test.inc (check_float_internal): Correct ulp
 	calculation for subnormal expected results.
diff --git a/NEWS b/NEWS
index 66b1041485..90662f4fa0 100644
--- a/NEWS
+++ b/NEWS
@@ -20,10 +20,11 @@ Version 2.16
   13583, 13592, 13618, 13637, 13656, 13658, 13673, 13691, 13695, 13704,
   13705, 13706, 13726, 13738, 13739, 13758, 13760, 13761, 13775, 13786,
   13787, 13792, 13806, 13824, 13840, 13841, 13844, 13846, 13851, 13852,
-  13854, 13871, 13872, 13873, 13879, 13883, 13886, 13892, 13895, 13908,
-  13910, 13911, 13912, 13913, 13915, 13916, 13917, 13918, 13919, 13920,
-  13921, 13922, 13924, 13926, 13927, 13928, 13938, 13941, 13942, 13963,
-  13967, 13970, 13973, 14027, 14033, 14034, 14040, 14049, 14055, 14064
+  13854, 13871, 13872, 13873, 13879, 13883, 13884, 13886, 13892, 13895,
+  13908, 13910, 13911, 13912, 13913, 13914, 13915, 13916, 13917, 13918,
+  13919, 13920, 13921, 13922, 13924, 13926, 13927, 13928, 13938, 13941,
+  13942, 13963, 13967, 13970, 13973, 14027, 14033, 14034, 14040, 14049,
+  14055, 14064
 
 * ISO C11 support:
 
diff --git a/math/libm-test.inc b/math/libm-test.inc
index d65c3e5844..5fe9e5a14e 100644
--- a/math/libm-test.inc
+++ b/math/libm-test.inc
@@ -3464,15 +3464,20 @@ 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);
+  TEST_f_f (exp10, -36, 1.0e-36L);
+#ifndef TEST_FLOAT
+  TEST_f_f (exp10, 305, 1.0e305L);
+  TEST_f_f (exp10, -305, 1.0e-305L);
+#endif
+#if defined TEST_LDOUBLE && LDBL_MAX_10_EXP >= 4932
+  TEST_f_f (exp10, 4932, 1.0e4932L);
+  TEST_f_f (exp10, -4932, 1.0e-4932L);
 #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);
   TEST_f_f (exp10, -max_value, 0);
-#endif
   TEST_f_f (exp10, 0.75L, 5.62341325190349080394951039776481231L);
 
   END (exp10);
diff --git a/sysdeps/i386/fpu/e_exp10l.S b/sysdeps/i386/fpu/e_exp10l.S
index 04ec8001d9..d843e2b5e8 100644
--- a/sysdeps/i386/fpu/e_exp10l.S
+++ b/sysdeps/i386/fpu/e_exp10l.S
@@ -1,39 +1,2 @@
-/*
- * Written by Ulrich Drepper <drepper@cygnus.com>.
- */
-
-#include <machine/asm.h>
-
-/* 10^x = 2^(x * log2l(10)) */
-ENTRY(__ieee754_exp10l)
-	fldt	4(%esp)
-/* I added the following ugly construct because expl(+-Inf) resulted
-   in NaN.  The ugliness results from the bright minds at Intel.
-   For the i686 the code can be written better.
-   -- drepper@cygnus.com.  */
-	fxam				/* Is NaN or +-Inf?  */
-	fstsw	%ax
-	movb	$0x45, %dh
-	andb	%ah, %dh
-	cmpb	$0x05, %dh
-	je	1f			/* Is +-Inf, jump.  */
-	fldl2t
-	fmulp				/* x * log2(10) */
-	fld	%st
-	frndint				/* int(x * log2(10)) */
-	fsubr	%st,%st(1)		/* fract(x * log2(10)) */
-	fxch
-	f2xm1				/* 2^(fract(x * log2(10))) - 1 */
-	fld1
-	faddp				/* 2^(fract(x * log2(10))) */
-	fscale				/* e^x */
-	fstp	%st(1)
-	ret
-
-1:	testl	$0x200, %eax		/* Test sign.  */
-	jz	2f			/* If positive, jump.  */
-	fstp	%st
-	fldz				/* Set result to 0.  */
-2:	ret
-END (__ieee754_exp10l)
-strong_alias (__ieee754_exp10l, __exp10l_finite)
+#define USE_AS_EXP10L
+#include <e_expl.S>
diff --git a/sysdeps/i386/fpu/e_expl.S b/sysdeps/i386/fpu/e_expl.S
index 45c4d07539..9adf2a489e 100644
--- a/sysdeps/i386/fpu/e_expl.S
+++ b/sysdeps/i386/fpu/e_expl.S
@@ -24,9 +24,29 @@
 
 #include <machine/asm.h>
 
+#ifdef USE_AS_EXP10L
+# define IEEE754_EXPL __ieee754_exp10l
+# define EXPL_FINITE __exp10l_finite
+# define FLDLOG fldl2t
+#else
+# define IEEE754_EXPL __ieee754_expl
+# define EXPL_FINITE __expl_finite
+# define FLDLOG fldl2e
+#endif
+
 	.section .rodata.cst16,"aM",@progbits,16
 
 	.p2align 4
+#ifdef USE_AS_EXP10L
+	ASM_TYPE_DIRECTIVE(c0,@object)
+c0:	.byte 0, 0, 0, 0, 0, 0, 0x9a, 0xd4, 0x00, 0x40
+	.byte 0, 0, 0, 0, 0, 0
+	ASM_SIZE_DIRECTIVE(c0)
+	ASM_TYPE_DIRECTIVE(c1,@object)
+c1:	.byte 0x58, 0x92, 0xfc, 0x15, 0x37, 0x9a, 0x97, 0xf0, 0xef, 0x3f
+	.byte 0, 0, 0, 0, 0, 0
+	ASM_SIZE_DIRECTIVE(c1)
+#else
 	ASM_TYPE_DIRECTIVE(c0,@object)
 c0:	.byte 0, 0, 0, 0, 0, 0, 0xaa, 0xb8, 0xff, 0x3f
 	.byte 0, 0, 0, 0, 0, 0
@@ -35,6 +55,7 @@ c0:	.byte 0, 0, 0, 0, 0, 0, 0xaa, 0xb8, 0xff, 0x3f
 c1:	.byte 0x20, 0xfa, 0xee, 0xc2, 0x5f, 0x70, 0xa5, 0xec, 0xed, 0x3f
 	.byte 0, 0, 0, 0, 0, 0
 	ASM_SIZE_DIRECTIVE(c1)
+#endif
 	ASM_TYPE_DIRECTIVE(csat,@object)
 csat:	.byte 0, 0, 0, 0, 0, 0, 0, 0x80, 0x0e, 0x40
 	.byte 0, 0, 0, 0, 0, 0
@@ -47,7 +68,7 @@ csat:	.byte 0, 0, 0, 0, 0, 0, 0, 0x80, 0x0e, 0x40
 #endif
 
 	.text
-ENTRY(__ieee754_expl)
+ENTRY(IEEE754_EXPL)
 	fldt	4(%esp)
 /* I added the following ugly construct because expl(+-Inf) resulted
    in NaN.  The ugliness results from the bright minds at Intel.
@@ -75,8 +96,8 @@ ENTRY(__ieee754_expl)
 	andb	$2, %ah
 	jz	3f
 	fchs
-3:	fldl2e			/* 1  log2(e)         */
-	fmul	%st(1), %st	/* 1  x log2(e)       */
+3:	FLDLOG			/* 1  log2(base)      */
+	fmul	%st(1), %st	/* 1  x log2(base)    */
 	frndint			/* 1  i               */
 	fld	%st(1)		/* 2  x               */
 	frndint			/* 2  xi              */
@@ -92,11 +113,11 @@ ENTRY(__ieee754_expl)
 	fldt	MO(c1)		/* 4                  */
 	fmul	%st(4), %st	/* 4  c1 * x          */
 	faddp	%st, %st(1)	/* 3  f = f + c1 * x  */
-	f2xm1			/* 3 2^(fract(x * log2(e))) - 1 */
+	f2xm1			/* 3 2^(fract(x * log2(base))) - 1 */
 	fld1			/* 4 1.0              */
-	faddp			/* 3 2^(fract(x * log2(e))) */
+	faddp			/* 3 2^(fract(x * log2(base))) */
 	fstp	%st(1)		/* 2  */
-	fscale			/* 2 scale factor is st(1); e^x */
+	fscale			/* 2 scale factor is st(1); base^x */
 	fstp	%st(1)		/* 1  */
 	fstp	%st(1)		/* 0  */
 	jmp	2f
@@ -105,5 +126,5 @@ ENTRY(__ieee754_expl)
 	fstp	%st
 	fldz			/* Set result to 0.  */
 2:	ret
-END(__ieee754_expl)
-strong_alias (__ieee754_expl, __expl_finite)
+END(IEEE754_EXPL)
+strong_alias (IEEE754_EXPL, EXPL_FINITE)
diff --git a/sysdeps/ieee754/dbl-64/e_exp10.c b/sysdeps/ieee754/dbl-64/e_exp10.c
new file mode 100644
index 0000000000..eeb63a6986
--- /dev/null
+++ b/sysdeps/ieee754/dbl-64/e_exp10.c
@@ -0,0 +1,48 @@
+/* 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>
+
+static const double log10_high = 0x2.4d7637p0;
+static const double log10_low = 0x7.6aaa2b05ba95cp-28;
+
+double
+__ieee754_exp10 (double arg)
+{
+  int32_t lx;
+  double arg_high, arg_low;
+  double exp_high, exp_low;
+
+  if (!__finite (arg))
+    return __ieee754_exp (arg);
+  if (arg < DBL_MIN_10_EXP - DBL_DIG - 10)
+    return DBL_MIN * DBL_MIN;
+  else if (arg > DBL_MAX_10_EXP + 1)
+    return DBL_MAX * DBL_MAX;
+
+  GET_LOW_WORD (lx, arg);
+  lx &= 0xf8000000;
+  arg_high = arg;
+  SET_LOW_WORD (arg_high, lx);
+  arg_low = arg - arg_high;
+  exp_high = arg_high * log10_high;
+  exp_low = arg_high * log10_low + arg_low * M_LN10;
+  return __ieee754_exp (exp_high) * __ieee754_exp (exp_low);
+}
+strong_alias (__ieee754_exp10, __exp10_finite)
diff --git a/sysdeps/ieee754/ldbl-128/e_exp10l.c b/sysdeps/ieee754/ldbl-128/e_exp10l.c
new file mode 100644
index 0000000000..503c1de3e8
--- /dev/null
+++ b/sysdeps/ieee754/ldbl-128/e_exp10l.c
@@ -0,0 +1,47 @@
+/* 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>
+
+static const long double log10_high = 0x2.4d763776aaa2bp0L;
+static const long double log10_low = 0x5.ba95b58ae0b4c28a38a3fb3e7698p-60L;
+
+long double
+__ieee754_exp10l (long double arg)
+{
+  ieee854_long_double_shape_type u;
+  long double arg_high, arg_low;
+  long double exp_high, exp_low;
+
+  if (!__finitel (arg))
+    return __ieee754_expl (arg);
+  if (arg < LDBL_MIN_10_EXP - LDBL_DIG - 10)
+    return LDBL_MIN * LDBL_MIN;
+  else if (arg > LDBL_MAX_10_EXP + 1)
+    return LDBL_MAX * LDBL_MAX;
+
+  u.value = arg;
+  u.parts64.lsw &= 0xfe00000000000000LL;
+  arg_high = u.value;
+  arg_low = arg - arg_high;
+  exp_high = arg_high * log10_high;
+  exp_low = arg_high * log10_low + arg_low * M_LN10l;
+  return __ieee754_expl (exp_high) * __ieee754_expl (exp_low);
+}
+strong_alias (__ieee754_exp10l, __exp10l_finite)
diff --git a/sysdeps/ieee754/ldbl-128ibm/e_exp10l.c b/sysdeps/ieee754/ldbl-128ibm/e_exp10l.c
new file mode 100644
index 0000000000..b19e6c7f0c
--- /dev/null
+++ b/sysdeps/ieee754/ldbl-128ibm/e_exp10l.c
@@ -0,0 +1,46 @@
+/* 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>
+
+static const long double log10_high = 0x2.4d763776aaa2cp0L;
+static const long double log10_low = -0xf.a456a4a751f4b3d75c75c04c18p-56L;
+
+long double
+__ieee754_exp10l (long double arg)
+{
+  union ibm_extended_long_double u;
+  long double arg_high, arg_low;
+  long double exp_high, exp_low;
+
+  if (!__finitel (arg))
+    return __ieee754_expl (arg);
+  if (arg < LDBL_MIN_10_EXP - LDBL_DIG - 10)
+    return LDBL_MIN * LDBL_MIN;
+  else if (arg > LDBL_MAX_10_EXP + 1)
+    return LDBL_MAX * LDBL_MAX;
+
+  u.d = arg;
+  arg_high = u.dd[0];
+  arg_low = u.dd[1];
+  exp_high = arg_high * log10_high;
+  exp_low = arg_high * log10_low + arg_low * M_LN10l;
+  return __ieee754_expl (exp_high) * __ieee754_expl (exp_low);
+}
+strong_alias (__ieee754_exp10l, __exp10l_finite)
diff --git a/sysdeps/x86_64/fpu/e_exp10l.S b/sysdeps/x86_64/fpu/e_exp10l.S
new file mode 100644
index 0000000000..d843e2b5e8
--- /dev/null
+++ b/sysdeps/x86_64/fpu/e_exp10l.S
@@ -0,0 +1,2 @@
+#define USE_AS_EXP10L
+#include <e_expl.S>
diff --git a/sysdeps/x86_64/fpu/e_expl.S b/sysdeps/x86_64/fpu/e_expl.S
index d497b28973..fd613f91d3 100644
--- a/sysdeps/x86_64/fpu/e_expl.S
+++ b/sysdeps/x86_64/fpu/e_expl.S
@@ -24,9 +24,29 @@
 
 #include <machine/asm.h>
 
+#ifdef USE_AS_EXP10L
+# define IEEE754_EXPL __ieee754_exp10l
+# define EXPL_FINITE __exp10l_finite
+# define FLDLOG fldl2t
+#else
+# define IEEE754_EXPL __ieee754_expl
+# define EXPL_FINITE __expl_finite
+# define FLDLOG fldl2e
+#endif
+
 	.section .rodata.cst16,"aM",@progbits,16
 
 	.p2align 4
+#ifdef USE_AS_EXP10L
+	ASM_TYPE_DIRECTIVE(c0,@object)
+c0:	.byte 0, 0, 0, 0, 0, 0, 0x9a, 0xd4, 0x00, 0x40
+	.byte 0, 0, 0, 0, 0, 0
+	ASM_SIZE_DIRECTIVE(c0)
+	ASM_TYPE_DIRECTIVE(c1,@object)
+c1:	.byte 0x58, 0x92, 0xfc, 0x15, 0x37, 0x9a, 0x97, 0xf0, 0xef, 0x3f
+	.byte 0, 0, 0, 0, 0, 0
+	ASM_SIZE_DIRECTIVE(c1)
+#else
 	ASM_TYPE_DIRECTIVE(c0,@object)
 c0:	.byte 0, 0, 0, 0, 0, 0, 0xaa, 0xb8, 0xff, 0x3f
 	.byte 0, 0, 0, 0, 0, 0
@@ -35,6 +55,7 @@ c0:	.byte 0, 0, 0, 0, 0, 0, 0xaa, 0xb8, 0xff, 0x3f
 c1:	.byte 0x20, 0xfa, 0xee, 0xc2, 0x5f, 0x70, 0xa5, 0xec, 0xed, 0x3f
 	.byte 0, 0, 0, 0, 0, 0
 	ASM_SIZE_DIRECTIVE(c1)
+#endif
 	ASM_TYPE_DIRECTIVE(csat,@object)
 csat:	.byte 0, 0, 0, 0, 0, 0, 0, 0x80, 0x0e, 0x40
 	.byte 0, 0, 0, 0, 0, 0
@@ -47,7 +68,7 @@ csat:	.byte 0, 0, 0, 0, 0, 0, 0, 0x80, 0x0e, 0x40
 #endif
 
 	.text
-ENTRY(__ieee754_expl)
+ENTRY(IEEE754_EXPL)
 	fldt	8(%rsp)
 /* I added the following ugly construct because expl(+-Inf) resulted
    in NaN.  The ugliness results from the bright minds at Intel.
@@ -72,8 +93,8 @@ ENTRY(__ieee754_expl)
 	andb	$2, %ah
 	jz	3f
 	fchs
-3:	fldl2e			/* 1  log2(e)         */
-	fmul	%st(1), %st	/* 1  x log2(e)       */
+3:	FLDLOG			/* 1  log2(base)      */
+	fmul	%st(1), %st	/* 1  x log2(base)    */
 	frndint			/* 1  i               */
 	fld	%st(1)		/* 2  x               */
 	frndint			/* 2  xi              */
@@ -89,11 +110,11 @@ ENTRY(__ieee754_expl)
 	fldt	MO(c1)		/* 4                  */
 	fmul	%st(4), %st	/* 4  c1 * x          */
 	faddp	%st, %st(1)	/* 3  f = f + c1 * x  */
-	f2xm1			/* 3 2^(fract(x * log2(e))) - 1 */
+	f2xm1			/* 3 2^(fract(x * log2(base))) - 1 */
 	fld1			/* 4 1.0              */
-	faddp			/* 3 2^(fract(x * log2(e))) */
+	faddp			/* 3 2^(fract(x * log2(base))) */
 	fstp	%st(1)		/* 2  */
-	fscale			/* 2 scale factor is st(1); e^x */
+	fscale			/* 2 scale factor is st(1); base^x */
 	fstp	%st(1)		/* 1  */
 	fstp	%st(1)		/* 0  */
 	jmp	2f
@@ -102,5 +123,5 @@ ENTRY(__ieee754_expl)
 	fstp	%st
 	fldz			/* Set result to 0.  */
 2:	ret
-END(__ieee754_expl)
-strong_alias (__ieee754_expl, __expl_finite)
+END(IEEE754_EXPL)
+strong_alias (IEEE754_EXPL, EXPL_FINITE)
diff --git a/sysdeps/x86_64/fpu/libm-test-ulps b/sysdeps/x86_64/fpu/libm-test-ulps
index 70b091f547..f33c07f326 100644
--- a/sysdeps/x86_64/fpu/libm-test-ulps
+++ b/sysdeps/x86_64/fpu/libm-test-ulps
@@ -1287,6 +1287,12 @@ idouble: 2
 ifloat: 1
 ildouble: 1
 ldouble: 1
+Test "exp10 (-305) == 1.0e-305":
+double: 1
+idouble: 1
+Test "exp10 (-36) == 1.0e-36":
+double: 1
+idouble: 1
 Test "exp10 (0.75) == 5.62341325190349080394951039776481231":
 double: 1
 float: 1
@@ -1301,6 +1307,9 @@ idouble: 6
 ifloat: 2
 ildouble: 8
 ldouble: 8
+Test "exp10 (36) == 1.0e36":
+double: 1
+idouble: 1
 
 # exp_downward
 Test "exp_downward (1) == e":