mathru/special/
error.rs

1//! Provides the [error](https://en.wikipedia.org/wiki/Error_function) and related functions
2use crate::algebra::abstr::Real;
3use crate::algebra::abstr::Sign;
4use crate::elementary::Trigonometry;
5use crate::special::gamma;
6use crate::special::gamma::Gamma;
7
8/// Provides [error, also called the Gauss error function](https://en.wikipedia.org/wiki/Error_function) related functions
9pub trait Error {
10    /// Error function
11    ///
12    /// ```math
13    /// \operatorname{erf}(z) = \frac{2}{\sqrt\pi}\int_0^z e^{-t^2}\,dt.
14    /// ```
15    ///
16    /// [For more information](https://en.wikipedia.org/wiki/Error_function)
17    ///
18    /// # Arguments
19    ///
20    /// self:
21    ///
22    fn erf(self) -> Self;
23
24    /// Complementary error function
25    ///
26    /// ```math
27    /// \operatorname{erfc}(z) = 1 - \operatorname{erf}(z),
28    /// ```
29    /// # Arguments
30    ///
31    /// self
32    fn erfc(self) -> Self;
33
34    /// Inverse error function
35    ///
36    /// ```math
37    /// \operatorname{erfinv}(x) = \operatorname{erf}^{-1}(x) \quad \text{for} \quad x \in (-1, 1)
38    /// ```
39    ///
40    /// [For more information](https://en.wikipedia.org/wiki/Error_function#Inverse_functions)
41    ///
42    /// # Arguments
43    ///
44    /// -1.0 < x < 1.0
45    ///
46    /// # Panics
47    ///
48    /// if x < -1.0 or x > 1.0
49    ///
50    fn erfinv(self) -> Self;
51
52    /// Inverse complementary error function
53    ///
54    /// ```math
55    /// \operatorname{erfc}^{-1}(1-z) = \operatorname{erf}^{-1}(z)
56    /// ```
57    /// # Arguments
58    ///
59    /// 0 <= x <= 2.0
60    ///
61    /// # Panics
62    ///
63    /// if x < 0.0 or x > 2.0
64    ///
65    fn erfcinv(self) -> Self;
66}
67
68impl Error for f64 {
69    fn erf(self) -> Self {
70        if self == 0.0 {
71            return self;
72        }
73        self.sign() * (1.0f64 - gamma::gamma_u(0.5f64, self * self) / Self::pi().sqrt())
74    }
75
76    fn erfc(self) -> Self {
77        1.0 - self.erf()
78    }
79
80    /// <https://en.wikipedia.org/wiki/Error_function#Inverse_functions>
81    /// Using the rational approximations tabulated in:
82    /// J. M. Blair, C. A. Edwards, and J. H. Johnson,
83    /// "Rational Chebyshev approximations for the inverse of the error function",
84    fn erfinv(self) -> Self {
85        let factors_leq075_p: [f64; 7] = [
86            1.603_049_558_440_662_3e1,
87            -9.078_495_926_296_033e1,
88            1.864_491_486_162_098_7e2,
89            -1.690_014_273_464_238e2,
90            6.545_466_284_794_487_e1,
91            -8.642_130_115_872_478,
92            0.176_058_782_139_059,
93        ];
94
95        let factors_leq075_q: [f64; 7] = [
96            1.478_064_707_151_383_1e1,
97            -9.137_416_702_426_032e1,
98            2.101_579_048_620_532e2,
99            -2.221_025_412_185_513_2e2,
100            1.076_045_391_605_512_4e2,
101            -2.060_107_303_282_654e1,
102            0.1e1,
103        ];
104
105        let factors_leg09375_p: [f64; 8] = [
106            -1.523_892_634_407_261_2e-2,
107            0.344_455_692_413_612_5,
108            -2.934_439_867_254_247_8,
109            1.176_350_570_521_782_8e1,
110            -2.265_529_282_310_110_3e1,
111            1.912_133_439_658_033e1,
112            -5.478_927_619_598_319,
113            0.23751_66890_24448,
114        ];
115
116        let factors_leg09375_q: [f64; 8] = [
117            -1.084_651_696_020_599_5e-2,
118            0.261_062_888_584_307_9,
119            -2.406_831_810_439_376,
120            1.069_512_997_338_701_5e1,
121            -2.371_671_552_159_658e1,
122            2.464_015_894_391_728_5e1,
123            -1.001_437_634_978_307e1,
124            0.1e1,
125        ];
126
127        let factors_p: [f64; 9] = [
128            1.050_131_152_373_343_8e-4,
129            1.053_261_131_423_333_8e-2,
130            0.269_878_027_362_432_8,
131            2.326_869_578_891_969,
132            7.167_854_794_910_8,
133            8.547_561_182_216_782,
134            6.873_808_807_354_384,
135            3.627_002_483_095_871,
136            0.886_062_739_296_515_5,
137        ];
138
139        let factors_q: [f64; 10] = [
140            1.050_126_668_703_033_8e-4,
141            1.053_286_230_093_332_7e-2,
142            0.270_198_623_737_515_6,
143            2.350_143_639_797_025_2,
144            7.607_802_878_580_127,
145            1.118_158_610_405_690_8e1,
146            1.194_878_791_843_539_6e1,
147            8.192_240_974_726_99,
148            4.099_387_907_636_801,
149            0.1e1,
150        ];
151
152        let a: f64 = self.abs();
153        if a >= 1.0 {
154            if self == 1.0 {
155                return f64::INFINITY;
156            } else if self == -1.0 {
157                return f64::NEG_INFINITY;
158            }
159
160            panic!("|self| has to be <= 1.0")
161        } else if a <= 0.75 {
162            let t: f64 = self * self - 0.5625;
163
164            self * horner(t, &factors_leq075_p) / horner(t, &factors_leq075_q)
165        } else if a <= 0.9375 {
166            let t: f64 = self * self - 0.87890625;
167
168            self * horner(t, &factors_leg09375_p) / horner(t, &factors_leg09375_q)
169        } else {
170            let t: f64 = 1.0 / (-(1.0 - a).ln()).sqrt();
171
172            horner(t, &factors_p) / (copysign(t, self) * horner(t, &factors_q))
173        }
174    }
175
176    fn erfcinv(self) -> Self {
177        if !(0.0f64..=2.0f64).contains(&self) {
178            panic!("self is > 2.0 or < 0.0");
179        }
180
181        (1.0 - self).erfinv()
182    }
183}
184
185impl Error for f32 {
186    fn erf(self) -> Self {
187        if self == 0.0 {
188            return self;
189        }
190        self.sign() * (1.0f32 - gamma::gamma_u(0.5f32, self * self) / Self::pi().sqrt())
191    }
192
193    fn erfc(self) -> Self {
194        1.0 - self.erf()
195    }
196
197    fn erfinv(self) -> Self {
198        let factors_leq075_p: [f32; 3] = [-1.309_599_7e1, 0.267_852_25e2, -9.289_057];
199
200        let factors_leq075_q: [f32; 4] = [-1.207_494_3e1, 3.096_061_5e1, -1.714_997_9e1, 0.1e1];
201
202        let factors_leg09375_p: [f32; 4] =
203            [-1.240_256_5e-1, 1.068_805_9, -1.959_455_6, 4.230_581_2e-1];
204
205        let factors_leg09375_q: [f32; 4] = [-8.827_698e-2, 8.900_743e-1, -2.175_703, 0.1e1];
206
207        let factors_p: [f32; 6] = [
208            1.550_47e-1,
209            1.382_719_6,
210            6.909_693_5e-1,
211            -1.128_081_4,
212            6.805_442_6e-1,
213            -1.644_415_7e-1,
214        ];
215
216        let factors_q: [f32; 3] = [1.550_248_6e-1, 1.385_228_2, 0.1e1];
217
218        let a: f32 = self.abs();
219        if a >= 1.0 {
220            if self == 1.0 {
221                return f32::INFINITY;
222            } else if self == -1.0 {
223                return f32::NEG_INFINITY;
224            }
225
226            panic!("|self| has to be <= 1.0")
227        } else if a <= 0.75 {
228            let t: f32 = self * self - 0.5625;
229
230            self * horner(t, &factors_leq075_p) / horner(t, &factors_leq075_q)
231        } else if a <= 0.9375 {
232            let t: f32 = self * self - 0.87890625;
233
234            self * horner(t, &factors_leg09375_p) / horner(t, &factors_leg09375_q)
235        } else {
236            let t: f32 = 1.0 / ((-((1.0 - a).ln())).sqrt());
237
238            horner(t, &factors_p) / (copysign(t, self) * horner(t, &factors_q))
239        }
240    }
241
242    fn erfcinv(self) -> Self {
243        if !(0.0f32..=2.0f32).contains(&self) {
244            panic!("self is > 2.0 or < 0.0");
245        }
246        (1.0 - self).erfinv()
247    }
248}
249
250/// Error Function
251///
252/// ```math
253/// \operatorname{erf}(z) = \frac{2}{\sqrt\pi}\int_0^z e^{-t^2}\,dt.
254/// ```
255///
256/// [For more information](https://en.wikipedia.org/wiki/Error_function)
257///
258/// # Example
259///
260/// ```
261/// use mathru::special::error;
262///
263/// let x: f64 = 0.0;
264/// let error: f64 = error::erf(x);
265/// ```
266pub fn erf<T>(x: T) -> T
267where
268    T: Real + Error + Gamma,
269{
270    x.erf()
271}
272
273/// Complementary error function
274///
275/// ```math
276/// \operatorname{erfc}(z) = 1 - \operatorname{erf} z,
277/// ```
278/// # Arguments
279///
280/// x
281///
282/// # Example
283///
284/// ```
285/// use mathru::special::error;
286///
287/// let x: f64 = 0.0;
288/// let error: f64 = error::erfc(x);
289/// ```
290pub fn erfc<T>(x: T) -> T
291where
292    T: Real + Error + Gamma,
293{
294    x.erfc()
295}
296
297/// Inverse error function
298///
299/// ```math
300/// \operatorname{erfinv}(x) = \operatorname{erf}^{-1}(x) \quad \text{for} \quad x \in (-1, 1)
301/// ```
302///
303/// # Arguments
304///
305/// -1.0 <= x <= 1.0
306///
307/// # Panics
308///
309/// if x < -1.0 or x > 1.0
310///
311/// # Example
312///
313/// ```
314/// use mathru::special::error;
315///
316/// let x: f64 = 0.0;
317/// let error: f64 = error::erfinv(x);
318/// ```
319pub fn erfinv<T>(x: T) -> T
320where
321    T: Real + Error,
322{
323    x.erfinv()
324}
325
326fn horner<T>(x: T, c: &[T]) -> T
327where
328    T: Real,
329{
330    c.iter()
331        .rev()
332        .skip(1)
333        .fold(*c.last().unwrap(), |s, c| s * x + *c)
334}
335
336fn copysign<T>(a: T, b: T) -> T
337where
338    T: Real,
339{
340    let a_neg: bool = a < T::zero();
341    let b_neg: bool = b < T::zero();
342    if (a_neg && !b_neg) || (!a_neg && b_neg) {
343        return -a;
344    }
345    a
346}
347
348/// Inverse complementary error function
349///
350/// ```math
351/// \operatorname{erfc}^{-1}(1-z) = \operatorname{erf}^{-1}(z)
352/// ```
353/// # Arguments
354///
355/// 0 <= x <= 2.0
356///
357/// # Panics
358///
359/// if x < 0.0 or x > 2.0
360///
361/// # Example
362///
363/// ```
364/// use mathru::special::error;
365///
366/// let x: f64 = 1.0;
367/// let error: f64 = error::erfcinv(x);
368/// ```
369pub fn erfcinv<T>(x: T) -> T
370where
371    T: Real + Error,
372{
373    x.erfcinv()
374}
mathru/special/error.rs

mathru/special/
error.rs
