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::polyeval::f_polyeval3;
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/// Logarithm of base 10
164///
165/// Max found ULP 0.5
166pub fn f_log10f(x: f32) -> f32 {
167    let mut x_u = x.to_bits();
168
169    const E_BIAS: u32 = (1u32 << (8 - 1u32)) - 1u32;
170
171    let mut m = -(E_BIAS as i32);
172    if x_u < f32::MIN_POSITIVE.to_bits() || x_u > f32::MAX.to_bits() {
173        if x == 0.0 {
174            return f32::NEG_INFINITY;
175        }
176        if x_u == 0x80000000u32 {
177            return f32::NEG_INFINITY;
178        }
179        if x.is_sign_negative() && !x.is_nan() {
180            return f32::NAN + x;
181        }
182        // x is +inf or nan
183        if x.is_nan() || x.is_infinite() {
184            return x + x;
185        }
186        // Normalize denormal inputs.
187        x_u = (x * f64::from_bits(0x4160000000000000) as f32).to_bits();
188        m -= 23;
189    }
190
191    m = m.wrapping_add(x_u.wrapping_shr(23) as i32);
192    let index = (x_u >> 16) & 0x7F;
193    x_u = set_exponent_f32(x_u, 0x7F);
194
195    let v;
196    let u = f32::from_bits(x_u);
197
198    #[cfg(any(
199        all(
200            any(target_arch = "x86", target_arch = "x86_64"),
201            target_feature = "fma"
202        ),
203        target_arch = "aarch64"
204    ))]
205    {
206        v = f_fmlaf(
207            u,
208            f32::from_bits(crate::logs::logf::LOG_REDUCTION_F32.0[index as usize]),
209            -1.0,
210        ) as f64; // Exact.
211    }
212    #[cfg(not(any(
213        all(
214            any(target_arch = "x86", target_arch = "x86_64"),
215            target_feature = "fma"
216        ),
217        target_arch = "aarch64"
218    )))]
219    {
220        v = f_fmla(
221            u as f64,
222            f32::from_bits(crate::logs::logf::LOG_REDUCTION_F32.0[index as usize]) as f64,
223            -1.0,
224        ); // Exact
225    }
226    // Degree-5 polynomial approximation of log10 generated by:
227    // > P = fpminimax(log10(1 + x)/x, 4, [|D...|], [-2^-8, 2^-7]);
228    const COEFFS: [u64; 5] = [
229        0x3fdbcb7b1526e2e5,
230        0xbfcbcb7b1528d43d,
231        0x3fc287a77eb4ca0d,
232        0xbfbbcb8110a181b5,
233        0x3fb60e7e3e747129,
234    ];
235    let v2 = v * v; // Exact
236    let p2 = f_fmla(v, f64::from_bits(COEFFS[4]), f64::from_bits(COEFFS[3]));
237    let p1 = f_fmla(v, f64::from_bits(COEFFS[2]), f64::from_bits(COEFFS[1]));
238    let p0 = f_fmla(
239        v,
240        f64::from_bits(COEFFS[0]),
241        f64::from_bits(LOG10_R[index as usize]),
242    );
243    const LOG_10_2: f64 = f64::from_bits(0x3fd34413509f79ff);
244    let r = f_fmla(m as f64, LOG_10_2, f_polyeval3(v2, p0, p1, p2));
245    r as f32
246}
247
248#[cfg(test)]
249mod tests {
250    use super::*;
251
252    #[test]
253    fn test_log10f() {
254        assert_eq!(f_log10f(0.35), -0.45593196f32);
255        assert_eq!(f_log10f(0.9), -4.5757502e-2);
256        assert_eq!(f_log10f(10.), 1.);
257        assert_eq!(f_log10f(100.), 2.);
258        assert_eq!(f_log10f(0.), f32::NEG_INFINITY);
259        assert!(f_log10f(-1.).is_nan());
260        assert!(f_log10f(f32::NAN).is_nan());
261        assert!(f_log10f(f32::NEG_INFINITY).is_nan());
262        assert_eq!(f_log10f(f32::INFINITY), f32::INFINITY);
263    }
264}
pxfm/logs/log10f.rs

pxfm/logs/
log10f.rs
