pxfm/exponents/
exp10.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::double_double::DoubleDouble;
30use crate::exponents::auxiliary::fast_ldexp;
31use crate::exponents::exp::{EXP_REDUCE_T0, EXP_REDUCE_T1, to_denormal};
32use crate::exponents::expf::{ExpfBackend, GenericExpfBackend};
33use crate::rounding::CpuRoundTiesEven;
34use std::hint::black_box;
35
36#[inline]
37fn exp10_poly_dd(z: DoubleDouble) -> DoubleDouble {
38    const C: [(u64, u64); 6] = [
39        (0xbcaf48ad494ea102, 0x40026bb1bbb55516),
40        (0xbcae2bfab318d399, 0x40053524c73cea69),
41        (0x3ca81f50779e162b, 0x4000470591de2ca4),
42        (0x3c931a5cc5d3d313, 0x3ff2bd7609fd98c4),
43        (0x3c8910de8c68a0c2, 0x3fe1429ffd336aa3),
44        (0xbc605e703d496537, 0x3fca7ed7086882b4),
45    ];
46
47    let mut r = DoubleDouble::quick_mul_add(
48        DoubleDouble::from_bit_pair(C[5]),
49        z,
50        DoubleDouble::from_bit_pair(C[4]),
51    );
52    r = DoubleDouble::quick_mul_add(r, z, DoubleDouble::from_bit_pair(C[3]));
53    r = DoubleDouble::quick_mul_add(r, z, DoubleDouble::from_bit_pair(C[2]));
54    r = DoubleDouble::quick_mul_add(r, z, DoubleDouble::from_bit_pair(C[1]));
55    DoubleDouble::quick_mul_add(r, z, DoubleDouble::from_bit_pair(C[0]))
56}
57
58#[cold]
59fn as_exp10_accurate(x: f64) -> f64 {
60    let mut ix = x.to_bits();
61    let t = (f64::from_bits(0x40ca934f0979a371) * x).cpu_round_ties_even();
62    let jt: i64 = unsafe {
63        t.to_int_unchecked::<i64>() // t is already integer, this is just a conversion
64    };
65    let i1 = jt & 0x3f;
66    let i0 = (jt >> 6) & 0x3f;
67    let ie = jt >> 12;
68    let t0 = DoubleDouble::from_bit_pair(EXP_REDUCE_T0[i0 as usize]);
69    let t1 = DoubleDouble::from_bit_pair(EXP_REDUCE_T1[i1 as usize]);
70    let dt = DoubleDouble::quick_mult(t0, t1);
71
72    const L0: f64 = f64::from_bits(0x3f13441350800000);
73    const L1: f64 = f64::from_bits(0xbd1f79fef311f12b);
74    const L2: f64 = f64::from_bits(0xb9aac0b7c917826b);
75
76    let dx = x - L0 * t;
77    let dx_dd = DoubleDouble::quick_mult_f64(DoubleDouble::new(L2, L1), t);
78    let dz = DoubleDouble::full_add_f64(dx_dd, dx);
79    let mut f = exp10_poly_dd(dz);
80    f = DoubleDouble::quick_mult(dz, f);
81
82    let mut zfh: f64;
83
84    if ix < 0xc0733a7146f72a42u64 {
85        if (jt & 0xfff) == 0 {
86            f = DoubleDouble::from_exact_add(f.hi, f.lo);
87            let zt = DoubleDouble::from_exact_add(dt.hi, f.hi);
88            f.hi = zt.lo;
89            f = DoubleDouble::from_exact_add(f.hi, f.lo);
90            ix = f.hi.to_bits();
91            if (ix.wrapping_shl(12)) == 0 {
92                let l = f.lo.to_bits();
93                let sfh: i64 = ((ix as i64) >> 63) ^ ((l as i64) >> 63);
94                ix = ix.wrapping_add(((1i64 << 51) ^ sfh) as u64);
95            }
96            zfh = zt.hi + f64::from_bits(ix);
97        } else {
98            f = DoubleDouble::quick_mult(f, dt);
99            f = DoubleDouble::add(dt, f);
100            f = DoubleDouble::from_exact_add(f.hi, f.lo);
101            zfh = f.hi;
102        }
103        zfh = fast_ldexp(zfh, ie as i32);
104    } else {
105        ix = (1u64.wrapping_sub(ie as u64)) << 52;
106        f = DoubleDouble::quick_mult(f, dt);
107        f = DoubleDouble::add(dt, f);
108
109        let zt = DoubleDouble::from_exact_add(f64::from_bits(ix), f.hi);
110        f.hi = zt.hi;
111        f.lo += zt.lo;
112
113        zfh = to_denormal(f.to_f64());
114    }
115    zfh
116}
117
118#[inline(always)]
119fn exp10_gen<B: ExpfBackend>(x: f64, backend: B) -> f64 {
120    let mut ix = x.to_bits();
121    let aix = ix & 0x7fff_ffff_ffff_ffff;
122    if aix > 0x40734413509f79feu64 {
123        // |x| > 0x40734413509f79fe
124        if aix > 0x7ff0000000000000u64 {
125            return x + x;
126        } // nan
127        if aix == 0x7ff0000000000000u64 {
128            return if (ix >> 63) != 0 { 0.0 } else { x };
129        }
130        if (ix >> 63) == 0 {
131            return f64::from_bits(0x7fe0000000000000) * 2.0; // x > 308.255
132        }
133        if aix > 0x407439b746e36b52u64 {
134            // x < -323.607
135            return black_box(f64::from_bits(0x0018000000000000))
136                * black_box(f64::from_bits(0x3c80000000000000));
137        }
138    }
139
140    // check x integer to avoid a spurious inexact exception
141    if ix.wrapping_shl(16) == 0 && (aix >> 48) <= 0x4036 {
142        let kx = backend.round_ties_even(x);
143        if kx == x {
144            let k = kx as i64;
145            if k >= 0 {
146                let mut r = 1.0;
147                for _ in 0..k {
148                    r *= 10.0;
149                }
150                return r;
151            }
152        }
153    }
154    /* avoid spurious underflow: for |x| <= 2.41082e-17
155    exp10(x) rounds to 1 to nearest */
156    if aix <= 0x3c7bcb7b1526e50eu64 {
157        return 1.0 + x; // |x| <= 2.41082e-17
158    }
159    let t = backend.round_ties_even(f64::from_bits(0x40ca934f0979a371) * x);
160    let jt: i64 = unsafe { t.to_int_unchecked::<i64>() }; // t is already integer this is just a conversion
161    let i1 = jt & 0x3f;
162    let i0 = (jt >> 6) & 0x3f;
163    let ie = jt >> 12;
164    let t00 = EXP_REDUCE_T0[i0 as usize];
165    let t01 = EXP_REDUCE_T1[i1 as usize];
166    let t0 = DoubleDouble::from_bit_pair(t00);
167    let t1 = DoubleDouble::from_bit_pair(t01);
168    let mut tz = backend.quick_mult(t0, t1);
169    const L0: f64 = f64::from_bits(0x3f13441350800000);
170    const L1: f64 = f64::from_bits(0x3d1f79fef311f12b);
171    let dx = backend.fma(-L1, t, backend.fma(-L0, t, x));
172    let dx2 = dx * dx;
173
174    const CH: [u64; 4] = [
175        0x40026bb1bbb55516,
176        0x40053524c73cea69,
177        0x4000470591fd74e1,
178        0x3ff2bd760a1f32a5,
179    ];
180
181    let p0 = backend.fma(dx, f64::from_bits(CH[1]), f64::from_bits(CH[0]));
182    let p1 = backend.fma(dx, f64::from_bits(CH[3]), f64::from_bits(CH[2]));
183
184    let p = backend.fma(dx2, p1, p0);
185
186    let mut fh = tz.hi;
187    let fx = tz.hi * dx;
188    let mut fl = backend.fma(fx, p, tz.lo);
189    const EPS: f64 = 1.63e-19;
190    if ix < 0xc0733a7146f72a42u64 {
191        // x > -307.653
192        // x > -0x1.33a7146f72a42p+8
193        let ub = fh + (fl + EPS);
194        let lb = fh + (fl - EPS);
195
196        if lb != ub {
197            return as_exp10_accurate(x);
198        }
199        fh = fast_ldexp(fh + fl, ie as i32);
200    } else {
201        // x <= -307.653: exp10(x) < 2^-1022
202        ix = 1u64.wrapping_sub(ie as u64).wrapping_shl(52);
203        tz = DoubleDouble::from_exact_add(f64::from_bits(ix), fh);
204        fl += tz.lo;
205
206        let ub = fh + (fl + EPS);
207        let lb = fh + (fl - EPS);
208
209        if lb != ub {
210            return as_exp10_accurate(x);
211        }
212
213        fh = to_denormal(fh + fl);
214    }
215    fh
216}
217
218#[cfg(any(target_arch = "x86", target_arch = "x86_64"))]
219#[target_feature(enable = "avx", enable = "fma")]
220unsafe fn exp10_fma_impl(x: f64) -> f64 {
221    use crate::exponents::expf::FmaBackend;
222    exp10_gen(x, FmaBackend {})
223}
224
225/// Computes exp10
226///
227/// Max found ULP 0.5
228pub fn f_exp10(x: f64) -> f64 {
229    #[cfg(not(any(target_arch = "x86", target_arch = "x86_64")))]
230    {
231        exp10_gen(x, GenericExpfBackend {})
232    }
233    #[cfg(any(target_arch = "x86", target_arch = "x86_64"))]
234    {
235        use std::sync::OnceLock;
236        static EXECUTOR: OnceLock<unsafe fn(f64) -> f64> = OnceLock::new();
237        let q = EXECUTOR.get_or_init(|| {
238            if std::arch::is_x86_feature_detected!("avx")
239                && std::arch::is_x86_feature_detected!("fma")
240            {
241                exp10_fma_impl
242            } else {
243                fn def_exp10(x: f64) -> f64 {
244                    exp10_gen(x, GenericExpfBackend {})
245                }
246                def_exp10
247            }
248        });
249        unsafe { q(x) }
250    }
251}
252
253#[cfg(test)]
254mod tests {
255    use super::*;
256
257    #[test]
258    fn test_exp10f() {
259        assert_eq!(f_exp10(-3.370739843267434), 0.00042585343701025656);
260        assert_eq!(f_exp10(1.), 10.0);
261        assert_eq!(f_exp10(2.), 100.0);
262        assert_eq!(f_exp10(3.), 1000.0);
263        assert_eq!(f_exp10(4.), 10000.0);
264        assert_eq!(f_exp10(5.), 100000.0);
265        assert_eq!(f_exp10(6.), 1000000.0);
266        assert_eq!(f_exp10(7.), 10000000.0);
267        assert_eq!(f_exp10(f64::INFINITY), f64::INFINITY);
268        assert_eq!(f_exp10(f64::NEG_INFINITY), 0.);
269        assert!(f_exp10(f64::NAN).is_nan());
270    }
271}
pxfm/exponents/exp10.rs

pxfm/exponents/
exp10.rs
