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

pxfm/logs/
log10f.rs
