/rust/registry/src/index.crates.io-1949cf8c6b5b557f/pxfm-0.1.25/src/logs/log10f.rs

Source
/*
 * // Copyright (c) Radzivon Bartoshyk 7/2025. All rights reserved.
 * //
 * // Redistribution and use in source and binary forms, with or without modification,
 * // are permitted provided that the following conditions are met:
 * //
 * // 1.  Redistributions of source code must retain the above copyright notice, this
 * // list of conditions and the following disclaimer.
 * //
 * // 2.  Redistributions in binary form must reproduce the above copyright notice,
 * // this list of conditions and the following disclaimer in the documentation
 * // and/or other materials provided with the distribution.
 * //
 * // 3.  Neither the name of the copyright holder nor the names of its
 * // contributors may be used to endorse or promote products derived from
 * // this software without specific prior written permission.
 * //
 * // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 * // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 * // DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
 * // FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
 * // DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
 * // SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
 * // CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
 * // OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
 * // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */
use crate::common::*;
use crate::polyeval::f_polyeval3;

static LOG10_R: [u64; 128] = [
    0x0000000000000000,
    0x3f6be76bd77b4fc3,
    0x3f7c03a80ae5e054,
    0x3f851824c7587eb0,
    0x3f8c3d0837784c41,
    0x3f91b85d6044e9ae,
    0x3f9559bd2406c3ba,
    0x3f9902c31d62a843,
    0x3f9cb38fccd8bfdb,
    0x3f9e8eeb09f2f6cb,
    0x3fa125d0432ea20e,
    0x3fa30838cdc2fbfd,
    0x3fa3faf7c663060e,
    0x3fa5e3966b7e9295,
    0x3fa7d070145f4fd7,
    0x3fa8c878eeb05074,
    0x3faabbcebd84fca0,
    0x3fabb7209d1e24e5,
    0x3fadb11ed766abf4,
    0x3faeafd05035bd3b,
    0x3fb0585283764178,
    0x3fb0d966cc6500fa,
    0x3fb1dd5460c8b16f,
    0x3fb2603072a25f82,
    0x3fb367ba3aaa1883,
    0x3fb3ec6ad5407868,
    0x3fb4f7aad9bbcbaf,
    0x3fb57e3d47c3af7b,
    0x3fb605735ee985f1,
    0x3fb715d0ce367afc,
    0x3fb79efb57b0f803,
    0x3fb828cfed29a215,
    0x3fb93e7de0fc3e80,
    0x3fb9ca5aa1729f45,
    0x3fba56e8325f5c87,
    0x3fbae4285509950b,
    0x3fbb721cd17157e3,
    0x3fbc902a19e65111,
    0x3fbd204698cb42bd,
    0x3fbdb11ed766abf4,
    0x3fbe42b4c16caaf3,
    0x3fbed50a4a26eafc,
    0x3fbffbfc2bbc7803,
    0x3fc0484e4942aa43,
    0x3fc093025a19976c,
    0x3fc0de1b56356b04,
    0x3fc1299a4fb3e306,
    0x3fc175805d1587c1,
    0x3fc1c1ce9955c0c6,
    0x3fc20e8624038fed,
    0x3fc25ba8215af7fc,
    0x3fc2a935ba5f1479,
    0x3fc2f7301cf4e87b,
    0x3fc345987bfeea91,
    0x3fc394700f7953fd,
    0x3fc3e3b8149739d4,
    0x3fc43371cde076c2,
    0x3fc4839e83506c87,
    0x3fc4d43f8275a483,
    0x3fc525561e9256ee,
    0x3fc576e3b0bde0a7,
    0x3fc5c8e998072fe2,
    0x3fc61b6939983048,
    0x3fc66e6400da3f77,
    0x3fc6c1db5f9bb336,
    0x3fc6c1db5f9bb336,
    0x3fc715d0ce367afc,
    0x3fc76a45cbb7e6ff,
    0x3fc7bf3bde099f30,
    0x3fc814b4921bd52b,
    0x3fc86ab17c10bc7f,
    0x3fc86ab17c10bc7f,
    0x3fc8c13437695532,
    0x3fc9183e673394fa,
    0x3fc96fd1b639fc09,
    0x3fc9c7efd734a2f9,
    0x3fca209a84fbcff8,
    0x3fca209a84fbcff8,
    0x3fca79d382bc21d9,
    0x3fcad39c9c2c6080,
    0x3fcb2df7a5c50299,
    0x3fcb2df7a5c50299,
    0x3fcb88e67cf97980,
    0x3fcbe46b087354bc,
    0x3fcc4087384f4f80,
    0x3fcc4087384f4f80,
    0x3fcc9d3d065c5b42,
    0x3fccfa8e765cbb72,
    0x3fccfa8e765cbb72,
    0x3fcd587d96494759,
    0x3fcdb70c7e96e7f3,
    0x3fcdb70c7e96e7f3,
    0x3fce163d527e68cf,
    0x3fce76124046b3f3,
    0x3fce76124046b3f3,
    0x3fced68d819191fc,
    0x3fcf37b15bab08d1,
    0x3fcf37b15bab08d1,
    0x3fcf99801fdb749d,
    0x3fcffbfc2bbc7803,
    0x3fcffbfc2bbc7803,
    0x3fd02f93f4c87101,
    0x3fd06182e84fd4ac,
    0x3fd06182e84fd4ac,
    0x3fd093cc32c90f84,
    0x3fd093cc32c90f84,
    0x3fd0c6711d6abd7a,
    0x3fd0f972f87ff3d6,
    0x3fd0f972f87ff3d6,
    0x3fd12cd31b9c99ff,
    0x3fd12cd31b9c99ff,
    0x3fd16092e5d3a9a6,
    0x3fd194b3bdef6b9e,
    0x3fd194b3bdef6b9e,
    0x3fd1c93712abc7ff,
    0x3fd1c93712abc7ff,
    0x3fd1fe1e5af2c141,
    0x3fd1fe1e5af2c141,
    0x3fd2336b161b3337,
    0x3fd2336b161b3337,
    0x3fd2691ecc29f042,
    0x3fd2691ecc29f042,
    0x3fd29f3b0e15584b,
    0x3fd29f3b0e15584b,
    0x3fd2d5c1760b86bb,
    0x3fd2d5c1760b86bb,
    0x3fd30cb3a7bb3625,
    0x3fd34413509f79ff,
];

