diff options
Diffstat (limited to 'sysdeps/ieee754/ldbl-128ibm/e_expl.c')
| -rw-r--r-- | sysdeps/ieee754/ldbl-128ibm/e_expl.c | 257 |
1 files changed, 257 insertions, 0 deletions
diff --git a/sysdeps/ieee754/ldbl-128ibm/e_expl.c b/sysdeps/ieee754/ldbl-128ibm/e_expl.c new file mode 100644 index 0000000..3c4088f --- /dev/null +++ b/sysdeps/ieee754/ldbl-128ibm/e_expl.c @@ -0,0 +1,257 @@ +/* Quad-precision floating point e^x. + Copyright (C) 1999,2004,2006 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Jakub Jelinek <jj@ultra.linux.cz> + Partly based on double-precision code + by Geoffrey Keating <geoffk@ozemail.com.au> + + 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, write to the Free + Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA + 02111-1307 USA. */ + +/* The basic design here is from + Abraham Ziv, "Fast Evaluation of Elementary Mathematical Functions with + Correctly Rounded Last Bit", ACM Trans. Math. Soft., 17 (3), September 1991, + pp. 410-423. + + We work with number pairs where the first number is the high part and + the second one is the low part. Arithmetic with the high part numbers must + be exact, without any roundoff errors. + + The input value, X, is written as + X = n * ln(2)_0 + arg1[t1]_0 + arg2[t2]_0 + x + - n * ln(2)_1 + arg1[t1]_1 + arg2[t2]_1 + xl + + where: + - n is an integer, 16384 >= n >= -16495; + - ln(2)_0 is the first 93 bits of ln(2), and |ln(2)_0-ln(2)-ln(2)_1| < 2^-205 + - t1 is an integer, 89 >= t1 >= -89 + - t2 is an integer, 65 >= t2 >= -65 + - |arg1[t1]-t1/256.0| < 2^-53 + - |arg2[t2]-t2/32768.0| < 2^-53 + - x + xl is whatever is left, |x + xl| < 2^-16 + 2^-53 + + Then e^x is approximated as + + e^x = 2^n_1 ( 2^n_0 e^(arg1[t1]_0 + arg1[t1]_1) e^(arg2[t2]_0 + arg2[t2]_1) + + 2^n_0 e^(arg1[t1]_0 + arg1[t1]_1) e^(arg2[t2]_0 + arg2[t2]_1) + * p (x + xl + n * ln(2)_1)) + where: + - p(x) is a polynomial approximating e(x)-1 + - e^(arg1[t1]_0 + arg1[t1]_1) is obtained from a table + - e^(arg2[t2]_0 + arg2[t2]_1) likewise + - n_1 + n_0 = n, so that |n_0| < -LDBL_MIN_EXP-1. + + If it happens that n_1 == 0 (this is the usual case), that multiplication + is omitted. + */ + +#ifndef _GNU_SOURCE +#define _GNU_SOURCE +#endif +#include <float.h> +#include <ieee754.h> +#include <math.h> +#include <fenv.h> +#include <inttypes.h> +#include <math_private.h> +#include <sysdeps/ieee754/ldbl-128/t_expl.h> + +static const long double C[] = { +/* Smallest integer x for which e^x overflows. */ +#define himark C[0] + 709.08956571282405153382846025171462914L, + +/* Largest integer x for which e^x underflows. */ +#define lomark C[1] +-709.08956571282405153382846025171462914L, + +/* 3x2^96 */ +#define THREEp96 C[2] + 59421121885698253195157962752.0L, + +/* 3x2^103 */ +#define THREEp103 C[3] + 30423614405477505635920876929024.0L, + +/* 3x2^111 */ +#define THREEp111 C[4] + 7788445287802241442795744493830144.0L, + +/* 1/ln(2) */ +#define M_1_LN2 C[5] + 1.44269504088896340735992468100189204L, + +/* first 93 bits of ln(2) */ +#define M_LN2_0 C[6] + 0.693147180559945309417232121457981864L, + +/* ln2_0 - ln(2) */ +#define M_LN2_1 C[7] +-1.94704509238074995158795957333327386E-31L, + +/* very small number */ +#define TINY C[8] + 1.0e-308L, + +/* 2^16383 */ +#define TWO1023 C[9] + 8.988465674311579538646525953945123668E+307L, + +/* 256 */ +#define TWO8 C[10] + 256.0L, + +/* 32768 */ +#define TWO15 C[11] + 32768.0L, + +/* Chebyshev polynom coeficients for (exp(x)-1)/x */ +#define P1 C[12] +#define P2 C[13] +#define P3 C[14] +#define P4 C[15] +#define P5 C[16] +#define P6 C[17] + 0.5L, + 1.66666666666666666666666666666666683E-01L, + 4.16666666666666666666654902320001674E-02L, + 8.33333333333333333333314659767198461E-03L, + 1.38888888889899438565058018857254025E-03L, + 1.98412698413981650382436541785404286E-04L, +}; + +long double +__ieee754_expl (long double x) +{ + /* Check for usual case. */ + if (isless (x, himark) && isgreater (x, lomark)) + { + int tval1, tval2, unsafe, n_i, exponent2; + long double x22, n, result, xl; + union ibm_extended_long_double ex2_u, scale_u; + fenv_t oldenv; + + feholdexcept (&oldenv); +#ifdef FE_TONEAREST + fesetround (FE_TONEAREST); +#endif + + n = roundl(x*M_1_LN2); + x = x-n*M_LN2_0; + xl = n*M_LN2_1; + + tval1 = roundl(x*TWO8); + x -= __expl_table[T_EXPL_ARG1+2*tval1]; + xl -= __expl_table[T_EXPL_ARG1+2*tval1+1]; + + tval2 = roundl(x*TWO15); + x -= __expl_table[T_EXPL_ARG2+2*tval2]; + xl -= __expl_table[T_EXPL_ARG2+2*tval2+1]; + + x = x + xl; + + /* Compute ex2 = 2^n_0 e^(argtable[tval1]) e^(argtable[tval2]). */ + ex2_u.d = __expl_table[T_EXPL_RES1 + tval1] + * __expl_table[T_EXPL_RES2 + tval2]; + n_i = (int)n; + /* 'unsafe' is 1 iff n_1 != 0. */ + unsafe = fabsl(n_i) >= -LDBL_MIN_EXP - 1; + ex2_u.ieee.exponent += n_i >> unsafe; + /* Fortunately, there are no subnormal lowpart doubles in + __expl_table, only normal values and zeros. + But after scaling it can be subnormal. */ + exponent2 = ex2_u.ieee.exponent2 + (n_i >> unsafe); + if (ex2_u.ieee.exponent2 == 0) + /* assert ((ex2_u.ieee.mantissa2|ex2_u.ieee.mantissa3) == 0) */; + else if (exponent2 > 0) + ex2_u.ieee.exponent2 = exponent2; + else if (exponent2 <= -54) + { + ex2_u.ieee.exponent2 = 0; + ex2_u.ieee.mantissa2 = 0; + ex2_u.ieee.mantissa3 = 0; + } + else + { + static const double + two54 = 1.80143985094819840000e+16, /* 4350000000000000 */ + twom54 = 5.55111512312578270212e-17; /* 3C90000000000000 */ + ex2_u.dd[1] *= two54; + ex2_u.ieee.exponent2 += n_i >> unsafe; + ex2_u.dd[1] *= twom54; + } + + /* Compute scale = 2^n_1. */ + scale_u.d = 1.0L; + scale_u.ieee.exponent += n_i - (n_i >> unsafe); + + /* Approximate e^x2 - 1, using a seventh-degree polynomial, + with maximum error in [-2^-16-2^-53,2^-16+2^-53] + less than 4.8e-39. */ + x22 = x + x*x*(P1+x*(P2+x*(P3+x*(P4+x*(P5+x*P6))))); + + /* Return result. */ + fesetenv (&oldenv); + + result = x22 * ex2_u.d + ex2_u.d; + + /* Now we can test whether the result is ultimate or if we are unsure. + In the later case we should probably call a mpn based routine to give + the ultimate result. + Empirically, this routine is already ultimate in about 99.9986% of + cases, the test below for the round to nearest case will be false + in ~ 99.9963% of cases. + Without proc2 routine maximum error which has been seen is + 0.5000262 ulp. + + union ieee854_long_double ex3_u; + + #ifdef FE_TONEAREST + fesetround (FE_TONEAREST); + #endif + ex3_u.d = (result - ex2_u.d) - x22 * ex2_u.d; + ex2_u.d = result; + ex3_u.ieee.exponent += LDBL_MANT_DIG + 15 + IEEE854_LONG_DOUBLE_BIAS + - ex2_u.ieee.exponent; + n_i = abs (ex3_u.d); + n_i = (n_i + 1) / 2; + fesetenv (&oldenv); + #ifdef FE_TONEAREST + if (fegetround () == FE_TONEAREST) + n_i -= 0x4000; + #endif + if (!n_i) { + return __ieee754_expl_proc2 (origx); + } + */ + if (!unsafe) + return result; + else + return result * scale_u.d; + } + /* Exceptional cases: */ + else if (isless (x, himark)) + { + if (__isinfl (x)) + /* e^-inf == 0, with no error. */ + return 0; + else + /* Underflow */ + return TINY * TINY; + } + else + /* Return x, if x is a NaN or Inf; or overflow, otherwise. */ + return TWO1023*x; +} |