/rust/registry/src/index.crates.io-1949cf8c6b5b557f/pxfm-0.1.30/src/exponents/exp10m1f.rs
Line | Count | Source |
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 | | */ |
29 | | use crate::exponents::exp10f::EXP10F_COEFFS; |
30 | | use crate::exponents::expf::{ExpfBackend, GenericExpfBackend}; |
31 | | |
32 | | #[cold] |
33 | | #[inline(always)] |
34 | 0 | fn exp10m1f_small<B: ExpfBackend>(x: f32, backend: &B) -> f32 { |
35 | 0 | let dx = x as f64; |
36 | 0 | let dx_sq = dx * dx; |
37 | 0 | let c0 = dx * f64::from_bits(EXP10F_COEFFS[0]); |
38 | 0 | let c1 = backend.fma( |
39 | 0 | dx, |
40 | 0 | f64::from_bits(EXP10F_COEFFS[2]), |
41 | 0 | f64::from_bits(EXP10F_COEFFS[1]), |
42 | | ); |
43 | 0 | let c2 = backend.fma( |
44 | 0 | dx, |
45 | 0 | f64::from_bits(EXP10F_COEFFS[4]), |
46 | 0 | 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 | 0 | backend.polyeval3(dx_sq, c0, c1, c2) as f32 |
51 | 0 | } Unexecuted instantiation: pxfm::exponents::exp10m1f::exp10m1f_small::<pxfm::exponents::expf::FmaBackend> Unexecuted instantiation: pxfm::exponents::exp10m1f::exp10m1f_small::<pxfm::exponents::expf::GenericExpfBackend> |
52 | | |
53 | | #[inline(always)] |
54 | 0 | fn exp10m1f_gen<B: ExpfBackend>(x: f32, backend: B) -> f32 { |
55 | 0 | let x_u = x.to_bits(); |
56 | 0 | let x_abs = x_u & 0x7fff_ffffu32; |
57 | | |
58 | | // When x >= log10(2^128), or x is nan |
59 | 0 | 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 | 0 | return x + f32::INFINITY; |
62 | 0 | } |
63 | | |
64 | 0 | if x_abs <= 0x3b9a_209bu32 { |
65 | | // |x| <= 0.004703594 |
66 | 0 | return exp10m1f_small(x, &backend); |
67 | 0 | } |
68 | | |
69 | | // When x <= log10(2^-25), or x is nan |
70 | 0 | if x_u >= 0xc0f0d2f1 { |
71 | | // exp10m1(-inf) = -1 |
72 | 0 | if x.is_infinite() { |
73 | 0 | return -1.0; |
74 | 0 | } |
75 | | // exp10m1(nan) = nan |
76 | 0 | if x.is_nan() { |
77 | 0 | return x; |
78 | 0 | } |
79 | | |
80 | 0 | if x_u == 0xc0f0d2f1 { |
81 | 0 | return f32::from_bits(0xbf7fffff); // -1.0f + 0x1.0p-24f |
82 | 0 | } |
83 | 0 | return -1.0; |
84 | 0 | } |
85 | | |
86 | | // Exact outputs when x = 1, 2, ..., 10. |
87 | | // Quick check mask: 0x800f'ffffU = ~(bits of 1.0f | ... | bits of 10.0f) |
88 | 0 | if x_u & 0x800f_ffffu32 == 0 { |
89 | 0 | match x_u { |
90 | 0 | 0x3f800000u32 => return 9.0, // x = 1.0f |
91 | 0 | 0x40000000u32 => return 99.0, // x = 2.0f |
92 | 0 | 0x40400000u32 => return 999.0, // x = 3.0f |
93 | 0 | 0x40800000u32 => return 9_999.0, // x = 4.0f |
94 | 0 | 0x40a00000u32 => return 99_999.0, // x = 5.0f |
95 | 0 | 0x40c00000u32 => return 999_999.0, // x = 6.0f |
96 | 0 | 0x40e00000u32 => return 9_999_999.0, // x = 7.0f |
97 | 0 | 0x41000000u32 => return 99_999_999.0, // x = 8.0f |
98 | 0 | 0x41100000u32 => return 999_999_999.0, // x = 9.0f |
99 | 0 | 0x41200000u32 => return 9_999_999_999.0, // x = 10.0f |
100 | 0 | _ => {} |
101 | | } |
102 | 0 | } |
103 | | |
104 | | // Range reduction: 10^x = 2^(mid + hi) * 10^lo |
105 | | // rr = (2^(mid + hi), lo) |
106 | 0 | 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 | 0 | let lo_sq = rr.lo * rr.lo; |
111 | 0 | let c0 = backend.fma(rr.lo, f64::from_bits(EXP10F_COEFFS[0]), 1.0); |
112 | 0 | let c1 = backend.fma( |
113 | 0 | rr.lo, |
114 | 0 | f64::from_bits(EXP10F_COEFFS[2]), |
115 | 0 | f64::from_bits(EXP10F_COEFFS[1]), |
116 | | ); |
117 | 0 | let c2 = backend.fma( |
118 | 0 | rr.lo, |
119 | 0 | f64::from_bits(EXP10F_COEFFS[4]), |
120 | 0 | f64::from_bits(EXP10F_COEFFS[3]), |
121 | | ); |
122 | 0 | 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 | 0 | backend.fma(exp10_lo, rr.hi, -1.0) as f32 |
126 | 0 | } Unexecuted instantiation: pxfm::exponents::exp10m1f::exp10m1f_gen::<pxfm::exponents::expf::FmaBackend> Unexecuted instantiation: pxfm::exponents::exp10m1f::exp10m1f_gen::<pxfm::exponents::expf::GenericExpfBackend> |
127 | | |
128 | | #[cfg(any(target_arch = "x86", target_arch = "x86_64"))] |
129 | | #[target_feature(enable = "avx", enable = "fma")] |
130 | 0 | unsafe fn exp10f_fma_impl(x: f32) -> f32 { |
131 | | use crate::exponents::expf::FmaBackend; |
132 | 0 | exp10m1f_gen(x, FmaBackend {}) |
133 | 0 | } |
134 | | |
135 | | /// Computes 10^x - 1 |
136 | | /// |
137 | | /// Max ULP 0.5 |
138 | | #[inline] |
139 | 0 | pub 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 | 0 | let q = EXECUTOR.get_or_init(|| { |
149 | 0 | if std::arch::is_x86_feature_detected!("avx") |
150 | 0 | && std::arch::is_x86_feature_detected!("fma") |
151 | | { |
152 | 0 | exp10f_fma_impl |
153 | | } else { |
154 | 0 | fn def_expf(x: f32) -> f32 { |
155 | 0 | exp10m1f_gen(x, GenericExpfBackend {}) |
156 | 0 | } |
157 | 0 | def_expf |
158 | | } |
159 | 0 | }); |
160 | 0 | unsafe { q(x) } |
161 | | } |
162 | 0 | } |
163 | | |
164 | | #[cfg(test)] |
165 | | mod 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 | | } |