/// Logarithm of base 10
///
/// Max found ULP 0.5
pub fn f_log10f(x: f32) -> f32 {
    let mut x_u = x.to_bits();

    const E_BIAS: u32 = (1u32 << (8 - 1u32)) - 1u32;

    let mut m = -(E_BIAS as i32);
    if x_u < f32::MIN_POSITIVE.to_bits() || x_u > f32::MAX.to_bits() {
        if x == 0.0 {
            return f32::NEG_INFINITY;
        }
        if x_u == 0x80000000u32 {
            return f32::NEG_INFINITY;
        }
        if x.is_sign_negative() && !x.is_nan() {
            return f32::NAN + x;
        }
        // x is +inf or nan
        if x.is_nan() || x.is_infinite() {
            return x + x;
        }
        // Normalize denormal inputs.
        x_u = (x * f64::from_bits(0x4160000000000000) as f32).to_bits();
        m -= 23;
    }

    m = m.wrapping_add(x_u.wrapping_shr(23) as i32);
    let index = (x_u >> 16) & 0x7F;
    x_u = set_exponent_f32(x_u, 0x7F);

    let v;
    let u = f32::from_bits(x_u);

    #[cfg(any(
        all(
            any(target_arch = "x86", target_arch = "x86_64"),
            target_feature = "fma"
        ),
        target_arch = "aarch64"
    ))]
    {
        v = f_fmlaf(
            u,
            f32::from_bits(crate::logs::logf::LOG_REDUCTION_F32.0[index as usize]),
            -1.0,
        ) as f64; // Exact.
    }
    #[cfg(not(any(
        all(
            any(target_arch = "x86", target_arch = "x86_64"),
            target_feature = "fma"
        ),
        target_arch = "aarch64"
    )))]
    {
        v = f_fmla(
            u as f64,
            f32::from_bits(crate::logs::logf::LOG_REDUCTION_F32.0[index as usize]) as f64,
            -1.0,
        ); // Exact
    }
    // Degree-5 polynomial approximation of log10 generated by:
    // > P = fpminimax(log10(1 + x)/x, 4, [|D...|], [-2^-8, 2^-7]);
    const COEFFS: [u64; 5] = [
        0x3fdbcb7b1526e2e5,
        0xbfcbcb7b1528d43d,
        0x3fc287a77eb4ca0d,
        0xbfbbcb8110a181b5,
        0x3fb60e7e3e747129,
    ];
    let v2 = v * v; // Exact
    let p2 = f_fmla(v, f64::from_bits(COEFFS[4]), f64::from_bits(COEFFS[3]));
    let p1 = f_fmla(v, f64::from_bits(COEFFS[2]), f64::from_bits(COEFFS[1]));
    let p0 = f_fmla(
        v,
        f64::from_bits(COEFFS[0]),
        f64::from_bits(LOG10_R[index as usize]),
    );
    const LOG_10_2: f64 = f64::from_bits(0x3fd34413509f79ff);
    let r = f_fmla(m as f64, LOG_10_2, f_polyeval3(v2, p0, p1, p2));
    r as f32
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn test_log10f() {
        assert_eq!(f_log10f(0.35), -0.45593196f32);
        assert_eq!(f_log10f(0.9), -4.5757502e-2);
        assert_eq!(f_log10f(10.), 1.);
        assert_eq!(f_log10f(100.), 2.);
        assert_eq!(f_log10f(0.), f32::NEG_INFINITY);
        assert!(f_log10f(-1.).is_nan());
        assert!(f_log10f(f32::NAN).is_nan());
        assert!(f_log10f(f32::NEG_INFINITY).is_nan());
        assert_eq!(f_log10f(f32::INFINITY), f32::INFINITY);
    }
}

