fast_math/
log.rs

1use core::f32 as f;
2use ieee754::Ieee754;
3
4/// Compute a fast approximation of the base-2 logarithm of `x`.
5///
6/// The maximum relative error across all positive f32s (including
7/// denormals) is less than 0.022. The maximum absolute error is less
8/// than 0.009.
9///
10/// If `x` is negative, or NaN, `log2` returns `NaN`.
11///
12/// See also `log2_raw` which only works on positive, finite,
13/// non-denormal floats, but is 30-40% faster.
14///
15///
16/// |               | Time (ns) |
17/// |--------------:|----------------|
18/// |    `x.log2()` | 14.3           |
19/// |     `log2(x)` | 4.0            |
20/// | `log2_raw(x)` | 2.7            |
21#[inline]
22pub fn log2(x: f32) -> f32 {
23    let (sign, exp, signif) = x.decompose_raw();
24    if sign {
25        f::NAN
26    } else if exp == 0 {
27        log2_exp_0(signif)
28    } else if exp == 0xFF {
29        if signif == 0 {
30            f::INFINITY
31        } else {
32            f::NAN
33        }
34    } else {
35        log2_raw(x)
36    }
37}
38
39#[inline(never)]
40fn log2_exp_0(signif: u32) -> f32 {
41    if signif == 0 {
42        f::NEG_INFINITY
43    } else {
44        // denormal
45        let zeros = signif.leading_zeros() - 9 + 1;
46        -126.0 - zeros as f32 + log2(f32::recompose_raw(false, 127, signif << zeros))
47    }
48}
49
50/// Compute a fast approximation of the base-2 logarithm of **positive,
51/// finite, non-denormal** `x`.
52///
53/// This will return unspecified nonsense if `x` is doesn't not
54/// satisfy those constraints. Use `log2` if correct handling is
55/// required (at the expense of some speed).
56///
57/// The maximum relative error across all valid input is less than
58/// 0.022. The maximum absolute error is less than 0.009.
59///
60/// |               | Time (ns) |
61/// |--------------:|----------------|
62/// |    `x.log2()` | 14.3           |
63/// |     `log2(x)` | 4.0            |
64/// | `log2_raw(x)` | 2.7            |
65#[inline]
66pub fn log2_raw(x: f32) -> f32 {
67    let (_sign, exp, signif) = x.decompose_raw();
68    debug_assert!(!_sign && 1 <= exp && exp <= 254);
69
70    let high_bit = ((signif >> 22) & 1) as u8;
71    let add_exp = (exp + high_bit) as i32 - 127;
72    let normalised = f32::recompose_raw(false, 0x7F ^ high_bit, signif) - 1.0;
73    const A: f32 = -0.6296735;
74    const B: f32 = 1.466967;
75    add_exp as f32 + normalised * (B + A * normalised)
76}
77
78#[cfg(test)]
79mod tests {
80    use super::*;
81    use quickcheck as qc;
82    use std::f32 as f;
83    use ieee754::Ieee754;
84
85    #[test]
86    fn log2_rel_err_qc() {
87        fn prop(x: f32) -> qc::TestResult {
88            if !(x > 0.0) { return qc::TestResult::discard() }
89
90            let e = log2(x);
91            let t = x.log2();
92
93            qc::TestResult::from_bool(e.rel_error(t).abs() < 0.025)
94        }
95        qc::quickcheck(prop as fn(f32) -> qc::TestResult)
96    }
97    const PREC: u32 = 1 << 20;
98    #[test]
99    fn log2_rel_err_exhaustive() {
100        let mut max = 0.0;
101        for i in 0..PREC + 1 {
102            for j in -5..6 {
103                let x = (1.0 + i as f32 / PREC as f32) * 2f32.powi(j * 20);
104                let e = log2(x);
105                let t = x.log2();
106                let rel = e.rel_error(t).abs();
107                if rel > max { max = rel }
108                assert!(rel < 0.025 && (e - t).abs() < 0.009,
109                        "{:.8}: {:.8}, {:.8}. {:.4}", x, e, t, rel);
110            }
111        }
112        println!("maximum {}", max);
113    }
114
115    #[test]
116    fn edge_cases() {
117        assert!(log2(f::NAN).is_nan());
118        assert!(log2(-1.0).is_nan());
119        assert!(log2(f::NEG_INFINITY).is_nan());
120        assert_eq!(log2(f::INFINITY), f::INFINITY);
121        assert_eq!(log2(0.0), f::NEG_INFINITY);
122        assert_eq!(log2(f32::recompose_raw(false, 0, 1)), -149.0);
123    }
124
125    #[test]
126    fn denormals() {
127        fn prop(x: u8, y: u16) -> bool {
128            let signif = ((x as u32) << 16) | (y as u32);
129            let mut x = f32::recompose_raw(false, 1, signif);
130
131            for _ in 0..23 {
132                assert!(x > 0.0);
133                let log = x.log2();
134                let e = log2(x);
135                let rel = e.rel_error(log).abs();
136                if rel >= 0.025 || (e - log).abs() > 0.009 {
137                    return false
138                }
139
140                x /= 2.0;
141            }
142            true
143        }
144        qc::quickcheck(prop as fn(u8, u16) -> bool)
145    }
146}
fast_math/log.rs

fast_math/
log.rs
