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::common::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;
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).round_ties_even_finite();
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/// Computes exp10
119///
120/// Max found ULP 0.5
121pub fn f_exp10(x: f64) -> f64 {
122    let mut ix = x.to_bits();
123    let aix = ix & 0x7fff_ffff_ffff_ffff;
124    if aix > 0x40734413509f79feu64 {
125        // |x| > 0x40734413509f79fe
126        if aix > 0x7ff0000000000000u64 {
127            return x + x;
128        } // nan
129        if aix == 0x7ff0000000000000u64 {
130            return if (ix >> 63) != 0 { 0.0 } else { x };
131        }
132        if (ix >> 63) == 0 {
133            return f64::from_bits(0x7fe0000000000000) * 2.0; // x > 308.255
134        }
135        if aix > 0x407439b746e36b52u64 {
136            // x < -323.607
137            return black_box(f64::from_bits(0x0018000000000000))
138                * black_box(f64::from_bits(0x3c80000000000000));
139        }
140    }
141
142    // check x integer to avoid a spurious inexact exception
143    if ix.wrapping_shl(16) == 0 && (aix >> 48) <= 0x4036 {
144        let kx = x.round_ties_even_finite();
145        if kx == x {
146            let k = kx as i64;
147            if k >= 0 {
148                let mut r = 1.0;
149                for _ in 0..k {
150                    r *= 10.0;
151                }
152                return r;
153            }
154        }
155    }
156    /* avoid spurious underflow: for |x| <= 2.41082e-17
157    exp10(x) rounds to 1 to nearest */
158    if aix <= 0x3c7bcb7b1526e50eu64 {
159        return 1.0 + x; // |x| <= 2.41082e-17
160    }
161    let t = (f64::from_bits(0x40ca934f0979a371) * x).round_ties_even_finite();
162    let jt: i64 = unsafe { t.to_int_unchecked::<i64>() }; // t is already integer this is just a conversion
163    let i1 = jt & 0x3f;
164    let i0 = (jt >> 6) & 0x3f;
165    let ie = jt >> 12;
166    let t00 = EXP_REDUCE_T0[i0 as usize];
167    let t01 = EXP_REDUCE_T1[i1 as usize];
168    let t0 = DoubleDouble::from_bit_pair(t00);
169    let t1 = DoubleDouble::from_bit_pair(t01);
170    let mut tz = DoubleDouble::quick_mult(t0, t1);
171    const L0: f64 = f64::from_bits(0x3f13441350800000);
172    const L1: f64 = f64::from_bits(0x3d1f79fef311f12b);
173    let dx = f_fmla(-L1, t, f_fmla(-L0, t, x));
174    let dx2 = dx * dx;
175
176    const CH: [u64; 4] = [
177        0x40026bb1bbb55516,
178        0x40053524c73cea69,
179        0x4000470591fd74e1,
180        0x3ff2bd760a1f32a5,
181    ];
182
183    let p0 = f_fmla(dx, f64::from_bits(CH[1]), f64::from_bits(CH[0]));
184    let p1 = f_fmla(dx, f64::from_bits(CH[3]), f64::from_bits(CH[2]));
185
186    let p = f_fmla(dx2, p1, p0);
187
188    let mut fh = tz.hi;
189    let fx = tz.hi * dx;
190    let mut fl = f_fmla(fx, p, tz.lo);
191    const EPS: f64 = 1.63e-19;
192    if ix < 0xc0733a7146f72a42u64 {
193        // x > -307.653
194        // x > -0x1.33a7146f72a42p+8
195        let ub = fh + (fl + EPS);
196        let lb = fh + (fl - EPS);
197
198        if lb != ub {
199            return as_exp10_accurate(x);
200        }
201        fh = fast_ldexp(fh + fl, ie as i32);
202    } else {
203        // x <= -307.653: exp10(x) < 2^-1022
204        ix = 1u64.wrapping_sub(ie as u64).wrapping_shl(52);
205        tz = DoubleDouble::from_exact_add(f64::from_bits(ix), fh);
206        fl += tz.lo;
207
208        let ub = fh + (fl + EPS);
209        let lb = fh + (fl - EPS);
210
211        if lb != ub {
212            return as_exp10_accurate(x);
213        }
214
215        fh = to_denormal(fh + fl);
216    }
217    fh
218}
219
220#[cfg(test)]
221mod tests {
222    use super::*;
223
224    #[test]
225    fn test_exp10f() {
226        assert_eq!(f_exp10(-3.370739843267434), 0.00042585343701025656);
227        assert_eq!(f_exp10(1.), 10.0);
228        assert_eq!(f_exp10(2.), 100.0);
229        assert_eq!(f_exp10(3.), 1000.0);
230        assert_eq!(f_exp10(4.), 10000.0);
231        assert_eq!(f_exp10(5.), 100000.0);
232        assert_eq!(f_exp10(6.), 1000000.0);
233        assert_eq!(f_exp10(7.), 10000000.0);
234        assert_eq!(f_exp10(f64::INFINITY), f64::INFINITY);
235        assert_eq!(f_exp10(f64::NEG_INFINITY), 0.);
236        assert!(f_exp10(f64::NAN).is_nan());
237    }
238}
pxfm/exponents/exp10.rs

pxfm/exponents/
exp10.rs
