-/* Emulation for expl.
- Contributed by Paolo Bonzini
+/* Exponential function.
+ Copyright (C) 2011-2014 Free Software Foundation, Inc.
- Copyright 2002, 2003 Free Software Foundation, Inc.
+ 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 file is part of gnulib.
+ 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.
- gnulib 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, or (at your option)
- any later version.
+ You should have received a copy of the GNU General Public License
+ along with this program. If not, see <http://www.gnu.org/licenses/>. */
- gnulib 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.
+#include <config.h>
- You should have received a copy of the GNU Lesser General Public License
- along with gnulib; see the file COPYING.LIB. If not, write to the Free
- Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307,
- USA.
- */
-
-#include <float.h>
+/* Specification. */
#include <math.h>
-#include "mathl.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};
+#if HAVE_SAME_LONG_DOUBLE_AS_DOUBLE
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;
+ return exp (x);
}
-#if 0
-int
-main ()
+#else
+
+# include <float.h>
+
+/* gl_expl_table[i] = exp((i - 128) * log(2)/256). */
+extern const long double gl_expl_table[257];
+
+/* A value slightly larger than log(2). */
+#define LOG2_PLUS_EPSILON 0.6931471805599454L
+
+/* Best possible approximation of log(2) as a 'long double'. */
+#define LOG2 0.693147180559945309417232121458176568075L
+
+/* Best possible approximation of 1/log(2) as a 'long double'. */
+#define LOG2_INVERSE 1.44269504088896340735992468100189213743L
+
+/* Best possible approximation of log(2)/256 as a 'long double'. */
+#define LOG2_BY_256 0.00270760617406228636491106297444600221904L
+
+/* Best possible approximation of 256/log(2) as a 'long double'. */
+#define LOG2_BY_256_INVERSE 369.329930467574632284140718336484387181L
+
+/* The upper 32 bits of log(2)/256. */
+#define LOG2_BY_256_HI_PART 0.0027076061733168899081647396087646484375L
+/* log(2)/256 - LOG2_HI_PART. */
+#define LOG2_BY_256_LO_PART \
+ 0.000000000000745396456746323365681353781544922399845L
+
+long double
+expl (long double x)
{
- 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));
+ if (isnanl (x))
+ return x;
+
+ if (x >= (long double) LDBL_MAX_EXP * LOG2_PLUS_EPSILON)
+ /* x > LDBL_MAX_EXP * log(2)
+ hence exp(x) > 2^LDBL_MAX_EXP, overflows to Infinity. */
+ return HUGE_VALL;
+
+ if (x <= (long double) (LDBL_MIN_EXP - 1 - LDBL_MANT_DIG) * LOG2_PLUS_EPSILON)
+ /* x < (LDBL_MIN_EXP - 1 - LDBL_MANT_DIG) * log(2)
+ hence exp(x) < 2^(LDBL_MIN_EXP-1-LDBL_MANT_DIG),
+ underflows to zero. */
+ return 0.0L;
+
+ /* Decompose x into
+ x = n * log(2) + m * log(2)/256 + y
+ where
+ n is an integer,
+ m is an integer, -128 <= m <= 128,
+ y is a number, |y| <= log(2)/512 + epsilon = 0.00135...
+ Then
+ exp(x) = 2^n * exp(m * log(2)/256) * exp(y)
+ The first factor is an ldexpl() call.
+ The second factor is a table lookup.
+ The third factor is computed
+ - either as sinh(y) + cosh(y)
+ where sinh(y) is computed through the power series:
+ sinh(y) = y + y^3/3! + y^5/5! + ...
+ and cosh(y) is computed as hypot(1, sinh(y)),
+ - or as exp(2*z) = (1 + tanh(z)) / (1 - tanh(z))
+ where z = y/2
+ and tanh(z) is computed through its power series:
+ tanh(z) = z
+ - 1/3 * z^3
+ + 2/15 * z^5
+ - 17/315 * z^7
+ + 62/2835 * z^9
+ - 1382/155925 * z^11
+ + 21844/6081075 * z^13
+ - 929569/638512875 * z^15
+ + ...
+ Since |z| <= log(2)/1024 < 0.0007, the relative contribution of the
+ z^13 term is < 0.0007^12 < 2^-120 <= 2^-LDBL_MANT_DIG, therefore we
+ can truncate the series after the z^11 term.
+
+ Given the usual bounds LDBL_MAX_EXP <= 16384, LDBL_MIN_EXP >= -16381,
+ LDBL_MANT_DIG <= 120, we can estimate x: -11440 <= x <= 11357.
+ This means, when dividing x by log(2), where we want x mod log(2)
+ to be precise to LDBL_MANT_DIG bits, we have to use an approximation
+ to log(2) that has 14+LDBL_MANT_DIG bits. */
+
+ {
+ long double nm = roundl (x * LOG2_BY_256_INVERSE); /* = 256 * n + m */
+ /* n has at most 15 bits, nm therefore has at most 23 bits, therefore
+ n * LOG2_HI_PART is computed exactly, and n * LOG2_LO_PART is computed
+ with an absolute error < 2^15 * 2e-10 * 2^-LDBL_MANT_DIG. */
+ long double y_tmp = x - nm * LOG2_BY_256_HI_PART;
+ long double y = y_tmp - nm * LOG2_BY_256_LO_PART;
+ long double z = 0.5L * y;
+
+/* Coefficients of the power series for tanh(z). */
+#define TANH_COEFF_1 1.0L
+#define TANH_COEFF_3 -0.333333333333333333333333333333333333334L
+#define TANH_COEFF_5 0.133333333333333333333333333333333333334L
+#define TANH_COEFF_7 -0.053968253968253968253968253968253968254L
+#define TANH_COEFF_9 0.0218694885361552028218694885361552028218L
+#define TANH_COEFF_11 -0.00886323552990219656886323552990219656886L
+#define TANH_COEFF_13 0.00359212803657248101692546136990581435026L
+#define TANH_COEFF_15 -0.00145583438705131826824948518070211191904L
+
+ long double z2 = z * z;
+ long double tanh_z =
+ (((((TANH_COEFF_11
+ * z2 + TANH_COEFF_9)
+ * z2 + TANH_COEFF_7)
+ * z2 + TANH_COEFF_5)
+ * z2 + TANH_COEFF_3)
+ * z2 + TANH_COEFF_1)
+ * z;
+
+ long double exp_y = (1.0L + tanh_z) / (1.0L - tanh_z);
+
+ int n = (int) roundl (nm * (1.0L / 256.0L));
+ int m = (int) nm - 256 * n;
+
+ return ldexpl (gl_expl_table[128 + m] * exp_y, n);
+ }
}
+
#endif
projects / gnulib.git / blobdiff