diff options
26 files changed, 885 insertions, 3672 deletions
@@ -16,7 +16,7 @@ Major new features: to set the install root if you wish to install into a non-default configured location. -* Optimized generic sinf, cosf, sincosf and tanf. +* Optimized generic exp, exp2, sinf, cosf, sincosf and tanf. * The reallocarray function is now declared under _DEFAULT_SOURCE, not just for _GNU_SOURCE, to match BSD environments. diff --git a/math/Makefile b/math/Makefile index bff726fe81..f1ba1a0c36 100644 --- a/math/Makefile +++ b/math/Makefile @@ -42,7 +42,7 @@ extra-libs-others = $(extra-libs) libm-support = s_lib_version s_matherr s_signgam \ fclrexcpt fgetexcptflg fraiseexcpt fsetexcptflg \ ftestexcept fegetround fesetround fegetenv feholdexcpt \ - fesetenv feupdateenv t_exp fedisblxcpt feenablxcpt \ + fesetenv feupdateenv fedisblxcpt feenablxcpt \ fegetexcept fesetexcept fetestexceptflag fegetmode \ fesetmode @@ -127,7 +127,7 @@ type-ldouble-yes := ldouble type-double-suffix := type-double-routines := branred doasin dosincos mpa mpatan2 \ k_rem_pio2 mpatan mpsqrt mptan sincos32 \ - sincostab + sincostab math_err e_exp_data # float support type-float-suffix := f diff --git a/sysdeps/aarch64/libm-test-ulps b/sysdeps/aarch64/libm-test-ulps index 6a89028459..585e5bbce7 100644 --- a/sysdeps/aarch64/libm-test-ulps +++ b/sysdeps/aarch64/libm-test-ulps @@ -716,9 +716,9 @@ ildouble: 2 ldouble: 2 Function: Imaginary part of "ccos_downward": -double: 2 +double: 3 float: 3 -idouble: 2 +idouble: 3 ifloat: 3 ildouble: 2 ldouble: 2 @@ -732,9 +732,9 @@ ildouble: 2 ldouble: 2 Function: Imaginary part of "ccos_towardzero": -double: 2 +double: 3 float: 3 -idouble: 2 +idouble: 3 ifloat: 3 ildouble: 2 ldouble: 2 @@ -772,33 +772,33 @@ ildouble: 1 ldouble: 1 Function: Real part of "ccosh_downward": -double: 1 +double: 2 float: 3 -idouble: 1 +idouble: 2 ifloat: 3 ildouble: 2 ldouble: 2 Function: Imaginary part of "ccosh_downward": -double: 2 +double: 3 float: 3 -idouble: 2 +idouble: 3 ifloat: 3 ildouble: 2 ldouble: 2 Function: Real part of "ccosh_towardzero": -double: 1 +double: 2 float: 3 -idouble: 1 +idouble: 2 ifloat: 3 ildouble: 2 ldouble: 2 Function: Imaginary part of "ccosh_towardzero": -double: 2 +double: 3 float: 3 -idouble: 2 +idouble: 3 ifloat: 3 ildouble: 2 ldouble: 2 @@ -836,33 +836,33 @@ ildouble: 1 ldouble: 1 Function: Real part of "cexp_downward": -double: 1 +double: 2 float: 2 -idouble: 1 +idouble: 2 ifloat: 2 ildouble: 2 ldouble: 2 Function: Imaginary part of "cexp_downward": -double: 1 +double: 3 float: 3 -idouble: 1 +idouble: 3 ifloat: 3 ildouble: 2 ldouble: 2 Function: Real part of "cexp_towardzero": -double: 1 +double: 2 float: 2 -idouble: 1 +idouble: 2 ifloat: 2 ildouble: 2 ldouble: 2 Function: Imaginary part of "cexp_towardzero": -double: 1 +double: 3 float: 3 -idouble: 1 +idouble: 3 ifloat: 3 ildouble: 2 ldouble: 2 @@ -876,9 +876,9 @@ ildouble: 3 ldouble: 3 Function: Imaginary part of "cexp_upward": -double: 1 +double: 3 float: 2 -idouble: 1 +idouble: 3 ifloat: 2 ildouble: 3 ldouble: 3 @@ -1052,25 +1052,25 @@ ildouble: 1 ldouble: 1 Function: "cosh_downward": -double: 1 +double: 2 float: 1 -idouble: 1 +idouble: 2 ifloat: 1 ildouble: 1 ldouble: 2 Function: "cosh_towardzero": -double: 1 +double: 2 float: 1 -idouble: 1 +idouble: 2 ifloat: 1 ildouble: 1 ldouble: 2 Function: "cosh_upward": -double: 1 +double: 2 float: 2 -idouble: 1 +idouble: 2 ifloat: 2 ildouble: 1 ldouble: 3 @@ -1090,9 +1090,9 @@ ildouble: 1 ldouble: 1 Function: Real part of "cpow_downward": -double: 4 +double: 5 float: 8 -idouble: 4 +idouble: 5 ifloat: 8 ildouble: 6 ldouble: 6 @@ -1106,9 +1106,9 @@ ildouble: 2 ldouble: 2 Function: Real part of "cpow_towardzero": -double: 4 +double: 5 float: 8 -idouble: 4 +idouble: 5 ifloat: 8 ildouble: 6 ldouble: 6 @@ -1150,9 +1150,9 @@ ildouble: 1 ldouble: 1 Function: Real part of "csin_downward": -double: 2 +double: 3 float: 3 -idouble: 2 +idouble: 3 ifloat: 3 ildouble: 2 ldouble: 2 @@ -1166,9 +1166,9 @@ ildouble: 2 ldouble: 2 Function: Real part of "csin_towardzero": -double: 2 +double: 3 float: 3 -idouble: 2 +idouble: 3 ifloat: 3 ildouble: 2 ldouble: 2 @@ -1220,9 +1220,9 @@ ildouble: 2 ldouble: 2 Function: Imaginary part of "csinh_downward": -double: 2 +double: 3 float: 3 -idouble: 2 +idouble: 3 ifloat: 3 ildouble: 2 ldouble: 2 @@ -1236,9 +1236,9 @@ ildouble: 2 ldouble: 2 Function: Imaginary part of "csinh_towardzero": -double: 2 +double: 3 float: 3 -idouble: 2 +idouble: 3 ifloat: 3 ildouble: 2 ldouble: 2 @@ -1492,9 +1492,9 @@ ildouble: 2 ldouble: 2 Function: "erfc_downward": -double: 3 +double: 4 float: 4 -idouble: 3 +idouble: 4 ifloat: 4 ildouble: 5 ldouble: 5 @@ -1508,9 +1508,9 @@ ildouble: 4 ldouble: 4 Function: "erfc_upward": -double: 3 +double: 4 float: 4 -idouble: 3 +idouble: 4 ifloat: 4 ildouble: 5 ldouble: 5 diff --git a/sysdeps/arm/libm-test-ulps b/sysdeps/arm/libm-test-ulps index e62cca3e03..798e3c465c 100644 --- a/sysdeps/arm/libm-test-ulps +++ b/sysdeps/arm/libm-test-ulps @@ -530,9 +530,9 @@ idouble: 1 ifloat: 1 Function: Imaginary part of "ccos_downward": -double: 2 +double: 3 float: 3 -idouble: 2 +idouble: 3 ifloat: 3 Function: Real part of "ccos_towardzero": @@ -542,9 +542,9 @@ idouble: 1 ifloat: 2 Function: Imaginary part of "ccos_towardzero": -double: 2 +double: 3 float: 3 -idouble: 2 +idouble: 3 ifloat: 3 Function: Real part of "ccos_upward": @@ -572,27 +572,27 @@ idouble: 1 ifloat: 1 Function: Real part of "ccosh_downward": -double: 1 +double: 2 float: 3 -idouble: 1 +idouble: 2 ifloat: 3 Function: Imaginary part of "ccosh_downward": -double: 2 +double: 3 float: 3 -idouble: 2 +idouble: 3 ifloat: 3 Function: Real part of "ccosh_towardzero": -double: 1 +double: 2 float: 3 -idouble: 1 +idouble: 2 ifloat: 3 Function: Imaginary part of "ccosh_towardzero": -double: 2 +double: 3 float: 3 -idouble: 2 +idouble: 3 ifloat: 3 Function: Real part of "ccosh_upward": @@ -620,27 +620,27 @@ idouble: 1 ifloat: 2 Function: Real part of "cexp_downward": -double: 1 +double: 2 float: 2 -idouble: 1 +idouble: 2 ifloat: 2 Function: Imaginary part of "cexp_downward": -double: 1 +double: 3 float: 3 -idouble: 1 +idouble: 3 ifloat: 3 Function: Real part of "cexp_towardzero": -double: 1 +double: 2 float: 2 -idouble: 1 +idouble: 2 ifloat: 2 Function: Imaginary part of "cexp_towardzero": -double: 1 +double: 3 float: 3 -idouble: 1 +idouble: 3 ifloat: 3 Function: Real part of "cexp_upward": @@ -650,9 +650,9 @@ idouble: 1 ifloat: 2 Function: Imaginary part of "cexp_upward": -double: 1 +double: 3 float: 2 -idouble: 1 +idouble: 3 ifloat: 2 Function: Real part of "clog": @@ -780,21 +780,21 @@ idouble: 1 ifloat: 1 Function: "cosh_downward": -double: 1 +double: 2 float: 1 -idouble: 1 +idouble: 2 ifloat: 1 Function: "cosh_towardzero": -double: 1 +double: 2 float: 1 -idouble: 1 +idouble: 2 ifloat: 1 Function: "cosh_upward": -double: 1 +double: 2 float: 2 -idouble: 1 +idouble: 2 ifloat: 2 Function: Real part of "cpow": @@ -808,9 +808,9 @@ float: 2 ifloat: 2 Function: Real part of "cpow_downward": -double: 4 +double: 5 float: 8 -idouble: 4 +idouble: 5 ifloat: 8 Function: Imaginary part of "cpow_downward": @@ -820,9 +820,9 @@ idouble: 1 ifloat: 2 Function: Real part of "cpow_towardzero": -double: 4 +double: 5 float: 8 -idouble: 4 +idouble: 5 ifloat: 8 Function: Imaginary part of "cpow_towardzero": @@ -850,9 +850,9 @@ idouble: 1 ifloat: 1 Function: Real part of "csin_downward": -double: 2 +double: 3 float: 3 -idouble: 2 +idouble: 3 ifloat: 3 Function: Imaginary part of "csin_downward": @@ -862,9 +862,9 @@ idouble: 1 ifloat: 1 Function: Real part of "csin_towardzero": -double: 2 +double: 3 float: 3 -idouble: 2 +idouble: 3 ifloat: 3 Function: Imaginary part of "csin_towardzero": @@ -902,9 +902,9 @@ idouble: 2 ifloat: 2 Function: Imaginary part of "csinh_downward": -double: 2 +double: 3 float: 3 -idouble: 2 +idouble: 3 ifloat: 3 Function: Real part of "csinh_towardzero": @@ -914,9 +914,9 @@ idouble: 2 ifloat: 2 Function: Imaginary part of "csinh_towardzero": -double: 2 +double: 3 float: 3 -idouble: 2 +idouble: 3 ifloat: 3 Function: Real part of "csinh_upward": @@ -1132,15 +1132,15 @@ double: 2 idouble: 2 Function: "exp10_downward": -double: 2 +double: 3 float: 1 -idouble: 2 +idouble: 3 ifloat: 1 Function: "exp10_towardzero": -double: 2 +double: 3 float: 1 -idouble: 2 +idouble: 3 ifloat: 1 Function: "exp10_upward": diff --git a/sysdeps/i386/fpu/e_exp_data.c b/sysdeps/i386/fpu/e_exp_data.c new file mode 100644 index 0000000000..1cc8931700 --- /dev/null +++ b/sysdeps/i386/fpu/e_exp_data.c @@ -0,0 +1 @@ +/* Not needed. */ diff --git a/sysdeps/i386/fpu/math_err.c b/sysdeps/i386/fpu/math_err.c new file mode 100644 index 0000000000..1cc8931700 --- /dev/null +++ b/sysdeps/i386/fpu/math_err.c @@ -0,0 +1 @@ +/* Not needed. */ diff --git a/sysdeps/i386/fpu/t_exp.c b/sysdeps/i386/fpu/t_exp.c deleted file mode 100644 index fd37963b05..0000000000 --- a/sysdeps/i386/fpu/t_exp.c +++ /dev/null @@ -1 +0,0 @@ -/* Empty. Not needed. */ diff --git a/sysdeps/ia64/fpu/e_exp_data.c b/sysdeps/ia64/fpu/e_exp_data.c new file mode 100644 index 0000000000..1cc8931700 --- /dev/null +++ b/sysdeps/ia64/fpu/e_exp_data.c @@ -0,0 +1 @@ +/* Not needed. */ diff --git a/sysdeps/ia64/fpu/math_err.c b/sysdeps/ia64/fpu/math_err.c new file mode 100644 index 0000000000..1cc8931700 --- /dev/null +++ b/sysdeps/ia64/fpu/math_err.c @@ -0,0 +1 @@ +/* Not needed. */ diff --git a/sysdeps/ia64/fpu/t_exp.c b/sysdeps/ia64/fpu/t_exp.c deleted file mode 100644 index 41254ae60a..0000000000 --- a/sysdeps/ia64/fpu/t_exp.c +++ /dev/null @@ -1 +0,0 @@ -/* Not needed. */ diff --git a/sysdeps/ieee754/dbl-64/e_exp.c b/sysdeps/ieee754/dbl-64/e_exp.c index 7d8b414034..209f20b972 100644 --- a/sysdeps/ieee754/dbl-64/e_exp.c +++ b/sysdeps/ieee754/dbl-64/e_exp.c @@ -1,48 +1,158 @@ -/* - * IBM Accurate Mathematical Library - * written by International Business Machines Corp. - * Copyright (C) 2001-2018 Free Software Foundation, Inc. - * - * This program 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. - * - * This program 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 this program; if not, see <http://www.gnu.org/licenses/>. - */ -/***************************************************************************/ -/* MODULE_NAME:uexp.c */ -/* */ -/* FUNCTION:uexp */ -/* exp1 */ -/* */ -/* FILES NEEDED:dla.h endian.h mpa.h mydefs.h uexp.h */ -/* */ -/* An ultimate exp routine. Given an IEEE double machine number x */ -/* it computes an almost correctly rounded (to nearest) value of e^x */ -/* Assumption: Machine arithmetic operations are performed in */ -/* round to nearest mode of IEEE 754 standard. */ -/* */ -/***************************************************************************/ +/* Double-precision e^x function. + Copyright (C) 2018 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 "endian.h" -#include "uexp.h" -#include "mydefs.h" -#include "MathLib.h" -#include "uexp.tbl" +#include <stdint.h> #include <math-barriers.h> -#include <math_private.h> -#include <fenv_private.h> -#include <fenv.h> -#include <float.h> -#include "eexp.tbl" +#include <math-narrow-eval.h> +#include "math_config.h" + +#define N (1 << EXP_TABLE_BITS) +#define InvLn2N __exp_data.invln2N +#define NegLn2hiN __exp_data.negln2hiN +#define NegLn2loN __exp_data.negln2loN +#define Shift __exp_data.shift +#define T __exp_data.tab +#define C2 __exp_data.poly[5 - EXP_POLY_ORDER] +#define C3 __exp_data.poly[6 - EXP_POLY_ORDER] +#define C4 __exp_data.poly[7 - EXP_POLY_ORDER] +#define C5 __exp_data.poly[8 - EXP_POLY_ORDER] + +/* Handle cases that may overflow or underflow when computing the result that + is scale*(1+TMP) without intermediate rounding. The bit representation of + scale is in SBITS, however it has a computed exponent that may have + overflown into the sign bit so that needs to be adjusted before using it as + a double. (int32_t)KI is the k used in the argument reduction and exponent + adjustment of scale, positive k here means the result may overflow and + negative k means the result may underflow. */ +static inline double +specialcase (double_t tmp, uint64_t sbits, uint64_t ki) +{ + double_t scale, y; + + if ((ki & 0x80000000) == 0) + { + /* k > 0, the exponent of scale might have overflowed by <= 460. */ + sbits -= 1009ull << 52; + scale = asdouble (sbits); + y = 0x1p1009 * (scale + scale * tmp); + return check_oflow (y); + } + /* k < 0, need special care in the subnormal range. */ + sbits += 1022ull << 52; + scale = asdouble (sbits); + y = scale + scale * tmp; + if (y < 1.0) + { + /* Round y to t |