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