1//! Exp approximation for a single-precision float.
2//!
3//! Method described at: <https://stackoverflow.com/a/6985769/2036035>
4
5use super::{EXPONENT_BIAS, F32};
6use core::f32::consts;
7
8impl F32 {
9 /// Returns `e^(self)`, (the exponential function).
10 #[inline]
11 pub fn exp(self) -> Self {
12 self.exp_ln2_approx(4)
13 }
14
15 /// Exp approximation for `f32`.
16 pub(crate) fn exp_ln2_approx(self, partial_iter: u32) -> Self {
17 if self == Self::ZERO {
18 return Self::ONE;
19 }
20
21 if (self - Self::ONE).abs() < f32::EPSILON {
22 return consts::E.into();
23 }
24
25 if (self - (-Self::ONE)).abs() < f32::EPSILON {
26 return Self::ONE / consts::E;
27 }
28
29 // log base 2(E) == 1/ln(2)
30 // x_fract + x_whole = x/ln2_recip
31 // ln2*(x_fract + x_whole) = x
32 let x_ln2recip = self * consts::LOG2_E;
33 let x_fract = x_ln2recip.fract();
34 let x_trunc = x_ln2recip.trunc();
35
36 //guaranteed to be 0 < x < 1.0
37 let x_fract = x_fract * consts::LN_2;
38 let fract_exp = x_fract.exp_smallx(partial_iter);
39
40 //need the 2^n portion, we can just extract that from the whole number exp portion
41 let fract_exponent: i32 = fract_exp
42 .extract_exponent_value()
43 .saturating_add(x_trunc.0 as i32);
44
45 if fract_exponent < -(EXPONENT_BIAS as i32) {
46 return Self::ZERO;
47 }
48
49 if fract_exponent > ((EXPONENT_BIAS + 1) as i32) {
50 return Self::INFINITY;
51 }
52
53 fract_exp.set_exponent(fract_exponent)
54 }
55
56 /// if x is between 0.0 and 1.0, we can approximate it with the a series
57 ///
58 /// Series from here:
59 /// <https://stackoverflow.com/a/6984495>
60 ///
61 /// e^x ~= 1 + x(1 + x/2(1 + (x?
62 #[inline]
63 pub(crate) fn exp_smallx(self, iter: u32) -> Self {
64 let mut total = 1.0;
65
66 for i in (1..=iter).rev() {
67 total = 1.0 + ((self.0 / (i as f32)) * total);
68 }
69
70 Self(total)
71 }
72}
73
74#[cfg(test)]
75mod tests {
76 use super::F32;
77
78 pub(crate) const MAX_ERROR: f32 = 0.001;
79
80 /// exp test vectors - `(input, output)`
81 pub(crate) const TEST_VECTORS: &[(f32, f32)] = &[
82 (1e-07, 1.0000001),
83 (1e-06, 1.000001),
84 (1e-05, 1.00001),
85 (1e-04, 1.0001),
86 (0.001, 1.0010005),
87 (0.01, 1.0100502),
88 (0.1, 1.105171),
89 (1.0, 2.7182817),
90 (10.0, 22026.465),
91 (-1e-08, 1.0),
92 (-1e-07, 0.9999999),
93 (-1e-06, 0.999999),
94 (-1e-05, 0.99999),
95 (-1e-04, 0.9999),
96 (-0.001, 0.9990005),
97 (-0.01, 0.99004984),
98 (-0.1, 0.9048374),
99 (-1.0, 0.36787945),
100 (-10.0, 4.539_993e-5),
101 ];
102
103 #[test]
104 fn sanity_check() {
105 assert_eq!(F32(-1000000.0).exp(), F32::ZERO);
106
107 for &(x, expected) in TEST_VECTORS {
108 let exp_x = F32(x).exp();
109 let relative_error = (exp_x - expected).abs() / expected;
110
111 assert!(
112 relative_error <= MAX_ERROR,
113 "relative_error {} too large for input {} : {} vs {}",
114 relative_error,
115 x,
116 exp_x,
117 expected
118 );
119 }
120 }
121}
122