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). |
72 | |
73 | /// Exponential, base 2 (f32) |
74 | /// |
75 | /// Calculate `2^x`, that is, 2 raised to the power `x`. |
76 | #[cfg_attr (all(test, assert_no_panic), no_panic::no_panic)] |
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(i!(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 | } |
136 | |