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::exponents::exp10f::EXP10F_COEFFS;
30use crate::exponents::expf::{ExpfBackend, GenericExpfBackend};
31
32#[cold]
33#[inline(always)]
34fn exp10m1f_small<B: ExpfBackend>(x: f32, backend: &B) -> 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 = backend.fma(
39 dx,
40 f64::from_bits(EXP10F_COEFFS[2]),
41 f64::from_bits(EXP10F_COEFFS[1]),
42 );
43 let c2 = backend.fma(
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 backend.polyeval3(dx_sq, c0, c1, c2) as f32
51}
52
53#[inline(always)]
54fn exp10m1f_gen<B: ExpfBackend>(x: f32, backend: B) -> f32 {
55 let x_u = x.to_bits();
56 let x_abs = x_u & 0x7fff_ffffu32;
57
58 // When x >= log10(2^128), or x is nan
59 if x.is_sign_positive() && x_u >= 0x421a_209bu32 {
60 // x >= log10(2^128) and 10^x - 1 rounds to +inf, or x is +inf or nan
61 return x + f32::INFINITY;
62 }
63
64 if x_abs <= 0x3b9a_209bu32 {
65 // |x| <= 0.004703594
66 return exp10m1f_small(x, &backend);
67 }
68
69 // When x <= log10(2^-25), or x is nan
70 if x_u >= 0xc0f0d2f1 {
71 // exp10m1(-inf) = -1
72 if x.is_infinite() {
73 return -1.0;
74 }
75 // exp10m1(nan) = nan
76 if x.is_nan() {
77 return x;
78 }
79
80 if x_u == 0xc0f0d2f1 {
81 return f32::from_bits(0xbf7fffff); // -1.0f + 0x1.0p-24f
82 }
83 return -1.0;
84 }
85
86 // Exact outputs when x = 1, 2, ..., 10.
87 // Quick check mask: 0x800f'ffffU = ~(bits of 1.0f | ... | bits of 10.0f)
88 if x_u & 0x800f_ffffu32 == 0 {
89 match x_u {
90 0x3f800000u32 => return 9.0, // x = 1.0f
91 0x40000000u32 => return 99.0, // x = 2.0f
92 0x40400000u32 => return 999.0, // x = 3.0f
93 0x40800000u32 => return 9_999.0, // x = 4.0f
94 0x40a00000u32 => return 99_999.0, // x = 5.0f
95 0x40c00000u32 => return 999_999.0, // x = 6.0f
96 0x40e00000u32 => return 9_999_999.0, // x = 7.0f
97 0x41000000u32 => return 99_999_999.0, // x = 8.0f
98 0x41100000u32 => return 999_999_999.0, // x = 9.0f
99 0x41200000u32 => return 9_999_999_999.0, // x = 10.0f
100 _ => {}
101 }
102 }
103
104 // Range reduction: 10^x = 2^(mid + hi) * 10^lo
105 // rr = (2^(mid + hi), lo)
106 let rr = crate::exponents::exp10f::exp_b_range_reduc(x, &backend);
107
108 // The low part is approximated by a degree-5 minimax polynomial.
109 // 10^lo ~ 1 + COEFFS[0] * lo + ... + COEFFS[4] * lo^5
110 let lo_sq = rr.lo * rr.lo;
111 let c0 = backend.fma(rr.lo, f64::from_bits(EXP10F_COEFFS[0]), 1.0);
112 let c1 = backend.fma(
113 rr.lo,
114 f64::from_bits(EXP10F_COEFFS[2]),
115 f64::from_bits(EXP10F_COEFFS[1]),
116 );
117 let c2 = backend.fma(
118 rr.lo,
119 f64::from_bits(EXP10F_COEFFS[4]),
120 f64::from_bits(EXP10F_COEFFS[3]),
121 );
122 let exp10_lo = backend.polyeval3(lo_sq, c0, c1, c2);
123 // 10^x - 1 = 2^(mid + hi) * 10^lo - 1
124 // ~ mh * exp10_lo - 1
125 backend.fma(exp10_lo, rr.hi, -1.0) as f32
126}
127
128#[cfg(any(target_arch = "x86", target_arch = "x86_64"))]
129#[target_feature(enable = "avx", enable = "fma")]
130unsafe fn exp10f_fma_impl(x: f32) -> f32 {
131 use crate::exponents::expf::FmaBackend;
132 exp10m1f_gen(x, FmaBackend {})
133}
134
135/// Computes 10^x - 1
136///
137/// Max ULP 0.5
138#[inline]
139pub fn f_exp10m1f(x: f32) -> f32 {
140 #[cfg(not(any(target_arch = "x86", target_arch = "x86_64")))]
141 {
142 exp10m1f_gen(x, GenericExpfBackend {})
143 }
144 #[cfg(any(target_arch = "x86", target_arch = "x86_64"))]
145 {
146 use std::sync::OnceLock;
147 static EXECUTOR: OnceLock<unsafe fn(f32) -> f32> = OnceLock::new();
148 let q = EXECUTOR.get_or_init(|| {
149 if std::arch::is_x86_feature_detected!("avx")
150 && std::arch::is_x86_feature_detected!("fma")
151 {
152 exp10f_fma_impl
153 } else {
154 fn def_expf(x: f32) -> f32 {
155 exp10m1f_gen(x, GenericExpfBackend {})
156 }
157 def_expf
158 }
159 });
160 unsafe { q(x) }
161 }
162}
163
164#[cfg(test)]
165mod tests {
166 use super::*;
167
168 #[test]
169 fn test_exp10m1f() {
170 assert_eq!(f_exp10m1f(0.0), 0.0);
171 assert_eq!(f_exp10m1f(1.0), 9.0);
172 assert_eq!(f_exp10m1f(1.5), 30.622776);
173 }
174}