1//! Natural log (ln) approximation for a single-precision float.
2//!
3//! Method described at: <https://stackoverflow.com/a/44232045/>
4//!
5//! Modified to not be restricted to int range and only values of x above 1.0.
6//! Also got rid of most of the slow conversions. Should work for all positive values of x.
7
8use super::{EXPONENT_MASK, F32};
9use core::f32::consts::LN_2;
10
11impl F32 {
12 /// Approximates the natural logarithm of the number.
13 // Note: excessive precision ignored because it hides the origin of the numbers used for the
14 // ln(1.0->2.0) polynomial
15 #[allow(clippy::excessive_precision)]
16 pub fn ln(self) -> Self {
17 // x may essentially be 1.0 but, as clippy notes, these kinds of
18 // floating point comparisons can fail when the bit pattern is not the sames
19 if (self - Self::ONE).abs() < f32::EPSILON {
20 return Self::ZERO;
21 }
22
23 let x_less_than_1 = self < 1.0;
24
25 // Note: we could use the fast inverse approximation here found in super::inv::inv_approx, but
26 // the precision of such an approximation is assumed not good enough.
27 let x_working = if x_less_than_1 { self.inv() } else { self };
28
29 // according to the SO post ln(x) = ln((2^n)*y)= ln(2^n) + ln(y) = ln(2) * n + ln(y)
30 // get exponent value
31 let base2_exponent = x_working.extract_exponent_value() as u32;
32 let divisor = f32::from_bits(x_working.to_bits() & EXPONENT_MASK);
33
34 // supposedly normalizing between 1.0 and 2.0
35 let x_working = x_working / divisor;
36
37 // approximate polynomial generated from maple in the post using Remez Algorithm:
38 // https://en.wikipedia.org/wiki/Remez_algorithm
39 let ln_1to2_polynomial = -1.741_793_9
40 + (2.821_202_6
41 + (-1.469_956_8 + (0.447_179_55 - 0.056_570_851 * x_working) * x_working)
42 * x_working)
43 * x_working;
44
45 // ln(2) * n + ln(y)
46 let result = (base2_exponent as f32) * LN_2 + ln_1to2_polynomial;
47
48 if x_less_than_1 {
49 -result
50 } else {
51 result
52 }
53 }
54}
55
56#[cfg(test)]
57mod tests {
58 use super::F32;
59
60 pub(crate) const MAX_ERROR: f32 = 0.001;
61
62 /// ln(x) test vectors - `(input, output)`
63 pub(crate) const TEST_VECTORS: &[(f32, f32)] = &[
64 (1e-20, -46.0517),
65 (1e-19, -43.749115),
66 (1e-18, -41.446533),
67 (1e-17, -39.143948),
68 (1e-16, -36.841362),
69 (1e-15, -34.538776),
70 (1e-14, -32.23619),
71 (1e-13, -29.933607),
72 (1e-12, -27.631021),
73 (1e-11, -25.328436),
74 (1e-10, -23.02585),
75 (1e-09, -20.723267),
76 (1e-08, -18.420681),
77 (1e-07, -16.118095),
78 (1e-06, -13.815511),
79 (1e-05, -11.512925),
80 (1e-04, -9.2103405),
81 (0.001, -6.9077554),
82 (0.01, -4.6051702),
83 (0.1, -2.3025851),
84 (10.0, 2.3025851),
85 (100.0, 4.6051702),
86 (1000.0, 6.9077554),
87 (10000.0, 9.2103405),
88 (100000.0, 11.512925),
89 (1000000.0, 13.815511),
90 (10000000.0, 16.118095),
91 (100000000.0, 18.420681),
92 (1000000000.0, 20.723267),
93 (10000000000.0, 23.02585),
94 (100000000000.0, 25.328436),
95 (1000000000000.0, 27.631021),
96 (10000000000000.0, 29.933607),
97 (100000000000000.0, 32.23619),
98 (1000000000000000.0, 34.538776),
99 (1e+16, 36.841362),
100 (1e+17, 39.143948),
101 (1e+18, 41.446533),
102 (1e+19, 43.749115),
103 ];
104
105 #[test]
106 fn sanity_check() {
107 assert_eq!(F32::ONE.ln(), F32::ZERO);
108 for &(x, expected) in TEST_VECTORS {
109 let ln_x = F32(x).ln();
110 let relative_error = (ln_x - expected).abs() / expected;
111
112 assert!(
113 relative_error <= MAX_ERROR,
114 "relative_error {} too large: {} vs {}",
115 relative_error,
116 ln_x,
117 expected
118 );
119 }
120 }
121}
122