fast_math/
atan.rs

1use core::f32::INFINITY;
2use core::f32::consts::{PI, FRAC_PI_2, FRAC_PI_4};
3use float::{flip_sign_nonnan};
4use ieee754::Ieee754;
5
6/// Compute a fast approximation of the inverse tangent for `|x| < 1`.
7///
8/// This will return unspecified nonsense if `x` is doesn't not
9/// satisfy those constraints. Use `atan2` if correct handling is
10/// required (at the expense of some speed).
11#[inline]
12pub fn atan_raw(x: f32) -> f32 {
13    // Quadratic approximation recommended in
14    // http://www-labs.iro.umontreal.ca/~mignotte/IFT2425/Documents/EfficientApproximationArctgFunction.pdf.
15    const N2: f32 = 0.273;
16    (FRAC_PI_4 + N2 - N2 * x.abs()) * x
17}
18
19/// Compute a fast approximation of the arctangent of `x`.
20///
21/// The maximum absolute error across all f32s is less than 0.0038.
22///
23/// See also `atan_raw` which only works on `|x| <= 1`, but is faster.
24#[inline]
25pub fn atan(x: f32) -> f32 {
26    if x.abs() > 1.0 {
27        // if x is NaN, abs(x) is NaN, so the comparison can't succeed
28        debug_assert!(!x.is_nan());
29        flip_sign_nonnan(x, FRAC_PI_2) - atan_raw(1./x)
30    } else {
31        atan_raw(x)
32    }
33}
34
35/// Compute a fast approximation of the four quadrant arctangent of `y` and `x`.
36///
37/// The maximum absolute error across all f32s is less than 0.0038.
38#[inline]
39pub fn atan2(y: f32, x: f32) -> f32 {
40    if y.abs() < x.abs() {
41        // x is not NaN and y is finite, so there should be no NaNs
42        // around
43        debug_assert!(!x.is_nan() && !y.is_nan() && !(y / x).is_nan());
44
45        let bias = if x > 0.0 { 0.0 } else { PI };
46        flip_sign_nonnan(y, bias) + atan_raw(y / x)
47    } else if x == 0. {
48        // x is non-NaN
49        if y == 0. {
50            let bias = if x.is_sign_positive() { 0.0 } else { PI };
51            flip_sign_nonnan(y, bias)
52        } else if y.is_nan() {
53            y
54        } else {
55            FRAC_PI_2.copy_sign(y)
56        }
57    } else if y.abs() == INFINITY && x.abs() == INFINITY {
58        // x and y are both infinite, meaning: not NaN, can't be
59        // divided, and the answer is statically obvious (some
60        // multiple of PI/4).
61        flip_sign_nonnan(y, FRAC_PI_2 - flip_sign_nonnan(x, FRAC_PI_4))
62    } else {
63        // Either one x or y is NaN (propogates through atan_raw
64        // properly), or |y| >= |x| (meaning |r| = |y / x| >= 1). Use
65        // `atan(1/r) == sign(r) * pi / 2 - atan(r)`, but inline the 0
66        // or PI `x` bias.
67        flip_sign_nonnan(y, FRAC_PI_2) - atan_raw(x / y)
68    }
69}
70
71#[cfg(test)]
72mod tests {
73    use super::*;
74    use quickcheck as qc;
75    use std::f32 as f;
76    use ieee754::Ieee754;
77
78    /// Maximal absolute error according to paper.
79    const TOL: f32 = 0.0038;
80
81    #[test]
82    fn atan_abs_err_qc() {
83        fn prop(x: f32) -> qc::TestResult {
84            let e = atan(x);
85            let t = x.atan();
86            let abs = (e - t).abs();
87
88            if x == 0.0 {
89                qc::TestResult::from_bool(e == 0.0)
90            } else {
91                qc::TestResult::from_bool(abs < TOL)
92            }
93        }
94        qc::quickcheck(prop as fn(f32) -> qc::TestResult)
95    }
96
97    const PREC: u32 = 1 << 20;
98    #[test]
99    fn atan_abs_err_exhaustive() {
100        for i in 0..PREC + 1 {
101            for j in -5..6 {
102                let x = (1.0 + i as f32 / PREC as f32) * 2f32.powi(j * 20);
103                let e = atan(x);
104                let t = x.atan();
105                let abs = (e - t).abs();
106
107                assert!((e == 0. && x == 0.) || abs < TOL,
108                        "{:.8}: {:.8}, {:.8}. {:.4}", x, e, t, abs);
109            }
110        }
111    }
112
113    #[test]
114    fn atan_edge_cases() {
115        assert!(atan(f::NAN).is_nan());
116        assert_eq!(atan(f::NEG_INFINITY), -PI / 2.);
117        assert_eq!(atan(0.), 0.);
118        assert_eq!(atan(f::INFINITY), PI / 2.);
119    }
120
121    #[test]
122    fn atan_denormals() {
123        fn prop(x: u8, y: u16) -> bool {
124            let signif = ((x as u32) << 16) | (y as u32);
125            let mut x = f32::recompose_raw(false, 1, signif);
126
127            for _ in 0..23 {
128                assert!(x > 0.0);
129                let e = atan(x);
130                let t = x.atan();
131                let abs = (e - t).abs();
132                if abs >= TOL {
133                    return false
134                }
135
136                x /= 2.0;
137            }
138            true
139        }
140        qc::quickcheck(prop as fn(u8, u16) -> bool)
141    }
142
143    #[test]
144    fn atan2_abs_err_qc() {
145        fn prop(y: f32, x: f32) -> qc::TestResult {
146            let e = atan2(y, x);
147            let t = y.atan2(x);
148            let abs = (e - t).abs();
149
150            qc::TestResult::from_bool(abs < 0.0038)
151        }
152        qc::quickcheck(prop as fn(f32, f32) -> qc::TestResult)
153    }
154
155    #[test]
156    fn atan2_edge_cases() {
157        let values = &[-2., -1., -0., 0., 1., 2., f::INFINITY, f::NEG_INFINITY, f::NAN];
158        for &x in values {
159            for &y in values {
160                let e = atan2(x, y);
161                let t = x.atan2(y);
162                assert_eq!(e.is_nan(), t.is_nan());
163                if !t.is_nan() {
164                    assert!((e - t).abs() < TOL ||
165                            (e - t - 2.0 * PI).abs() < TOL ||
166                            (e - t + 2.0 * PI).abs() < TOL);
167                }
168            }
169        }
170    }
171}
fast_math/atan.rs

fast_math/
atan.rs
