]>
Commit | Line | Data |
---|---|---|
8faf50e0 XL |
1 | // origin: FreeBSD /usr/src/lib/msun/src/s_exp2f.c |
2 | //- | |
3 | // Copyright (c) 2005 David Schultz <das@FreeBSD.ORG> | |
4 | // All rights reserved. | |
5 | // | |
6 | // Redistribution and use in source and binary forms, with or without | |
7 | // modification, are permitted provided that the following conditions | |
8 | // are met: | |
9 | // 1. Redistributions of source code must retain the above copyright | |
10 | // notice, this list of conditions and the following disclaimer. | |
11 | // 2. Redistributions in binary form must reproduce the above copyright | |
12 | // notice, this list of conditions and the following disclaimer in the | |
13 | // documentation and/or other materials provided with the distribution. | |
14 | // | |
15 | // THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND | |
16 | // ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE | |
17 | // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE | |
18 | // ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE | |
19 | // FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL | |
20 | // DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS | |
21 | // OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) | |
22 | // HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT | |
23 | // LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY | |
24 | // OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF | |
25 | // SUCH DAMAGE. | |
26 | ||
27 | const TBLSIZE: usize = 16; | |
28 | ||
29 | static EXP2FT: [u64; TBLSIZE] = [ | |
30 | 0x3fe6a09e667f3bcd, | |
31 | 0x3fe7a11473eb0187, | |
32 | 0x3fe8ace5422aa0db, | |
33 | 0x3fe9c49182a3f090, | |
34 | 0x3feae89f995ad3ad, | |
35 | 0x3fec199bdd85529c, | |
36 | 0x3fed5818dcfba487, | |
37 | 0x3feea4afa2a490da, | |
38 | 0x3ff0000000000000, | |
39 | 0x3ff0b5586cf9890f, | |
40 | 0x3ff172b83c7d517b, | |
41 | 0x3ff2387a6e756238, | |
42 | 0x3ff306fe0a31b715, | |
43 | 0x3ff3dea64c123422, | |
44 | 0x3ff4bfdad5362a27, | |
45 | 0x3ff5ab07dd485429, | |
46 | ]; | |
47 | ||
48 | // exp2f(x): compute the base 2 exponential of x | |
49 | // | |
50 | // Accuracy: Peak error < 0.501 ulp; location of peak: -0.030110927. | |
51 | // | |
52 | // Method: (equally-spaced tables) | |
53 | // | |
54 | // Reduce x: | |
55 | // x = k + y, for integer k and |y| <= 1/2. | |
56 | // Thus we have exp2f(x) = 2**k * exp2(y). | |
57 | // | |
58 | // Reduce y: | |
59 | // y = i/TBLSIZE + z for integer i near y * TBLSIZE. | |
60 | // Thus we have exp2(y) = exp2(i/TBLSIZE) * exp2(z), | |
61 | // with |z| <= 2**-(TBLSIZE+1). | |
62 | // | |
63 | // We compute exp2(i/TBLSIZE) via table lookup and exp2(z) via a | |
64 | // degree-4 minimax polynomial with maximum error under 1.4 * 2**-33. | |
65 | // Using double precision for everything except the reduction makes | |
66 | // roundoff error insignificant and simplifies the scaling step. | |
67 | // | |
68 | // This method is due to Tang, but I do not use his suggested parameters: | |
69 | // | |
70 | // Tang, P. Table-driven Implementation of the Exponential Function | |
71 | // in IEEE Floating-Point Arithmetic. TOMS 15(2), 144-157 (1989). | |
dc9dc135 XL |
72 | |
73 | /// Exponential, base 2 (f32) | |
74 | /// | |
75 | /// Calculate `2^x`, that is, 2 raised to the power `x`. | |
48663c56 | 76 | #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)] |
8faf50e0 XL |
77 | pub fn exp2f(mut x: f32) -> f32 { |
78 | let redux = f32::from_bits(0x4b400000) / TBLSIZE as f32; | |
79 | let p1 = f32::from_bits(0x3f317218); | |
80 | let p2 = f32::from_bits(0x3e75fdf0); | |
81 | let p3 = f32::from_bits(0x3d6359a4); | |
82 | let p4 = f32::from_bits(0x3c1d964e); | |
83 | ||
84 | // double_t t, r, z; | |
85 | // uint32_t ix, i0, k; | |
86 | ||
87 | let x1p127 = f32::from_bits(0x7f000000); | |
88 | ||
89 | /* Filter out exceptional cases. */ | |
90 | let ui = f32::to_bits(x); | |
91 | let ix = ui & 0x7fffffff; | |
92 | if ix > 0x42fc0000 { | |
93 | /* |x| > 126 */ | |
94 | if ix > 0x7f800000 { | |
95 | /* NaN */ | |
96 | return x; | |
97 | } | |
98 | if ui >= 0x43000000 && ui < 0x80000000 { | |
99 | /* x >= 128 */ | |
100 | x *= x1p127; | |
101 | return x; | |
102 | } | |
103 | if ui >= 0x80000000 { | |
104 | /* x < -126 */ | |
105 | if ui >= 0xc3160000 || (ui & 0x0000ffff != 0) { | |
106 | force_eval!(f32::from_bits(0x80000001) / x); | |
107 | } | |
108 | if ui >= 0xc3160000 { | |
109 | /* x <= -150 */ | |
110 | return 0.0; | |
111 | } | |
112 | } | |
113 | } else if ix <= 0x33000000 { | |
114 | /* |x| <= 0x1p-25 */ | |
115 | return 1.0 + x; | |
116 | } | |
117 | ||
118 | /* Reduce x, computing z, i0, and k. */ | |
119 | let ui = f32::to_bits(x + redux); | |
120 | let mut i0 = ui; | |
121 | i0 += TBLSIZE as u32 / 2; | |
122 | let k = i0 / TBLSIZE as u32; | |
123 | let ukf = f64::from_bits(((0x3ff + k) as u64) << 52); | |
124 | i0 &= TBLSIZE as u32 - 1; | |
125 | let mut uf = f32::from_bits(ui); | |
126 | uf -= redux; | |
127 | let z: f64 = (x - uf) as f64; | |
128 | /* Compute r = exp2(y) = exp2ft[i0] * p(z). */ | |
129 | let r: f64 = f64::from_bits(EXP2FT[i0 as usize]); | |
130 | let t: f64 = r as f64 * z; | |
131 | let r: f64 = r + t * (p1 as f64 + z * p2 as f64) + t * (z * z) * (p3 as f64 + z * p4 as f64); | |
132 | ||
133 | /* Scale by 2**k */ | |
134 | (r * ukf) as f32 | |
135 | } |