fastmaths 0.1.0

High-performance no_std math routines targeting glibc parity and strict f64 accuracy.
Documentation 
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

//! atan2(y,x) implementation.
//!
//! Handles quadrants, signed zeros, and infinities explicitly, then reduces to
//! atan(y/x) using the atan polynomial core. Matches glibc sign/edge behavior.

use super::{PIO2_HI, atan, fma_internal};

const PI: f64 = core::f64::consts::PI;
const PI_LO: f64 = 1.224_646_799_147_353_207_2e-16;
const PIO4: f64 = core::f64::consts::FRAC_PI_4;

#[inline(always)]
fn div_hi_lo(n: f64, d: f64) -> (f64, f64) {
    let r0 = n / d;
    if !r0.is_finite() {
        return (r0, 0.0);
    }
    if n.is_subnormal() || d.is_subnormal() {
        return (r0, 0.0);
    }
    let err = fma_internal(-r0, d, n);
    (r0, err / d)
}

#[inline(always)]
fn atan_with_correction(r: f64, r_lo: f64) -> f64 {
    let z = atan(r);
    if !r.is_finite() {
        return z;
    }
    if r_lo == 0.0 {
        return z;
    }
    let denom = 1.0 + r * r;
    z + r_lo / denom
}

#[inline]
pub fn atan2(y: f64, x: f64) -> f64 {
    if y.is_nan() || x.is_nan() {
        return f64::NAN;
    }

    if x.is_infinite() {
        if y.is_infinite() {
            return if x.is_sign_positive() {
                if y.is_sign_negative() { -PIO4 } else { PIO4 }
            } else if y.is_sign_negative() {
                -3.0 * PIO4
            } else {
                3.0 * PIO4
            };
        }
        return if x.is_sign_positive() {
            0.0_f64.copysign(y)
        } else if y.is_sign_negative() {
            -PI
        } else {
            PI
        };
    }

    if y.is_infinite() {
        return if y.is_sign_negative() {
            -PIO2_HI
        } else {
            PIO2_HI
        };
    }

    if x == 0.0 {
        if y == 0.0 {
            if x.is_sign_negative() {
                return if y.is_sign_negative() { -PI } else { PI };
            }
            return y;
        }
        return if y.is_sign_negative() {
            -PIO2_HI
        } else {
            PIO2_HI
        };
    }

    if y == 0.0 {
        if x.is_sign_positive() {
            return y;
        }
        return if y.is_sign_negative() { -PI } else { PI };
    }

    if x.is_sign_positive() {
        let (r, r_lo) = div_hi_lo(y, x);
        return atan_with_correction(r, r_lo);
    }

    let (r, r_lo) = div_hi_lo(y, x);
    let z = atan_with_correction(r, r_lo);
    if y.is_sign_negative() {
        (z - PI_LO) - PI
    } else {
        (z + PI_LO) + PI
    }
}