1 | /* origin: FreeBSD /usr/src/lib/msun/src/s_expm1f.c */ |
2 | /* |
3 | * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. |
4 | */ |
5 | /* |
6 | * ==================================================== |
7 | * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. |
8 | * |
9 | * Developed at SunPro, a Sun Microsystems, Inc. business. |
10 | * Permission to use, copy, modify, and distribute this |
11 | * software is freely granted, provided that this notice |
12 | * is preserved. |
13 | * ==================================================== |
14 | */ |
15 | |
16 | const O_THRESHOLD: f32 = 8.8721679688e+01; /* 0x42b17180 */ |
17 | const LN2_HI: f32 = 6.9313812256e-01; /* 0x3f317180 */ |
18 | const LN2_LO: f32 = 9.0580006145e-06; /* 0x3717f7d1 */ |
19 | const INV_LN2: f32 = 1.4426950216e+00; /* 0x3fb8aa3b */ |
20 | /* |
21 | * Domain [-0.34568, 0.34568], range ~[-6.694e-10, 6.696e-10]: |
22 | * |6 / x * (1 + 2 * (1 / (exp(x) - 1) - 1 / x)) - q(x)| < 2**-30.04 |
23 | * Scaled coefficients: Qn_here = 2**n * Qn_for_q (see s_expm1.c): |
24 | */ |
25 | const Q1: f32 = -3.3333212137e-2; /* -0x888868.0p-28 */ |
26 | const Q2: f32 = 1.5807170421e-3; /* 0xcf3010.0p-33 */ |
27 | |
28 | /// Exponential, base *e*, of x-1 (f32) |
29 | /// |
30 | /// Calculates the exponential of `x` and subtract 1, that is, *e* raised |
31 | /// to the power `x` minus 1 (where *e* is the base of the natural |
32 | /// system of logarithms, approximately 2.71828). |
33 | /// The result is accurate even for small values of `x`, |
34 | /// where using `exp(x)-1` would lose many significant digits. |
35 | #[cfg_attr (all(test, assert_no_panic), no_panic::no_panic)] |
36 | pub fn expm1f(mut x: f32) -> f32 { |
37 | let x1p127 = f32::from_bits(0x7f000000); // 0x1p127f === 2 ^ 127 |
38 | |
39 | let mut hx = x.to_bits(); |
40 | let sign = (hx >> 31) != 0; |
41 | hx &= 0x7fffffff; |
42 | |
43 | /* filter out huge and non-finite argument */ |
44 | if hx >= 0x4195b844 { |
45 | /* if |x|>=27*ln2 */ |
46 | if hx > 0x7f800000 { |
47 | /* NaN */ |
48 | return x; |
49 | } |
50 | if sign { |
51 | return -1.; |
52 | } |
53 | if x > O_THRESHOLD { |
54 | x *= x1p127; |
55 | return x; |
56 | } |
57 | } |
58 | |
59 | let k: i32; |
60 | let hi: f32; |
61 | let lo: f32; |
62 | let mut c = 0f32; |
63 | /* argument reduction */ |
64 | if hx > 0x3eb17218 { |
65 | /* if |x| > 0.5 ln2 */ |
66 | if hx < 0x3F851592 { |
67 | /* and |x| < 1.5 ln2 */ |
68 | if !sign { |
69 | hi = x - LN2_HI; |
70 | lo = LN2_LO; |
71 | k = 1; |
72 | } else { |
73 | hi = x + LN2_HI; |
74 | lo = -LN2_LO; |
75 | k = -1; |
76 | } |
77 | } else { |
78 | k = (INV_LN2 * x + (if sign { -0.5 } else { 0.5 })) as i32; |
79 | let t = k as f32; |
80 | hi = x - t * LN2_HI; /* t*ln2_hi is exact here */ |
81 | lo = t * LN2_LO; |
82 | } |
83 | x = hi - lo; |
84 | c = (hi - x) - lo; |
85 | } else if hx < 0x33000000 { |
86 | /* when |x|<2**-25, return x */ |
87 | if hx < 0x00800000 { |
88 | force_eval!(x * x); |
89 | } |
90 | return x; |
91 | } else { |
92 | k = 0; |
93 | } |
94 | |
95 | /* x is now in primary range */ |
96 | let hfx = 0.5 * x; |
97 | let hxs = x * hfx; |
98 | let r1 = 1. + hxs * (Q1 + hxs * Q2); |
99 | let t = 3. - r1 * hfx; |
100 | let mut e = hxs * ((r1 - t) / (6. - x * t)); |
101 | if k == 0 { |
102 | /* c is 0 */ |
103 | return x - (x * e - hxs); |
104 | } |
105 | e = x * (e - c) - c; |
106 | e -= hxs; |
107 | /* exp(x) ~ 2^k (x_reduced - e + 1) */ |
108 | if k == -1 { |
109 | return 0.5 * (x - e) - 0.5; |
110 | } |
111 | if k == 1 { |
112 | if x < -0.25 { |
113 | return -2. * (e - (x + 0.5)); |
114 | } |
115 | return 1. + 2. * (x - e); |
116 | } |
117 | let twopk = f32::from_bits(((0x7f + k) << 23) as u32); /* 2^k */ |
118 | if (k < 0) || (k > 56) { |
119 | /* suffice to return exp(x)-1 */ |
120 | let mut y = x - e + 1.; |
121 | if k == 128 { |
122 | y = y * 2. * x1p127; |
123 | } else { |
124 | y = y * twopk; |
125 | } |
126 | return y - 1.; |
127 | } |
128 | let uf = f32::from_bits(((0x7f - k) << 23) as u32); /* 2^-k */ |
129 | if k < 23 { |
130 | (x - e + (1. - uf)) * twopk |
131 | } else { |
132 | (x - (e + uf) + 1.) * twopk |
133 | } |
134 | } |
135 | |