pxfm/exponents/
exp10m1f.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::exponents::exp10f::EXP10F_COEFFS;
31use crate::polyeval::f_polyeval3;
32
33#[cold]
34fn exp10m1f_small(x: f32) -> f32 {
35    let dx = x as f64;
36    let dx_sq = dx * dx;
37    let c0 = dx * f64::from_bits(EXP10F_COEFFS[0]);
38    let c1 = f_fmla(
39        dx,
40        f64::from_bits(EXP10F_COEFFS[2]),
41        f64::from_bits(EXP10F_COEFFS[1]),
42    );
43    let c2 = f_fmla(
44        dx,
45        f64::from_bits(EXP10F_COEFFS[4]),
46        f64::from_bits(EXP10F_COEFFS[3]),
47    );
48    // 10^dx - 1 ~ (1 + COEFFS[0] * dx + ... + COEFFS[4] * dx^5) - 1
49    //           = COEFFS[0] * dx + ... + COEFFS[4] * dx^5
50    f_polyeval3(dx_sq, c0, c1, c2) as f32
51}
52
53/// Computes 10^x - 1
54///
55/// Max ULP 0.5
56#[inline]
57pub fn f_exp10m1f(x: f32) -> f32 {
58    let x_u = x.to_bits();
59    let x_abs = x_u & 0x7fff_ffffu32;
60
61    // When x >= log10(2^128), or x is nan
62    if x.is_sign_positive() && x_u >= 0x421a_209bu32 {
63        // x >= log10(2^128) and 10^x - 1 rounds to +inf, or x is +inf or nan
64        return x + f32::INFINITY;
65    }
66
67    if x_abs <= 0x3b9a_209bu32 {
68        // |x| <= 0.004703594
69        return exp10m1f_small(x);
70    }
71
72    // When x <= log10(2^-25), or x is nan
73    if x_u >= 0xc0f0d2f1 {
74        // exp10m1(-inf) = -1
75        if x.is_infinite() {
76            return -1.0;
77        }
78        // exp10m1(nan) = nan
79        if x.is_nan() {
80            return x;
81        }
82
83        if x_u == 0xc0f0d2f1 {
84            return f32::from_bits(0xbf7fffff); // -1.0f + 0x1.0p-24f
85        }
86        return -1.0;
87    }
88
89    // Exact outputs when x = 1, 2, ..., 10.
90    // Quick check mask: 0x800f'ffffU = ~(bits of 1.0f | ... | bits of 10.0f)
91    if x_u & 0x800f_ffffu32 == 0 {
92        match x_u {
93            0x3f800000u32 => return 9.0,             // x = 1.0f
94            0x40000000u32 => return 99.0,            // x = 2.0f
95            0x40400000u32 => return 999.0,           // x = 3.0f
96            0x40800000u32 => return 9_999.0,         // x = 4.0f
97            0x40a00000u32 => return 99_999.0,        // x = 5.0f
98            0x40c00000u32 => return 999_999.0,       // x = 6.0f
99            0x40e00000u32 => return 9_999_999.0,     // x = 7.0f
100            0x41000000u32 => return 99_999_999.0,    // x = 8.0f
101            0x41100000u32 => return 999_999_999.0,   // x = 9.0f
102            0x41200000u32 => return 9_999_999_999.0, // x = 10.0f
103            _ => {}
104        }
105    }
106
107    // Range reduction: 10^x = 2^(mid + hi) * 10^lo
108    //   rr = (2^(mid + hi), lo)
109    let rr = crate::exponents::exp10f::exp_b_range_reduc(x);
110
111    // The low part is approximated by a degree-5 minimax polynomial.
112    // 10^lo ~ 1 + COEFFS[0] * lo + ... + COEFFS[4] * lo^5
113    let lo_sq = rr.lo * rr.lo;
114    let c0 = f_fmla(rr.lo, f64::from_bits(EXP10F_COEFFS[0]), 1.0);
115    let c1 = f_fmla(
116        rr.lo,
117        f64::from_bits(EXP10F_COEFFS[2]),
118        f64::from_bits(EXP10F_COEFFS[1]),
119    );
120    let c2 = f_fmla(
121        rr.lo,
122        f64::from_bits(EXP10F_COEFFS[4]),
123        f64::from_bits(EXP10F_COEFFS[3]),
124    );
125    let exp10_lo = f_polyeval3(lo_sq, c0, c1, c2);
126    // 10^x - 1 = 2^(mid + hi) * 10^lo - 1
127    //          ~ mh * exp10_lo - 1
128    f_fmla(exp10_lo, rr.hi, -1.0) as f32
129}
130
131#[cfg(test)]
132mod tests {
133    use super::*;
134
135    #[test]
136    fn test_exp10m1f() {
137        assert_eq!(f_exp10m1f(0.0), 0.0);
138        assert_eq!(f_exp10m1f(1.0), 9.0);
139        assert_eq!(f_exp10m1f(1.5), 30.622776);
140    }
141}
pxfm/exponents/exp10m1f.rs

pxfm/exponents/
exp10m1f.rs
