diff options
author | Richard Maw <richard.maw@codethink.co.uk> | 2012-01-17 14:43:55 +0000 |
---|---|---|
committer | Richard Maw <richard.maw@codethink.co.uk> | 2012-01-17 14:43:55 +0000 |
commit | 2de9abc5c9d40b3c716307d67d16146f823fd554 (patch) | |
tree | 6979db67934ddc8b564150b465846a383b428ff8 /gnulib/lib/expl.c | |
parent | 33cc1c6fda6e72a7bae1401e9b2cec495a4d3ff1 (diff) | |
download | patch-baserock/bootstrap-pass2.tar.gz |
add the output of bootstrapbaserock/bootstrap-pass2baserock/bootstrap
Diffstat (limited to 'gnulib/lib/expl.c')
m--------- | gnulib | 0 | ||||
-rw-r--r-- | gnulib/lib/expl.c | 147 |
2 files changed, 147 insertions, 0 deletions
diff --git a/gnulib b/gnulib deleted file mode 160000 -Subproject 443bc5ffcf7429e557f4a371b0661abe98ddbc1 diff --git a/gnulib/lib/expl.c b/gnulib/lib/expl.c new file mode 100644 index 0000000..49fdefa --- /dev/null +++ b/gnulib/lib/expl.c @@ -0,0 +1,147 @@ +/* Emulation for expl. + Contributed by Paolo Bonzini + + Copyright 2002-2003, 2007, 2009-2011 Free Software Foundation, Inc. + + This file is part of gnulib. + + This program is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 3 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 General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program. If not, see <http://www.gnu.org/licenses/>. */ + +#include <config.h> + +/* Specification. */ +#include <math.h> + +#if HAVE_SAME_LONG_DOUBLE_AS_DOUBLE + +long double +expl (long double x) +{ + return exp (x); +} + +#else + +# include <float.h> + +static const long double C[] = { +/* Chebyshev polynom coeficients for (exp(x)-1)/x */ +# define P1 C[0] +# define P2 C[1] +# define P3 C[2] +# define P4 C[3] +# define P5 C[4] +# define P6 C[5] + 0.5L, + 1.66666666666666666666666666666666683E-01L, + 4.16666666666666666666654902320001674E-02L, + 8.33333333333333333333314659767198461E-03L, + 1.38888888889899438565058018857254025E-03L, + 1.98412698413981650382436541785404286E-04L, + +/* Smallest integer x for which e^x overflows. */ +# define himark C[6] + 11356.523406294143949491931077970765L, + +/* Largest integer x for which e^x underflows. */ +# define lomark C[7] +-11433.4627433362978788372438434526231L, + +/* very small number */ +# define TINY C[8] + 1.0e-4900L, + +/* 2^16383 */ +# define TWO16383 C[9] + 5.94865747678615882542879663314003565E+4931L}; + +long double +expl (long double x) +{ + /* Check for usual case. */ + if (x < himark && x > lomark) + { + int exponent; + long double t, x22; + int k = 1; + long double result = 1.0; + + /* Compute an integer power of e with a granularity of 0.125. */ + exponent = (int) floorl (x * 8.0L); + x -= exponent / 8.0L; + + if (x > 0.0625) + { + exponent++; + x -= 0.125L; + } + + if (exponent < 0) + { + t = 0.8824969025845954028648921432290507362220L; /* e^-0.25 */ + exponent = -exponent; + } + else + t = 1.1331484530668263168290072278117938725655L; /* e^0.25 */ + + while (exponent) + { + if (exponent & k) + { + result *= t; + exponent ^= k; + } + t *= t; + k <<= 1; + } + + /* Approximate (e^x - 1)/x, 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 + result * x22; + } + /* Exceptional cases: */ + else if (x < himark) + { + if (x + x == 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; +} + +#endif + +#if 0 +int +main (void) +{ + printf ("%.16Lg\n", expl (1.0L)); + printf ("%.16Lg\n", expl (-1.0L)); + printf ("%.16Lg\n", expl (2.0L)); + printf ("%.16Lg\n", expl (4.0L)); + printf ("%.16Lg\n", expl (-2.0L)); + printf ("%.16Lg\n", expl (-4.0L)); + printf ("%.16Lg\n", expl (0.0625L)); + printf ("%.16Lg\n", expl (0.3L)); + printf ("%.16Lg\n", expl (0.6L)); +} +#endif |