fast_math/
exp.rs

1use core::f32;
2use core::f32::consts as f;
3use float;
4use ieee754::Ieee754;
5
6#[derive(Clone, Copy)]
7enum Base {
8    E,
9    Two,
10}
11impl Base {
12    #[inline(always)]
13    fn log2(self) -> f32 {
14        match self {
15            Base::E => f::LOG2_E,
16            Base::Two => 1.0,
17        }
18    }
19
20    #[inline(always)]
21    fn upper_limit(self) -> f32 {
22        128.0 / self.log2()
23    }
24
25    #[inline(always)]
26    fn lower_limit(self) -> f32 {
27        -127.0 / self.log2()
28    }
29}
30
31#[inline(always)]
32fn exp_raw_impl(x: f32, base: Base) -> f32 {
33    const A: f32 = (1 << float::SIGNIF) as f32;
34    const MASK: i32 = 0xff800000u32 as i32;
35    const EXP2_23: f32 = 1.1920929e-7;
36    const C0: f32 = 0.3371894346 * EXP2_23 * EXP2_23;
37    const C1: f32 = 0.657636276 * EXP2_23;
38    const C2: f32 = 1.00172476;
39
40    let a = A * base.log2();
41    let mul = (a * x) as i32;
42    let floor = mul & MASK;
43    let frac = (mul - floor) as f32;
44
45    let approx = (C0 * frac + C1) * frac + C2;
46    f32::from_bits(approx.bits().wrapping_add(floor as u32))
47}
48
49#[inline(always)]
50fn exp_impl(x: f32, base: Base) -> f32 {
51    if x <= base.lower_limit() {
52        0.0
53    } else if x < base.upper_limit() {
54        exp_raw_impl(x, base)
55    } else {
56        // too big, or NaN, so lets overflow to infinity with some
57        // arithmetic to propagate the NaN.
58        x + f32::INFINITY
59    }
60}
61
62/// Compute a fast approximation to 2<sup><code>x</code></sup> for
63/// -151 &le; `x` &le; 151.
64///
65/// This will return unspecified nonsense if `x` does not satisfy
66/// those requirements. Use `exp2` if correct handling is required (at
67/// the expense of some speed).
68///
69/// The maximum relative error for inputs for which the result is
70/// normal (`x` &ge; -128) is less than 0.011. For `x` < -128, the
71/// relative error in the (subnormal) result can be as large as 1.
72#[inline]
73pub fn exp2_raw(x: f32) -> f32 {
74    exp_raw_impl(x, Base::Two)
75}
76
77/// Compute a fast approximation to 2<sup><code>x</code></sup>.
78///
79/// The maximum relative error for inputs for which the result is
80/// normal (`x` &ge; -128) is less than 0.011. For `x` < -128, the
81/// relative error in the (subnormal) result can be as large as 1.
82///
83/// If `x` is NaN, `exp2` returns NaN.
84///
85/// See also `exp2_raw` which only works on -151 &le; `x` &le; 151,
86/// but is % faster.
87#[inline]
88pub fn exp2(x: f32) -> f32 {
89    exp_impl(x, Base::Two)
90}
91
92/// Compute a fast approximation to *e*<sup><code>x</code></sup> for
93/// -104 &le; `x` &le; 104.
94///
95/// This will return unspecified nonsense if `x` does not satisfy
96/// those requirements. Use `exp` if correct handling is required (at
97/// the expense of some speed).
98///
99/// The maximum relative error for inputs for which the result is
100/// normal (`x` &ge; -128 ln(2) &approx; -88.7) is less than
101/// 0.011. For `x` < -128 ln(2), the relative error in the (subnormal)
102/// result can be as large as 1.
103#[inline]
104pub fn exp_raw(x: f32) -> f32 {
105    exp_raw_impl(x, Base::E)
106}
107
108/// Compute a fast approximation to *e*<sup><code>x</code></sup>.
109///
110/// The maximum relative error for inputs for which the result is
111/// normal (`x` &ge; -128 ln 2 &approx; -88.7) is less than
112/// 0.011. For `x` < -128 ln 2, the relative error in the (subnormal)
113/// result can be as large as 1.
114///
115/// If `x` is NaN, `exp` returns NaN.
116///
117/// See also `exp_raw` which only works on -104 &le; `x` &le; 104,
118/// but is % faster.
119#[inline]
120pub fn exp(x: f32) -> f32 {
121    exp_impl(x, Base::E)
122}
123
124#[cfg(test)]
125mod tests {
126    use super::*;
127    use std::{f32, num};
128
129    const PREC: u32 = 1 << 19;
130
131    #[test]
132    fn exp_rel_err_exhaustive() {
133        let mut max = 0.0;
134        for i in 0..PREC + 1 {
135            for j in -5..6 {
136                for &sign in &[-1.0, 1.0] {
137                    let x = sign * (1.0 + i as f32 / PREC as f32) * 2f32.powi(j * 2);
138                    let e = exp(x);
139                    let t = x.exp();
140                    let rel = e.rel_error(t).abs();
141
142                    if t.classify() == num::FpCategory::Subnormal {
143                        // subnormal should be approximately right
144                        assert!(rel <= 1.0,
145                                "{:.8}: e = {:.8e}, t = {:.8e}. {:.4}", x, e, t, rel);                  } else {
146                        if rel > max { max = rel }
147                        // e == t handles the infinity case
148                        assert!(rel <= 0.002,
149                                "{:.8}: e = {:.8e}, t = {:.8e}. {:.4}", x, e, t, rel);
150                    }
151                }
152            }
153        }
154        println!("maximum {}", max);
155    }
156
157    #[test]
158    fn exp2_rel_err_exhaustive() {
159        let mut max = 0.0;
160        for i in 0..PREC + 1 {
161            for j in -5..6 {
162                for &sign in &[-1.0, 1.0] {
163                    let x = sign * (1.0 + i as f32 / PREC as f32) * 2f32.powi(j * 2);
164                    let e = exp2(x);
165                    let t = x.exp2();
166                    let rel = e.rel_error(t).abs();
167                    if t.classify() == num::FpCategory::Subnormal {
168                        // subnormal should be approximately right
169                        assert!(rel <= 1.0,
170                                "{:.8}: e = {:.8e}, t = {:.8e}. {:.4}", x, e, t, rel);                  } else {
171                        if rel > max { max = rel }
172                        // e == t handles the infinity case
173                        assert!(rel <= 0.002,
174                                "{:.8}: e = {:.8e}, t = {:.8e}. {:.4}", x, e, t, rel);
175                    }
176                }
177            }
178        }
179        println!("maximum {}", max);
180    }
181
182    #[test]
183    fn exp_edge_cases() {
184        assert!(exp(f32::NAN).is_nan());
185        assert_eq!(exp(f32::NEG_INFINITY), 0.0);
186        assert!((exp(0.0) - 1.0).abs() < 0.002);
187        assert_eq!(exp(f32::INFINITY), f32::INFINITY);
188    }
189
190    #[test]
191    fn exp2_edge_cases() {
192        assert!(exp2(f32::NAN).is_nan());
193        assert_eq!(exp2(f32::NEG_INFINITY), 0.0);
194        assert!((exp2(0.0) - 1.0).abs() < 0.002);
195        assert_eq!(exp2(f32::INFINITY), f32::INFINITY);
196    }
197}
fast_math/exp.rs

fast_math/
exp.rs
