This is the mail archive of the
libc-alpha@sourceware.org
mailing list for the glibc project.
Fix exp10 inaccuracy and exceptions (bugs 13884, 13914)
- From: "Joseph S. Myers" <joseph at codesourcery dot com>
- To: libc-alpha at sourceware dot org
- Date: Sun, 6 May 2012 01:12:21 +0000 (UTC)
- Subject: Fix exp10 inaccuracy and exceptions (bugs 13884, 13914)
Bug 13884 is that exp10 functions are inaccurate, and bug 13914 is
that the x86 version raises spurious exceptions for large arguments.
Inaccuracy results from not using larger intermediate precision for
internal calculations. For x86 and x86_64, it's natural to use
exactly the same logic and code as for expl, in e_expl.S, for the
higher precision, just with different constants for log2(10). This
patch accordingly makes e_expl.S usable for exp10l if USE_AS_EXP10L is
defined, so replacing the previous e_exp10l.S for x86 and providing
one for x86_64 that includes e_expl.S where previously a generic,
inaccurate C version would have been provided.
That covers the libc targets using ldbl-96; ports targets using
ldbl-96 appear to have their own versions of exp10l. flt-32 was fixed
previously (using double for internal calculations). This patch
provides C versions for dbl-64, ldbl-128 and ldbl-128ibm, all using
similar basic logic of splitting both the argument and log(10) into
high and low parts.
Tested x86 and x86_64 and ulps updated accordingly. I also built and
did limited testing for mips64 and powerpc to verify that the ldbl-128
and ldbl-128ibm versions built and did not produce wild results. For
ldbl-128ibm I see a failure
Failure: Test: exp10 (-305) == 1.0e-305
Result:
is: 9.99999999999999996282e-306 0x1.c16c5c52535750700000p-1014
should be: 9.99999999999999996282e-306 0x1.c16c5c52535750800000p-1014
difference: 4.94065645841246544177e-324 0x0.00000000000010000000p-1022
ulp : 35184372088832.0000
max.ulp : 0.0000
but the large ulps value looks like a bug in how libm-test.inc
calculates ulps for subnormals (filed as bug 14064), not an excessive
error from exp10l.
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.
diff --git a/math/libm-test.inc b/math/libm-test.inc
index 59192ed..3b692b9 100644
--- a/math/libm-test.inc
+++ b/math/libm-test.inc
@@ -3463,15 +3463,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 04ec800..d843e2b 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 45c4d07..9adf2a4 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 0000000..eeb63a6
--- /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 0000000..503c1de
--- /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 0000000..b19e6c7
--- /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 0000000..d843e2b
--- /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 d497b28..fd613f9 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 9a1e80e..44eeaff 100644
--- a/sysdeps/x86_64/fpu/libm-test-ulps
+++ b/sysdeps/x86_64/fpu/libm-test-ulps
@@ -1275,6 +1275,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
@@ -1289,6 +1295,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":
--
Joseph S. Myers
joseph@codesourcery.com