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
70static const double
71half[2] = {0.5,-0.5},
72ln2hi = 6.93147180369123816490e-01, /* 0x3fe62e42, 0xfee00000 */
73ln2lo = 1.90821492927058770002e-10, /* 0x3dea39ef, 0x35793c76 */
74invln2 = 1.44269504088896338700e+00, /* 0x3ff71547, 0x652b82fe */
75P1   =  1.66666666666666019037e-01, /* 0x3FC55555, 0x5555553E */
76P2   = -2.77777777770155933842e-03, /* 0xBF66C16C, 0x16BEBD93 */
77P3   =  6.61375632143793436117e-05, /* 0x3F11566A, 0xAF25DE2C */
78P4   = -1.65339022054652515390e-06, /* 0xBEBBBD41, 0xC5D26BF1 */
79P5   =  4.13813679705723846039e-08; /* 0x3E663769, 0x72BEA4D0 */
80
81double 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}
135