1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
use super::{EXPONENT_MASK, F32};
use core::f32::consts::LN_2;
impl F32 {
#[allow(clippy::excessive_precision)]
pub fn ln(self) -> Self {
if (self - Self::ONE).abs() < f32::EPSILON {
return Self::ZERO;
}
let x_less_than_1 = self < 1.0;
let x_working = if x_less_than_1 { self.inv() } else { self };
let base2_exponent = x_working.extract_exponent_value() as u32;
let divisor = f32::from_bits(x_working.to_bits() & EXPONENT_MASK);
let x_working = x_working / divisor;
let ln_1to2_polynomial = -1.741_793_9
+ (2.821_202_6
+ (-1.469_956_8 + (0.447_179_55 - 0.056_570_851 * x_working) * x_working)
* x_working)
* x_working;
let result = (base2_exponent as f32) * LN_2 + ln_1to2_polynomial;
if x_less_than_1 {
-result
} else {
result
}
}
}
#[cfg(test)]
mod tests {
use super::F32;
pub(crate) const MAX_ERROR: f32 = 0.001;
pub(crate) const TEST_VECTORS: &[(f32, f32)] = &[
(1e-20, -46.0517),
(1e-19, -43.749115),
(1e-18, -41.446533),
(1e-17, -39.143948),
(1e-16, -36.841362),
(1e-15, -34.538776),
(1e-14, -32.23619),
(1e-13, -29.933607),
(1e-12, -27.631021),
(1e-11, -25.328436),
(1e-10, -23.02585),
(1e-09, -20.723267),
(1e-08, -18.420681),
(1e-07, -16.118095),
(1e-06, -13.815511),
(1e-05, -11.512925),
(1e-04, -9.2103405),
(0.001, -6.9077554),
(0.01, -4.6051702),
(0.1, -2.3025851),
(10.0, 2.3025851),
(100.0, 4.6051702),
(1000.0, 6.9077554),
(10000.0, 9.2103405),
(100000.0, 11.512925),
(1000000.0, 13.815511),
(10000000.0, 16.118095),
(100000000.0, 18.420681),
(1000000000.0, 20.723267),
(10000000000.0, 23.02585),
(100000000000.0, 25.328436),
(1000000000000.0, 27.631021),
(10000000000000.0, 29.933607),
(100000000000000.0, 32.23619),
(1000000000000000.0, 34.538776),
(1e+16, 36.841362),
(1e+17, 39.143948),
(1e+18, 41.446533),
(1e+19, 43.749115),
];
#[test]
fn sanity_check() {
assert_eq!(F32::ONE.ln(), F32::ZERO);
for &(x, expected) in TEST_VECTORS {
let ln_x = F32(x).ln();
let relative_error = (ln_x - expected).abs() / expected;
assert!(
relative_error <= MAX_ERROR,
"relative_error {} too large: {} vs {}",
relative_error,
ln_x,
expected
);
}
}
}