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

This commit is contained in:
Joseph Myers 2012-05-06 18:23:44 +00:00
parent 6c23e11c4d
commit d8b82cad1b
11 changed files with 255 additions and 62 deletions

View File

@ -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.

9
NEWS
View File

@ -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:

View File

@ -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);

View File

@ -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>

View File

@ -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)

View File

@ -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)

View File

@ -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)

View File

@ -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)

View File

@ -0,0 +1,2 @@
#define USE_AS_EXP10L
#include <e_expl.S>

View File

@ -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)

View File

@ -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":