src/math/exp.c

   1 /* origin: FreeBSD /usr/src/lib/msun/src/e_exp.c */
   2 /*
   3  * ====================================================
   4  * Copyright (C) 2004 by Sun Microsystems, Inc. All rights reserved.
   5  *
   6  * Permission to use, copy, modify, and distribute this
   7  * software is freely granted, provided that this notice
   8  * is preserved.
   9  * ====================================================
  10  */
  11 /* exp(x)
  12  * Returns the exponential of x.
  13  *
  14  * Method
  15  *   1. Argument reduction:
  16  *      Reduce x to an r so that |r| <= 0.5*ln2 ~ 0.34658.
  17  *      Given x, find r and integer k such that
  18  *
  19  *               x = k*ln2 + r,  |r| <= 0.5*ln2.
  20  *
  21  *      Here r will be represented as r = hi-lo for better
  22  *      accuracy.
  23  *
  24  *   2. Approximation of exp(r) by a special rational function on
  25  *      the interval [0,0.34658]:
  26  *      Write
  27  *          R(r**2) = r*(exp(r)+1)/(exp(r)-1) = 2 + r*r/6 - r**4/360 + ...
  28  *      We use a special Remez algorithm on [0,0.34658] to generate
  29  *      a polynomial of degree 5 to approximate R. The maximum error
  30  *      of this polynomial approximation is bounded by 2**-59. In
  31  *      other words,
  32  *          R(z) ~ 2.0 + P1*z + P2*z**2 + P3*z**3 + P4*z**4 + P5*z**5
  33  *      (where z=r*r, and the values of P1 to P5 are listed below)
  34  *      and
  35  *          |                  5          |     -59
  36  *          | 2.0+P1*z+...+P5*z   -  R(z) | <= 2
  37  *          |                             |
  38  *      The computation of exp(r) thus becomes
  39  *                              2*r
  40  *              exp(r) = 1 + ----------
  41  *                            R(r) - r
  42  *                                 r*c(r)
  43  *                     = 1 + r + ----------- (for better accuracy)
  44  *                                2 - c(r)
  45  *      where
  46  *                              2       4             10
  47  *              c(r) = r - (P1*r  + P2*r  + ... + P5*r   ).
  48  *
  49  *   3. Scale back to obtain exp(x):
  50  *      From step 1, we have
  51  *         exp(x) = 2^k * exp(r)
  52  *
  53  * Special cases:
  54  *      exp(INF) is INF, exp(NaN) is NaN;
  55  *      exp(-INF) is 0, and
  56  *      for finite argument, only exp(0)=1 is exact.
  57  *
  58  * Accuracy:
  59  *      according to an error analysis, the error is always less than
  60  *      1 ulp (unit in the last place).
  61  *
  62  * Misc. info.
  63  *      For IEEE double
  64  *          if x >  709.782712893383973096 then exp(x) overflows
  65  *          if x < -745.133219101941108420 then exp(x) underflows
  66  */
  67
  68 #include "libm.h"
  69
  70 static const double
  71 half[2] = {0.5,-0.5},
  72 ln2hi = 6.93147180369123816490e-01, /* 0x3fe62e42, 0xfee00000 */
  73 ln2lo = 1.90821492927058770002e-10, /* 0x3dea39ef, 0x35793c76 */
  74 invln2 = 1.44269504088896338700e+00, /* 0x3ff71547, 0x652b82fe */
  75 P1   =  1.66666666666666019037e-01, /* 0x3FC55555, 0x5555553E */
  76 P2   = -2.77777777770155933842e-03, /* 0xBF66C16C, 0x16BEBD93 */
  77 P3   =  6.61375632143793436117e-05, /* 0x3F11566A, 0xAF25DE2C */
  78 P4   = -1.65339022054652515390e-06, /* 0xBEBBBD41, 0xC5D26BF1 */
  79 P5   =  4.13813679705723846039e-08; /* 0x3E663769, 0x72BEA4D0 */
  80
  81 double exp(double x)
  82 {
  83         double_t hi, lo, c, xx, y;
  84         int k, sign;
  85         uint32_t hx;
  86
  87         GET_HIGH_WORD(hx, x);
  88         sign = hx>>31;
  89         hx &= 0x7fffffff;  /* high word of |x| */
  90
  91         /* special cases */
  92         if (hx >= 0x4086232b) {  /* if |x| >= 708.39... */
  93                 if (isnan(x))
  94                         return x;
  95                 if (x > 709.782712893383973096) {
  96                         /* overflow if x!=inf */
  97                         x *= 0x1p1023;
  98                         return x;
  99                 }
 100                 if (x < -708.39641853226410622) {
 101                         /* underflow if x!=-inf */
 102                         FORCE_EVAL((float)(-0x1p-149/x));
 103                         if (x < -745.13321910194110842)
 104                                 return 0;
 105                 }
 106         }
 107
 108         /* argument reduction */
 109         if (hx > 0x3fd62e42) {  /* if |x| > 0.5 ln2 */
 110                 if (hx >= 0x3ff0a2b2)  /* if |x| >= 1.5 ln2 */
 111                         k = (int)(invln2*x + half[sign]);
 112                 else
 113                         k = 1 - sign - sign;
 114                 hi = x - k*ln2hi;  /* k*ln2hi is exact here */
 115                 lo = k*ln2lo;
 116                 x = hi - lo;
 117         } else if (hx > 0x3e300000)  {  /* if |x| > 2**-28 */
 118                 k = 0;
 119                 hi = x;
 120                 lo = 0;
 121         } else {
 122                 /* inexact if x!=0 */
 123                 FORCE_EVAL(0x1p1023 + x);
 124                 return 1 + x;
 125         }
 126
 127         /* x is now in primary range */
 128         xx = x*x;
 129         c = x - xx*(P1+xx*(P2+xx*(P3+xx*(P4+xx*P5))));
 130         y = 1 + (x*c/(2-c) - lo + hi);
 131         if (k == 0)
 132                 return y;
 133         return scalbn(y, k);
 134 }