pxfm/exponents/
exp2.rs

1/*
2 * // Copyright (c) Radzivon Bartoshyk 7/2025. All rights reserved.
3 * //
4 * // Redistribution and use in source and binary forms, with or without modification,
5 * // are permitted provided that the following conditions are met:
6 * //
7 * // 1.  Redistributions of source code must retain the above copyright notice, this
8 * // list of conditions and the following disclaimer.
9 * //
10 * // 2.  Redistributions in binary form must reproduce the above copyright notice,
11 * // this list of conditions and the following disclaimer in the documentation
12 * // and/or other materials provided with the distribution.
13 * //
14 * // 3.  Neither the name of the copyright holder nor the names of its
15 * // contributors may be used to endorse or promote products derived from
16 * // this software without specific prior written permission.
17 * //
18 * // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
19 * // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
20 * // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
21 * // DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
22 * // FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
23 * // DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
24 * // SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
25 * // CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
26 * // OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
27 * // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
28 */
29use crate::common::{dd_fmla, f_fmla};
30use crate::double_double::DoubleDouble;
31use crate::exponents::auxiliary::fast_ldexp;
32use crate::exponents::exp::{EXP_REDUCE_T0, EXP_REDUCE_T1, to_denormal};
33use crate::round_ties_even::RoundTiesEven;
34
35#[inline]
36fn exp2_poly_dd(z: f64) -> DoubleDouble {
37    const C: [(u64, u64); 6] = [
38        (0x3bbabc9e3b39873e, 0x3f262e42fefa39ef),
39        (0xbae5e43a53e44950, 0x3e4ebfbdff82c58f),
40        (0xba0d3a15710d3d83, 0x3d6c6b08d704a0c0),
41        (0x3914dd5d2a5e025a, 0x3c83b2ab6fba4e77),
42        (0xb83dc47e47beb9dd, 0x3b95d87fe7a66459),
43        (0xb744fcd51fcb7640, 0x3aa430912f9fb79d),
44    ];
45
46    let mut r = DoubleDouble::quick_mul_f64_add(
47        DoubleDouble::from_bit_pair(C[5]),
48        z,
49        DoubleDouble::from_bit_pair(C[4]),
50    );
51    r = DoubleDouble::quick_mul_f64_add(r, z, DoubleDouble::from_bit_pair(C[3]));
52    r = DoubleDouble::quick_mul_f64_add(r, z, DoubleDouble::from_bit_pair(C[2]));
53    r = DoubleDouble::quick_mul_f64_add(r, z, DoubleDouble::from_bit_pair(C[1]));
54    DoubleDouble::quick_mul_f64_add(r, z, DoubleDouble::from_bit_pair(C[0]))
55}
56
57#[cold]
58fn exp2_accurate(x: f64) -> f64 {
59    let mut ix = x.to_bits();
60    let sx = 4096.0 * x;
61    let fx = sx.round_ties_even_finite();
62    let z = sx - fx;
63    let k: i64 = unsafe {
64        fx.to_int_unchecked::<i64>() // this is already finite here
65    };
66    let i1 = k & 0x3f;
67    let i0 = (k >> 6) & 0x3f;
68    let ie = k >> 12;
69
70    let t0 = DoubleDouble::from_bit_pair(EXP_REDUCE_T0[i0 as usize]);
71    let t1 = DoubleDouble::from_bit_pair(EXP_REDUCE_T1[i1 as usize]);
72    let dt = DoubleDouble::quick_mult(t0, t1);
73
74    let mut f = exp2_poly_dd(z);
75    f = DoubleDouble::quick_mult_f64(f, z);
76    if ix <= 0xc08ff00000000000u64 {
77        // x >= -1022
78        // for -0x1.71547652b82fep-54 <= x <= 0x1.71547652b82fdp-53,
79        // exp2(x) round to x to nearest
80        if f64::from_bits(0xbc971547652b82fe) <= x && x <= f64::from_bits(0x3ca71547652b82fd) {
81            return dd_fmla(x, 0.5, 1.0);
82        } else if (k & 0xfff) == 0 {
83            // 4096*x rounds to 4096*integer
84            let zf = DoubleDouble::from_exact_add(dt.hi, f.hi);
85            let zfl = DoubleDouble::from_exact_add(zf.lo, f.lo);
86            f.hi = zf.hi;
87            f.lo = zfl.hi;
88            ix = zfl.hi.to_bits();
89            if ix & 0x000fffffffffffff == 0 {
90                // fl is a power of 2
91                if ((ix >> 52) & 0x7ff) != 0 {
92                    // |fl| is Inf
93                    let v = zfl.lo.to_bits();
94                    let d: i64 = ((((ix as i64) >> 63) ^ ((v as i64) >> 63)) as u64)
95                        .wrapping_shl(1)
96                        .wrapping_add(1) as i64;
97                    ix = ix.wrapping_add(d as u64);
98                    f.lo = f64::from_bits(ix);
99                }
100            }
101        } else {
102            f = DoubleDouble::quick_mult(f, dt);
103            f = DoubleDouble::add(dt, f);
104        }
105        let hf = DoubleDouble::from_exact_add(f.hi, f.lo);
106
107        fast_ldexp(hf.hi, ie as i32)
108    } else {
109        ix = 1u64.wrapping_sub(ie as u64).wrapping_shl(52);
110        f = DoubleDouble::quick_mult(f, dt);
111        f = DoubleDouble::add(dt, f);
112        let zve = DoubleDouble::from_exact_add(f64::from_bits(ix), f.hi);
113        f.hi = zve.hi;
114        f.lo += zve.lo;
115
116        to_denormal(f.to_f64())
117    }
118}
119
120/// Computes exp2
121///
122/// Max found ULP 0.5
123pub fn f_exp2(x: f64) -> f64 {
124    let mut ix = x.to_bits();
125    let ax = ix.wrapping_shl(1);
126    if ax == 0 {
127        return 1.0;
128    }
129    if ax >= 0x8120000000000000u64 {
130        // |x| >= 1024
131        if ax > 0xffe0000000000000u64 {
132            return x + x; // nan
133        }
134        if ax == 0xffe0000000000000u64 {
135            return if (ix >> 63) != 0 { 0.0 } else { x };
136        }
137        // +/-inf
138        if (ix >> 63) != 0 {
139            // x <= -1024
140            if ix >= 0xc090cc0000000000u64 {
141                // x <= -1075
142                const Z: f64 = f64::from_bits(0x0010000000000000);
143                return Z * Z;
144            }
145        } else {
146            // x >= 1024
147            return f64::from_bits(0x7fe0000000000000) * x;
148        }
149    }
150
151    // for |x| <= 0x1.71547652b82fep-54, 2^x rounds to 1 to nearest
152    // this avoids a spurious underflow in z*z below
153    if ax <= 0x792e2a8eca5705fcu64 {
154        return 1.0 + f64::copysign(f64::from_bits(0x3c90000000000000), x);
155    }
156
157    let m = ix.wrapping_shl(12);
158    let ex = (ax >> 53).wrapping_sub(0x3ff);
159    let frac = ex >> 63 | m << (ex & 63);
160    let sx = 4096.0 * x;
161    let fx = sx.round_ties_even_finite();
162    let z = sx - fx;
163    let z2 = z * z;
164    let k = unsafe {
165        fx.to_int_unchecked::<i64>() // this already finite here
166    };
167    let i1 = k & 0x3f;
168    let i0 = (k >> 6) & 0x3f;
169    let ie = k >> 12;
170    let t0 = DoubleDouble::from_bit_pair(EXP_REDUCE_T0[i0 as usize]);
171    let t1 = DoubleDouble::from_bit_pair(EXP_REDUCE_T1[i1 as usize]);
172    let ti0 = DoubleDouble::quick_mult(t0, t1);
173    const C: [u64; 4] = [
174        0x3f262e42fefa39ef,
175        0x3e4ebfbdff82c58f,
176        0x3d6c6b08d73b3e01,
177        0x3c83b2ab6fdda001,
178    ];
179    let tz = ti0.hi * z;
180    let mut fh = ti0.hi;
181
182    let p0 = f_fmla(z, f64::from_bits(C[1]), f64::from_bits(C[0]));
183    let p1 = f_fmla(z, f64::from_bits(C[3]), f64::from_bits(C[2]));
184    let p2 = f_fmla(z2, p1, p0);
185
186    let mut fl = f_fmla(tz, p2, ti0.lo);
187
188    const EPS: f64 = f64::from_bits(0x3c0833beace2b6fe);
189
190    if ix <= 0xc08ff00000000000u64 {
191        // x >= -1022
192        if frac != 0 {
193            let ub = fh + (fl + EPS);
194            fh += fl - EPS;
195            if ub != fh {
196                return exp2_accurate(x);
197            }
198        }
199        fh = fast_ldexp(fh, ie as i32);
200    } else {
201        // subnormal case
202        ix = 1u64.wrapping_sub(ie as u64).wrapping_shl(52);
203        let rs = DoubleDouble::from_exact_add(f64::from_bits(ix), fh);
204        fl += rs.lo;
205        fh = rs.hi;
206        if frac != 0 {
207            let ub = fh + (fl + EPS);
208            fh += fl - EPS;
209            if ub != fh {
210                return exp2_accurate(x);
211            }
212        }
213        // when 2^x is exact, no underflow should be raised
214        fh = to_denormal(fh);
215    }
216    fh
217}
218
219#[cfg(test)]
220mod tests {
221    use super::*;
222
223    #[test]
224    fn test_exp2d() {
225        assert_eq!(f_exp2(2.0), 4.0);
226        assert_eq!(f_exp2(3.0), 8.0);
227        assert_eq!(f_exp2(4.0), 16.0);
228        assert_eq!(f_exp2(0.35f64), 1.2745606273192622);
229        assert_eq!(f_exp2(-0.6f64), 0.6597539553864471);
230        assert_eq!(f_exp2(f64::INFINITY), f64::INFINITY);
231        assert_eq!(f_exp2(f64::NEG_INFINITY), 0.);
232        assert!(f_exp2(f64::NAN).is_nan());
233    }
234}
pxfm/exponents/exp2.rs

pxfm/exponents/
exp2.rs