Coverage Report

Created: 2025-10-14 06:57

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::common::*;
30		use crate::polyeval::f_polyeval3;
31
32		static 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
166	0	pub fn f_log10f(x: f32) -> f32 {
167	0	let mut x_u = x.to_bits();
168
169		const E_BIAS: u32 = (1u32 << (8 - 1u32)) - 1u32;
170
171	0	let mut m = -(E_BIAS as i32);
172	0	if x_u < f32::MIN_POSITIVE.to_bits() \|\| x_u > f32::MAX.to_bits() {
173	0	if x == 0.0 {
174	0	return f32::NEG_INFINITY;
175	0	}
176	0	if x_u == 0x80000000u32 {
177	0	return f32::NEG_INFINITY;
178	0	}
179	0	if x.is_sign_negative() && !x.is_nan() {
180	0	return f32::NAN + x;
181	0	}
182		// x is +inf or nan
183	0	if x.is_nan() \|\| x.is_infinite() {
184	0	return x + x;
185	0	}
186		// Normalize denormal inputs.
187	0	x_u = (x * f64::from_bits(0x4160000000000000) as f32).to_bits();
188	0	m -= 23;
189	0	}
190
191	0	m = m.wrapping_add(x_u.wrapping_shr(23) as i32);
192	0	let index = (x_u >> 16) & 0x7F;
193	0	x_u = set_exponent_f32(x_u, 0x7F);
194
195		let v;
196	0	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	0	{
220	0	v = f_fmla(
221	0	u as f64,
222	0	f32::from_bits(crate::logs::logf::LOG_REDUCTION_F32.0[index as usize]) as f64,
223	0	-1.0,
224	0	); // Exact
225	0	}
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	0	let v2 = v * v; // Exact
236	0	let p2 = f_fmla(v, f64::from_bits(COEFFS[4]), f64::from_bits(COEFFS[3]));
237	0	let p1 = f_fmla(v, f64::from_bits(COEFFS[2]), f64::from_bits(COEFFS[1]));
238	0	let p0 = f_fmla(
239	0	v,
240	0	f64::from_bits(COEFFS[0]),
241	0	f64::from_bits(LOG10_R[index as usize]),
242		);
243		const LOG_10_2: f64 = f64::from_bits(0x3fd34413509f79ff);
244	0	let r = f_fmla(m as f64, LOG_10_2, f_polyeval3(v2, p0, p1, p2));
245	0	r as f32
246	0	}
247
248		#[cfg(test)]
249		mod 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		}