special/
error.rs

1use crate::primitive::Primitive;
2
3/// Error functions.
4pub trait Error {
5    /// Compute the error function.
6    fn error(self) -> Self;
7
8    /// Compute the complementary error function.
9    fn compl_error(self) -> Self;
10
11    /// Compute the inverse of the error function.
12    ///
13    /// The code is based on a [C implementation][1] by Alijah Ahmed.
14    ///
15    /// [1]: https://scistatcalc.blogspot.com/2013/09/numerical-estimate-of-inverse-error.html
16    fn inv_error(self) -> Self;
17}
18
19macro_rules! implement {
20    ($kind:ty) => {
21        impl Error for $kind {
22            #[inline]
23            fn error(self) -> Self {
24                self.erf()
25            }
26
27            #[inline]
28            fn compl_error(self) -> Self {
29                self.erfc()
30            }
31
32            fn inv_error(self) -> Self {
33                const SQRT_PI: $kind = 1.772453850905515881919427556567825376987457275391;
34
35                let mut w: $kind = -((1.0 - self) * (1.0 + self)).ln();
36                let mut p: $kind;
37
38                if w < 5.0 {
39                    w -= 2.5;
40                    p = 2.81022636e-08;
41                    p = 3.43273939e-07 + p * w;
42                    p = -3.5233877e-06 + p * w;
43                    p = -4.39150654e-06 + p * w;
44                    p = 0.00021858087 + p * w;
45                    p = -0.00125372503 + p * w;
46                    p = -0.00417768164 + p * w;
47                    p = 0.246640727 + p * w;
48                    p = 1.50140941 + p * w;
49                } else {
50                    w = w.sqrt() - 3.0;
51                    p = -0.000200214257;
52                    p = 0.000100950558 + p * w;
53                    p = 0.00134934322 + p * w;
54                    p = -0.00367342844 + p * w;
55                    p = 0.00573950773 + p * w;
56                    p = -0.0076224613 + p * w;
57                    p = 0.00943887047 + p * w;
58                    p = 1.00167406 + p * w;
59                    p = 2.83297682 + p * w;
60                }
61
62                let res_ra = p * self;
63
64                let fx: Self = res_ra.error() - self;
65                let df = 2.0 / SQRT_PI as $kind * (-(res_ra * res_ra)).exp();
66                let d2f = -2.0 * res_ra * df;
67
68                res_ra - (2.0 * fx * df) / ((2.0 * df * df) - (fx * d2f))
69            }
70        }
71    };
72}
73
74implement!(f32);
75implement!(f64);
76
77#[cfg(test)]
78mod tests {
79    use assert;
80
81    use super::*;
82
83    #[test]
84    fn inv_error_negative() {
85        assert::close(-0.99.inv_error(), -1.8213863677184492, 1e-12);
86        assert::close(-0.999.inv_error(), -2.3267537655135242, 1e-12);
87    }
88
89    #[test]
90    fn inv_error_positive() {
91        assert::close(0.5.inv_error(), 0.47693627620446982, 1e-12);
92        assert::close(0.121.inv_error(), 0.10764782605515244, 1e-12);
93    }
94
95    #[test]
96    fn inv_error_zero() {
97        assert::close(0.0.inv_error(), 0.0, 1e-12);
98    }
99}
special/error.rs

special/
error.rs
