pxfm/logs/
log10f.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::*;
30use crate::logs::logf::{GenLogfBackend, LogfBackend};
31
32static LOG10_R: [u64; 128] = [
33    0x0000000000000000,
34    0x3f6be76bd77b4fc3,
35    0x3f7c03a80ae5e054,
36    0x3f851824c7587eb0,
37    0x3f8c3d0837784c41,
38    0x3f91b85d6044e9ae,
39    0x3f9559bd2406c3ba,
40    0x3f9902c31d62a843,
41    0x3f9cb38fccd8bfdb,
42    0x3f9e8eeb09f2f6cb,
43    0x3fa125d0432ea20e,
44    0x3fa30838cdc2fbfd,
45    0x3fa3faf7c663060e,
46    0x3fa5e3966b7e9295,
47    0x3fa7d070145f4fd7,
48    0x3fa8c878eeb05074,
49    0x3faabbcebd84fca0,
50    0x3fabb7209d1e24e5,
51    0x3fadb11ed766abf4,
52    0x3faeafd05035bd3b,
53    0x3fb0585283764178,
54    0x3fb0d966cc6500fa,
55    0x3fb1dd5460c8b16f,
56    0x3fb2603072a25f82,
57    0x3fb367ba3aaa1883,
58    0x3fb3ec6ad5407868,
59    0x3fb4f7aad9bbcbaf,
60    0x3fb57e3d47c3af7b,
61    0x3fb605735ee985f1,
62    0x3fb715d0ce367afc,
63    0x3fb79efb57b0f803,
64    0x3fb828cfed29a215,
65    0x3fb93e7de0fc3e80,
66    0x3fb9ca5aa1729f45,
67    0x3fba56e8325f5c87,
68    0x3fbae4285509950b,
69    0x3fbb721cd17157e3,
70    0x3fbc902a19e65111,
71    0x3fbd204698cb42bd,
72    0x3fbdb11ed766abf4,
73    0x3fbe42b4c16caaf3,
74    0x3fbed50a4a26eafc,
75    0x3fbffbfc2bbc7803,
76    0x3fc0484e4942aa43,
77    0x3fc093025a19976c,
78    0x3fc0de1b56356b04,
79    0x3fc1299a4fb3e306,
80    0x3fc175805d1587c1,
81    0x3fc1c1ce9955c0c6,
82    0x3fc20e8624038fed,
83    0x3fc25ba8215af7fc,
84    0x3fc2a935ba5f1479,
85    0x3fc2f7301cf4e87b,
86    0x3fc345987bfeea91,
87    0x3fc394700f7953fd,
88    0x3fc3e3b8149739d4,
89    0x3fc43371cde076c2,
90    0x3fc4839e83506c87,
91    0x3fc4d43f8275a483,
92    0x3fc525561e9256ee,
93    0x3fc576e3b0bde0a7,
94    0x3fc5c8e998072fe2,
95    0x3fc61b6939983048,
96    0x3fc66e6400da3f77,
97    0x3fc6c1db5f9bb336,
98    0x3fc6c1db5f9bb336,
99    0x3fc715d0ce367afc,
100    0x3fc76a45cbb7e6ff,
101    0x3fc7bf3bde099f30,
102    0x3fc814b4921bd52b,
103    0x3fc86ab17c10bc7f,
104    0x3fc86ab17c10bc7f,
105    0x3fc8c13437695532,
106    0x3fc9183e673394fa,
107    0x3fc96fd1b639fc09,
108    0x3fc9c7efd734a2f9,
109    0x3fca209a84fbcff8,
110    0x3fca209a84fbcff8,
111    0x3fca79d382bc21d9,
112    0x3fcad39c9c2c6080,
113    0x3fcb2df7a5c50299,
114    0x3fcb2df7a5c50299,
115    0x3fcb88e67cf97980,
116    0x3fcbe46b087354bc,
117    0x3fcc4087384f4f80,
118    0x3fcc4087384f4f80,
119    0x3fcc9d3d065c5b42,
120    0x3fccfa8e765cbb72,
121    0x3fccfa8e765cbb72,
122    0x3fcd587d96494759,
123    0x3fcdb70c7e96e7f3,
124    0x3fcdb70c7e96e7f3,
125    0x3fce163d527e68cf,
126    0x3fce76124046b3f3,
127    0x3fce76124046b3f3,
128    0x3fced68d819191fc,
129    0x3fcf37b15bab08d1,
130    0x3fcf37b15bab08d1,
131    0x3fcf99801fdb749d,
132    0x3fcffbfc2bbc7803,
133    0x3fcffbfc2bbc7803,
134    0x3fd02f93f4c87101,
135    0x3fd06182e84fd4ac,
136    0x3fd06182e84fd4ac,
137    0x3fd093cc32c90f84,
138    0x3fd093cc32c90f84,
139    0x3fd0c6711d6abd7a,
140    0x3fd0f972f87ff3d6,
141    0x3fd0f972f87ff3d6,
142    0x3fd12cd31b9c99ff,
143    0x3fd12cd31b9c99ff,
144    0x3fd16092e5d3a9a6,
145    0x3fd194b3bdef6b9e,
146    0x3fd194b3bdef6b9e,
147    0x3fd1c93712abc7ff,
148    0x3fd1c93712abc7ff,
149    0x3fd1fe1e5af2c141,
150    0x3fd1fe1e5af2c141,
151    0x3fd2336b161b3337,
152    0x3fd2336b161b3337,
153    0x3fd2691ecc29f042,
154    0x3fd2691ecc29f042,
155    0x3fd29f3b0e15584b,
156    0x3fd29f3b0e15584b,
157    0x3fd2d5c1760b86bb,
158    0x3fd2d5c1760b86bb,
159    0x3fd30cb3a7bb3625,
160    0x3fd34413509f79ff,
161];
162
163#[inline(always)]
164fn log10f_gen<B: LogfBackend>(x: f32, backend: B) -> f32 {
165    let mut x_u = x.to_bits();
166
167    const E_BIAS: u32 = (1u32 << (8 - 1u32)) - 1u32;
168
169    let mut m = -(E_BIAS as i32);
170    if x_u < f32::MIN_POSITIVE.to_bits() || x_u > f32::MAX.to_bits() {
171        if x == 0.0 {
172            return f32::NEG_INFINITY;
173        }
174        if x_u == 0x80000000u32 {
175            return f32::NEG_INFINITY;
176        }
177        if x.is_sign_negative() && !x.is_nan() {
178            return f32::NAN + x;
179        }
180        // x is +inf or nan
181        if x.is_nan() || x.is_infinite() {
182            return x + x;
183        }
184        // Normalize denormal inputs.
185        x_u = (x * f64::from_bits(0x4160000000000000) as f32).to_bits();
186        m -= 23;
187    }
188
189    m = m.wrapping_add(x_u.wrapping_shr(23) as i32);
190    let index = (x_u >> 16) & 0x7F;
191    x_u = set_exponent_f32(x_u, 0x7F);
192
193    let u = f32::from_bits(x_u);
194
195    let v = if B::HAS_FMA {
196        backend.fmaf(
197            u,
198            f32::from_bits(crate::logs::logf::LOG_REDUCTION_F32.0[index as usize]),
199            -1.0,
200        ) as f64 // Exact.
201    } else {
202        backend.fma(
203            u as f64,
204            f32::from_bits(crate::logs::logf::LOG_REDUCTION_F32.0[index as usize]) as f64,
205            -1.0,
206        ) // Exact
207    };
208    // Degree-5 polynomial approximation of log10 generated by:
209    // > P = fpminimax(log10(1 + x)/x, 4, [|D...|], [-2^-8, 2^-7]);
210    const COEFFS: [u64; 5] = [
211        0x3fdbcb7b1526e2e5,
212        0xbfcbcb7b1528d43d,
213        0x3fc287a77eb4ca0d,
214        0xbfbbcb8110a181b5,
215        0x3fb60e7e3e747129,
216    ];
217    let v2 = v * v; // Exact
218    let p2 = backend.fma(v, f64::from_bits(COEFFS[4]), f64::from_bits(COEFFS[3]));
219    let p1 = backend.fma(v, f64::from_bits(COEFFS[2]), f64::from_bits(COEFFS[1]));
220    let p0 = backend.fma(
221        v,
222        f64::from_bits(COEFFS[0]),
223        f64::from_bits(LOG10_R[index as usize]),
224    );
225    const LOG_10_2: f64 = f64::from_bits(0x3fd34413509f79ff);
226    let r = backend.fma(m as f64, LOG_10_2, backend.polyeval3(v2, p0, p1, p2));
227    r as f32
228}
229
230#[cfg(any(target_arch = "x86", target_arch = "x86_64"))]
231#[target_feature(enable = "avx", enable = "fma")]
232unsafe fn log10f_fma_impl(x: f32) -> f32 {
233    use crate::logs::logf::FmaLogfBackend;
234    log10f_gen(x, FmaLogfBackend {})
235}
236
237/// Logarithm of base 10
238///
239/// Max found ULP 0.5
240pub fn f_log10f(x: f32) -> f32 {
241    #[cfg(not(any(target_arch = "x86", target_arch = "x86_64")))]
242    {
243        log10f_gen(x, GenLogfBackend {})
244    }
245    #[cfg(any(target_arch = "x86", target_arch = "x86_64"))]
246    {
247        use std::sync::OnceLock;
248        static EXECUTOR: OnceLock<unsafe fn(f32) -> f32> = OnceLock::new();
249        let q = EXECUTOR.get_or_init(|| {
250            if std::arch::is_x86_feature_detected!("avx")
251                && std::arch::is_x86_feature_detected!("fma")
252            {
253                log10f_fma_impl
254            } else {
255                fn def_log10f(x: f32) -> f32 {
256                    log10f_gen(x, GenLogfBackend {})
257                }
258                def_log10f
259            }
260        });
261        unsafe { q(x) }
262    }
263}
264
265#[cfg(test)]
266mod tests {
267    use super::*;
268
269    #[test]
270    fn test_log10f() {
271        assert_eq!(f_log10f(0.35), -0.45593196f32);
272        assert_eq!(f_log10f(0.9), -4.5757502e-2);
273        assert_eq!(f_log10f(10.), 1.);
274        assert_eq!(f_log10f(100.), 2.);
275        assert_eq!(f_log10f(0.), f32::NEG_INFINITY);
276        assert!(f_log10f(-1.).is_nan());
277        assert!(f_log10f(f32::NAN).is_nan());
278        assert!(f_log10f(f32::NEG_INFINITY).is_nan());
279        assert_eq!(f_log10f(f32::INFINITY), f32::INFINITY);
280    }
281}
pxfm/logs/log10f.rs

pxfm/logs/
log10f.rs
