nsz Git - libm/blob - src/math/expl.c

   1 /* origin: OpenBSD /usr/src/lib/libm/src/ld80/e_expl.c */
   2 /*
   3  * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
   4  *
   5  * Permission to use, copy, modify, and distribute this software for any
   6  * purpose with or without fee is hereby granted, provided that the above
   7  * copyright notice and this permission notice appear in all copies.
   8  *
   9  * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
  10  * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
  11  * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
  12  * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
  13  * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
  14  * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
  15  * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
  16  */
  17 /*
  18  *      Exponential function, long double precision
  19  *
  20  *
  21  * SYNOPSIS:
  22  *
  23  * long double x, y, expl();
  24  *
  25  * y = expl( x );
  26  *
  27  *
  28  * DESCRIPTION:
  29  *
  30  * Returns e (2.71828...) raised to the x power.
  31  *
  32  * Range reduction is accomplished by separating the argument
  33  * into an integer k and fraction f such that
  34  *
  35  *     x    k  f
  36  *    e  = 2  e.
  37  *
  38  * A Pade' form of degree 2/3 is used to approximate exp(f) - 1
  39  * in the basic range [-0.5 ln 2, 0.5 ln 2].
  40  *
  41  *
  42  * ACCURACY:
  43  *
  44  *                      Relative error:
  45  * arithmetic   domain     # trials      peak         rms
  46  *    IEEE      +-10000     50000       1.12e-19    2.81e-20
  47  *
  48  *
  49  * Error amplification in the exponential function can be
  50  * a serious matter.  The error propagation involves
  51  * exp( X(1+delta) ) = exp(X) ( 1 + X*delta + ... ),
  52  * which shows that a 1 lsb error in representing X produces
  53  * a relative error of X times 1 lsb in the function.
  54  * While the routine gives an accurate result for arguments
  55  * that are exactly represented by a long double precision
  56  * computer number, the result contains amplified roundoff
  57  * error for large arguments not exactly represented.
  58  *
  59  *
  60  * ERROR MESSAGES:
  61  *
  62  *   message         condition      value returned
  63  * exp underflow    x < MINLOG         0.0
  64  * exp overflow     x > MAXLOG         MAXNUM
  65  *
  66  */
  67
  68 #include "libm.h"
  69
  70 #if LD64
  71 long double expl(long double x)
  72 {
  73         return x;
  74 }
  75 #elif LD80
  76
  77 static long double P[3] = {
  78  1.2617719307481059087798E-4L,
  79  3.0299440770744196129956E-2L,
  80  9.9999999999999999991025E-1L,
  81 };
  82 static long double Q[4] = {
  83  3.0019850513866445504159E-6L,
  84  2.5244834034968410419224E-3L,
  85  2.2726554820815502876593E-1L,
  86  2.0000000000000000000897E0L,
  87 };
  88 static const long double
  89 C1 = 6.9314575195312500000000E-1L,
  90 C2 = 1.4286068203094172321215E-6L,
  91 MAXLOGL = 1.1356523406294143949492E4L,
  92 MINLOGL = -1.13994985314888605586758E4L,
  93 LOG2EL = 1.4426950408889634073599E0L;
  94
  95 long double expl(long double x)
  96 {
  97         long double px, xx;
  98         int n;
  99
 100         if (isnan(x))
 101                 return x;
 102         if (x > MAXLOGL)
 103                 return INFINITY;
 104         if (x < MINLOGL)
 105                 return 0.0L;
 106
 107         /* Express e**x = e**g 2**n
 108          *   = e**g e**(n loge(2))
 109          *   = e**(g + n loge(2))
 110          */
 111         px = floorl(LOG2EL * x + 0.5L); /* floor() truncates toward -infinity. */
 112         n = px;
 113         x -= px * C1;
 114         x -= px * C2;
 115
 116         /* rational approximation for exponential
 117          * of the fractional part:
 118          * e**x =  1 + 2x P(x**2)/(Q(x**2) - P(x**2))
 119          */
 120         xx = x * x;
 121         px = x * __polevll(xx, P, 2);
 122         x =  px/(__polevll(xx, Q, 3) - px);
 123         x = 1.0L + ldexpl(x, 1);
 124         x = ldexpl(x, n);
 125         return x;
 126 }
 127 #endif