1 | /* origin: FreeBSD /usr/src/lib/msun/src/s_expm1.c */ |
2 | /* |
3 | * ==================================================== |
4 | * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. |
5 | * |
6 | * Developed at SunPro, a Sun Microsystems, Inc. business. |
7 | * Permission to use, copy, modify, and distribute this |
8 | * software is freely granted, provided that this notice |
9 | * is preserved. |
10 | * ==================================================== |
11 | */ |
12 | |
13 | use core::f64; |
14 | |
15 | const O_THRESHOLD: f64 = 7.09782712893383973096e+02; /* 0x40862E42, 0xFEFA39EF */ |
16 | const LN2_HI: f64 = 6.93147180369123816490e-01; /* 0x3fe62e42, 0xfee00000 */ |
17 | const LN2_LO: f64 = 1.90821492927058770002e-10; /* 0x3dea39ef, 0x35793c76 */ |
18 | const INVLN2: f64 = 1.44269504088896338700e+00; /* 0x3ff71547, 0x652b82fe */ |
19 | /* Scaled Q's: Qn_here = 2**n * Qn_above, for R(2*z) where z = hxs = x*x/2: */ |
20 | const Q1: f64 = -3.33333333333331316428e-02; /* BFA11111 111110F4 */ |
21 | const Q2: f64 = 1.58730158725481460165e-03; /* 3F5A01A0 19FE5585 */ |
22 | const Q3: f64 = -7.93650757867487942473e-05; /* BF14CE19 9EAADBB7 */ |
23 | const Q4: f64 = 4.00821782732936239552e-06; /* 3ED0CFCA 86E65239 */ |
24 | const Q5: f64 = -2.01099218183624371326e-07; /* BE8AFDB7 6E09C32D */ |
25 | |
26 | /// Exponential, base *e*, of x-1 (f64) |
27 | /// |
28 | /// Calculates the exponential of `x` and subtract 1, that is, *e* raised |
29 | /// to the power `x` minus 1 (where *e* is the base of the natural |
30 | /// system of logarithms, approximately 2.71828). |
31 | /// The result is accurate even for small values of `x`, |
32 | /// where using `exp(x)-1` would lose many significant digits. |
33 | #[cfg_attr (all(test, assert_no_panic), no_panic::no_panic)] |
34 | pub fn expm1(mut x: f64) -> f64 { |
35 | let hi: f64; |
36 | let lo: f64; |
37 | let k: i32; |
38 | let c: f64; |
39 | let mut t: f64; |
40 | let mut y: f64; |
41 | |
42 | let mut ui = x.to_bits(); |
43 | let hx = ((ui >> 32) & 0x7fffffff) as u32; |
44 | let sign = (ui >> 63) as i32; |
45 | |
46 | /* filter out huge and non-finite argument */ |
47 | if hx >= 0x4043687A { |
48 | /* if |x|>=56*ln2 */ |
49 | if x.is_nan() { |
50 | return x; |
51 | } |
52 | if sign != 0 { |
53 | return -1.0; |
54 | } |
55 | if x > O_THRESHOLD { |
56 | x *= f64::from_bits(0x7fe0000000000000); |
57 | return x; |
58 | } |
59 | } |
60 | |
61 | /* argument reduction */ |
62 | if hx > 0x3fd62e42 { |
63 | /* if |x| > 0.5 ln2 */ |
64 | if hx < 0x3FF0A2B2 { |
65 | /* and |x| < 1.5 ln2 */ |
66 | if sign == 0 { |
67 | hi = x - LN2_HI; |
68 | lo = LN2_LO; |
69 | k = 1; |
70 | } else { |
71 | hi = x + LN2_HI; |
72 | lo = -LN2_LO; |
73 | k = -1; |
74 | } |
75 | } else { |
76 | k = (INVLN2 * x + if sign != 0 { -0.5 } else { 0.5 }) as i32; |
77 | t = k as f64; |
78 | hi = x - t * LN2_HI; /* t*ln2_hi is exact here */ |
79 | lo = t * LN2_LO; |
80 | } |
81 | x = hi - lo; |
82 | c = (hi - x) - lo; |
83 | } else if hx < 0x3c900000 { |
84 | /* |x| < 2**-54, return x */ |
85 | if hx < 0x00100000 { |
86 | force_eval!(x); |
87 | } |
88 | return x; |
89 | } else { |
90 | c = 0.0; |
91 | k = 0; |
92 | } |
93 | |
94 | /* x is now in primary range */ |
95 | let hfx = 0.5 * x; |
96 | let hxs = x * hfx; |
97 | let r1 = 1.0 + hxs * (Q1 + hxs * (Q2 + hxs * (Q3 + hxs * (Q4 + hxs * Q5)))); |
98 | t = 3.0 - r1 * hfx; |
99 | let mut e = hxs * ((r1 - t) / (6.0 - x * t)); |
100 | if k == 0 { |
101 | /* c is 0 */ |
102 | return x - (x * e - hxs); |
103 | } |
104 | e = x * (e - c) - c; |
105 | e -= hxs; |
106 | /* exp(x) ~ 2^k (x_reduced - e + 1) */ |
107 | if k == -1 { |
108 | return 0.5 * (x - e) - 0.5; |
109 | } |
110 | if k == 1 { |
111 | if x < -0.25 { |
112 | return -2.0 * (e - (x + 0.5)); |
113 | } |
114 | return 1.0 + 2.0 * (x - e); |
115 | } |
116 | ui = ((0x3ff + k) as u64) << 52; /* 2^k */ |
117 | let twopk = f64::from_bits(ui); |
118 | if k < 0 || k > 56 { |
119 | /* suffice to return exp(x)-1 */ |
120 | y = x - e + 1.0; |
121 | if k == 1024 { |
122 | y = y * 2.0 * f64::from_bits(0x7fe0000000000000); |
123 | } else { |
124 | y = y * twopk; |
125 | } |
126 | return y - 1.0; |
127 | } |
128 | ui = ((0x3ff - k) as u64) << 52; /* 2^-k */ |
129 | let uf = f64::from_bits(ui); |
130 | if k < 20 { |
131 | y = (x - e + (1.0 - uf)) * twopk; |
132 | } else { |
133 | y = (x - (e + uf) + 1.0) * twopk; |
134 | } |
135 | y |
136 | } |
137 | |
138 | #[cfg (test)] |
139 | mod tests { |
140 | #[test ] |
141 | fn sanity_check() { |
142 | assert_eq!(super::expm1(1.1), 2.0041660239464334); |
143 | } |
144 | } |
145 | |