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 5/6 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 LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
  71 long double expl(long double x)
  72 {
  73         return exp(x);
  74 }
  75 #elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
  76
  77 static const long double P[3] = {
  78  1.2617719307481059087798E-4L,
  79  3.0299440770744196129956E-2L,
  80  9.9999999999999999991025E-1L,
  81 };
  82 static const 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 LN2HI = 6.9314575195312500000000E-1L,
  90 LN2LO = 1.4286068203094172321215E-6L,
  91 LOG2E = 1.4426950408889634073599E0L;
  92
  93 long double expl(long double x)
  94 {
  95         long double px, xx;
  96         int k;
  97
  98         if (isnan(x))
  99                 return x;
 100         if (x > 11356.5234062941439488L) /* x > ln(2^16384 - 0.5) */
 101                 return x * 0x1p16383L;
 102         if (x < -11399.4985314888605581L) /* x < ln(2^-16446) */
 103                 return 0x1p-10000L * 0x1p-10000L;
 104
 105         /* Express e**x = e**f 2**k
 106          *   = e**(f + k ln(2))
 107          */
 108         px = floorl(LOG2E * x + 0.5);
 109         k = px;
 110         x -= px * LN2HI;
 111         x -= px * LN2LO;
 112
 113         /* rational approximation of the fractional part:
 114          * e**x =  1 + 2x P(x**2)/(Q(x**2) - x P(x**2))
 115          */
 116         xx = x * x;
 117         px = x * __polevll(xx, P, 2);
 118         x = px/(__polevll(xx, Q, 3) - px);
 119         x = 1.0 + 2.0 * x;
 120         return scalbnl(x, k);
 121 }
 122 #endif