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
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
// SPDX-License-Identifier: MIT
// Copyright 2023 IROX Contributors
//

//!
//! A collection of utilities for the f32 built-in
//!

impl crate::f64::FloatExt for f32 {
    type Type = f32;

    fn trunc(self) -> f32 {
        (self as u64) as f32
    }

    fn fract(self) -> f32 {
        self - self.trunc()
    }

    fn abs(self) -> f32 {
        f32::from_bits(self.to_bits() & 0x7FFF_FFFF)
    }
    fn round(self) -> f32 {
        (self + 0.5 * self.signum()).trunc()
    }

    fn floor(self) -> f32 {
        if self.is_sign_negative() {
            return (self - 1.0).trunc();
        }
        self.trunc()
    }

    fn ceil(self) -> f32 {
        if self.is_sign_positive() {
            return (self + 1.0).trunc();
        }
        self.trunc()
    }

    fn signum(self) -> f32 {
        if self.is_nan() {
            return f32::NAN;
        }
        if self.is_sign_negative() {
            return -1.0;
        }
        1.0
    }
    ///
    /// Implementation of Exponential Function from NIST DTMF eq 4.2.19: `<https://dlmf.nist.gov/4.2.E19>`
    fn exp(self) -> Self::Type {
        let mut out = 1.0;
        let i = self;
        let mut z = self;
        let mut exp = 1.0;
        let mut idx = 1;
        let mut next = self;

        while next.abs() > f32::EPSILON {
            out += next;
            idx += 1;
            z *= i;
            if z.is_infinite() {
                break;
            }
            exp *= idx as Self::Type;
            if exp.is_infinite() {
                break;
            }
            next = z / exp;
            if next.is_infinite() {
                break;
            }
        }

        out
    }

    ///
    /// Implementation of Natural Logarithm using NIST DLMF eq 4.6.4: `<https://dlmf.nist.gov/4.6.E4>`
    fn ln(self) -> Self::Type {
        let z = self as f64;
        let iter = (z - 1.) / (z + 1.);
        let mut out = 0.0f64;
        let mut next = iter;
        let mut base = iter;
        let mut idx = 1u64;
        while next.abs() > f64::EPSILON {
            out += next;
            idx += 2;
            base *= iter * iter;
            next = base / idx as f64;
        }
        (out * 2.0) as f32
    }

    ///
    /// Implementation of general power function using NIST DLMF eq 4.2.26: `<https://dlmf.nist.gov/4.2.E26>`
    fn powf(self, a: Self::Type) -> Self::Type {
        let z = self;

        (a * z.ln()).exp()
    }

    /// Naive implementation of integer power fn.  Will do something smarter later.
    fn powi(self, val: i32) -> Self::Type {
        let mut out = self;
        let i = self;
        for _ in 0..val.abs() {
            out *= i;
        }
        out
    }

    fn sqrt(self) -> Self::Type {
        self.powf(0.5)
    }
}

#[cfg(all(test, not(feature = "std")))]
mod tests {
    #[test]
    pub fn test_ln() {
        assert_eq!(0.0, crate::f64::FloatExt::ln(1.0_f32));
        assert_eq_eps!(1.0, crate::f64::FloatExt::ln(core::f32::consts::E), 1e-6);
        assert_eq_eps!(4.6051702, crate::f64::FloatExt::ln(100f32), 1e-6);
        assert_eq_eps!(
            11.09033963004403,
            crate::f64::FloatExt::ln(u16::MAX as f32),
            1e-6
        );
    }

    #[test]
    pub fn test_exp() {
        assert_eq_eps!(1.0, crate::f64::FloatExt::exp(0.0f32), 1e-6);
        assert_eq_eps!(
            core::f32::consts::E,
            crate::f64::FloatExt::exp(1.0f32),
            1e-6
        );
        assert_eq_eps!(7.389056098930649, crate::f64::FloatExt::exp(2.0f32), 1e-6);
        assert_eq_eps!(
            15.154261,
            crate::f64::FloatExt::exp(core::f32::consts::E),
            1e-6
        );
    }
}