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 | use super::scalbn; |
69 | |
70 | const HALF: [f64; 2] = [0.5, -0.5]; |
71 | const LN2HI: f64 = 6.93147180369123816490e-01; /* 0x3fe62e42, 0xfee00000 */ |
72 | const LN2LO: f64 = 1.90821492927058770002e-10; /* 0x3dea39ef, 0x35793c76 */ |
73 | const INVLN2: f64 = 1.44269504088896338700e+00; /* 0x3ff71547, 0x652b82fe */ |
74 | const P1: f64 = 1.66666666666666019037e-01; /* 0x3FC55555, 0x5555553E */ |
75 | const P2: f64 = -2.77777777770155933842e-03; /* 0xBF66C16C, 0x16BEBD93 */ |
76 | const P3: f64 = 6.61375632143793436117e-05; /* 0x3F11566A, 0xAF25DE2C */ |
77 | const P4: f64 = -1.65339022054652515390e-06; /* 0xBEBBBD41, 0xC5D26BF1 */ |
78 | const P5: f64 = 4.13813679705723846039e-08; /* 0x3E663769, 0x72BEA4D0 */ |
79 | |
80 | /// Exponential, base *e* (f64) |
81 | /// |
82 | /// Calculate the exponential of `x`, that is, *e* raised to the power `x` |
83 | /// (where *e* is the base of the natural system of logarithms, approximately 2.71828). |
84 | #[cfg_attr (all(test, assert_no_panic), no_panic::no_panic)] |
85 | pub fn exp(mut x: f64) -> f64 { |
86 | let x1p1023 = f64::from_bits(0x7fe0000000000000); // 0x1p1023 === 2 ^ 1023 |
87 | let x1p_149 = f64::from_bits(0x36a0000000000000); // 0x1p-149 === 2 ^ -149 |
88 | |
89 | let hi: f64; |
90 | let lo: f64; |
91 | let c: f64; |
92 | let xx: f64; |
93 | let y: f64; |
94 | let k: i32; |
95 | let sign: i32; |
96 | let mut hx: u32; |
97 | |
98 | hx = (x.to_bits() >> 32) as u32; |
99 | sign = (hx >> 31) as i32; |
100 | hx &= 0x7fffffff; /* high word of |x| */ |
101 | |
102 | /* special cases */ |
103 | if hx >= 0x4086232b { |
104 | /* if |x| >= 708.39... */ |
105 | if x.is_nan() { |
106 | return x; |
107 | } |
108 | if x > 709.782712893383973096 { |
109 | /* overflow if x!=inf */ |
110 | x *= x1p1023; |
111 | return x; |
112 | } |
113 | if x < -708.39641853226410622 { |
114 | /* underflow if x!=-inf */ |
115 | force_eval!((-x1p_149 / x) as f32); |
116 | if x < -745.13321910194110842 { |
117 | return 0.; |
118 | } |
119 | } |
120 | } |
121 | |
122 | /* argument reduction */ |
123 | if hx > 0x3fd62e42 { |
124 | /* if |x| > 0.5 ln2 */ |
125 | if hx >= 0x3ff0a2b2 { |
126 | /* if |x| >= 1.5 ln2 */ |
127 | k = (INVLN2 * x + i!(HALF, sign as usize)) as i32; |
128 | } else { |
129 | k = 1 - sign - sign; |
130 | } |
131 | hi = x - k as f64 * LN2HI; /* k*ln2hi is exact here */ |
132 | lo = k as f64 * LN2LO; |
133 | x = hi - lo; |
134 | } else if hx > 0x3e300000 { |
135 | /* if |x| > 2**-28 */ |
136 | k = 0; |
137 | hi = x; |
138 | lo = 0.; |
139 | } else { |
140 | /* inexact if x!=0 */ |
141 | force_eval!(x1p1023 + x); |
142 | return 1. + x; |
143 | } |
144 | |
145 | /* x is now in primary range */ |
146 | xx = x * x; |
147 | c = x - xx * (P1 + xx * (P2 + xx * (P3 + xx * (P4 + xx * P5)))); |
148 | y = 1. + (x * c / (2. - c) - lo + hi); |
149 | if k == 0 { |
150 | y |
151 | } else { |
152 | scalbn(y, k) |
153 | } |
154 | } |
155 | |