From 652df710ba782d47d40ff47fa699bd3d7fca57d6 Mon Sep 17 00:00:00 2001 From: Ulrich Drepper Date: Fri, 6 Aug 1999 17:20:19 +0000 Subject: Update. 1999-08-06 Jakub Jelinek * sysdeps/ieee754/ldbl-128/e_expl.c: New file. * sysdeps/ieee754/ldbl-128/t_expl.h: New file. * sysdeps/ieee754/ldbl-128/Dist: Add t_expl.h. --- sysdeps/ieee754/ldbl-128/e_expl.c | 248 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 248 insertions(+) create mode 100644 sysdeps/ieee754/ldbl-128/e_expl.c (limited to 'sysdeps/ieee754/ldbl-128/e_expl.c') diff --git a/sysdeps/ieee754/ldbl-128/e_expl.c b/sysdeps/ieee754/ldbl-128/e_expl.c new file mode 100644 index 0000000000..959d958736 --- /dev/null +++ b/sysdeps/ieee754/ldbl-128/e_expl.c @@ -0,0 +1,248 @@ +/* Quad-precision floating point e^x. + Copyright (C) 1999 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Jakub Jelinek + Partly based on double-precision code + by Geoffrey Keating + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 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 + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. 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 +#include +#include +#include +#include +#include +#include "t_expl.h" + +static const long double C[] = { +/* Smallest integer x for which e^x overflows. */ +#define himark C[0] + 11356.523406294143949491931077970765L, + +/* Largest integer x for which e^x underflows. */ +#define lomark C[1] +-11433.4627433362978788372438434526231L, + +/* 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-4900L, + +/* 2^16383 */ +#define TWO16383 C[9] + 5.94865747678615882542879663314003565E+4931L, + +/* 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; + long double x22, n, t, result, xl; + union ieee854_long_double ex2_u, scale_u; + fenv_t oldenv; + + feholdexcept (&oldenv); +#ifdef FE_TONEAREST + fesetround (FE_TONEAREST); +#endif + + /* Calculate n. */ + n = x * M_1_LN2 + THREEp111; + n -= THREEp111; + x = x - n * M_LN2_0; + xl = n * M_LN2_1; + + /* Calculate t/256. */ + t = x + THREEp103; + t -= THREEp103; + + /* Compute tval1 = t. */ + tval1 = (int) (t * TWO8); + + x -= __expl_table[T_EXPL_ARG1+2*tval1]; + xl -= __expl_table[T_EXPL_ARG1+2*tval1+1]; + + /* Calculate t/32768. */ + t = x + THREEp96; + t -= THREEp96; + + /* Compute tval2 = t. */ + tval2 = (int) (t * 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 = abs(n_i) >= -LDBL_MIN_EXP - 1; + ex2_u.ieee.exponent += n_i >> unsafe; + + /* 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 TWO16383*x; +} -- cgit 1.4.1