arr_rs/math/operations/
special.rs

1use crate::{
2    core::prelude::*,
3    errors::prelude::*,
4    numeric::prelude::*,
5};
6
7/// `ArrayTrait` - Array Math Special functions
8pub trait ArrayMathSpecial<N: NumericOps> where Self: Sized + Clone {
9
10    /// Modified Bessel function of the first kind, order 0
11    ///
12    /// # Examples
13    ///
14    /// ```
15    /// use arr_rs::prelude::*;
16    ///
17    /// let arr = Array::flat(vec![-2., 0., 3.5]);
18    /// assert_eq!(Array::flat(vec![2.2795851066228696, 0.9999999999999997, 7.378203432225479]), arr.i0());
19    /// ```
20    ///
21    /// # Errors
22    ///
23    /// may returns `ArrayError`
24    fn i0(&self) -> Result<Array<N>, ArrayError>;
25
26    /// Return the normalized sinc function
27    ///
28    /// # Examples
29    ///
30    /// ```
31    /// use arr_rs::prelude::*;
32    ///
33    /// let arr = Array::flat(vec![-1., 0., 1.]);
34    /// assert_eq!(Array::flat(vec![3.898_171_832_519_375_5e-17, 1.0, 3.898_171_832_519_375_5e-17]), arr.sinc());
35    /// ```
36    ///
37    /// # Errors
38    ///
39    /// may returns `ArrayError`
40    fn sinc(&self) -> Result<Array<N>, ArrayError>;
41}
42
43impl <N: NumericOps> ArrayMathSpecial<N> for Array<N> {
44
45    fn i0(&self) -> Result<Self, ArrayError> {
46        self.map(|&x| N::from(i0(x.to_f64())))
47    }
48
49    fn sinc(&self) -> Result<Self, ArrayError> {
50        self.map(|&x| {
51            let y = std::f64::consts::PI *
52                if x == N::zero() { 1.0e-20 }
53                else { x.to_f64() };
54            N::from(y.sin() / y)
55        })
56    }
57}
58
59impl <N: NumericOps> ArrayMathSpecial<N> for Result<Array<N>, ArrayError> {
60
61    fn i0(&self) -> Self {
62        self.clone()?.i0()
63    }
64
65    fn sinc(&self) -> Self {
66        self.clone()?.sinc()
67    }
68}
69
70const VALUES_I0_A: &[f64] = &[
71    -4.415_341_646_479_339E-18,
72    3.330_794_518_822_238E-17,
73    -2.431_279_846_547_954E-16,
74    1.715_391_285_555_133E-15,
75    -1.168_533_287_799_345E-14,
76    7.676_185_498_604_935E-14,
77    -4.856_446_783_111_929E-13,
78    2.955_052_663_129_639E-12,
79    -1.726_826_291_441_555E-11,
80    9.675_809_035_373_236E-11,
81    -5.189_795_601_635_262E-10,
82    2.659_823_724_682_386E-9,
83    -1.300_025_009_986_248E-8,
84    6.046_995_022_541_918E-8,
85    -2.670_793_853_940_611E-7,
86    1.117_387_539_120_103E-6,
87    -4.416_738_358_458_75E-6,
88    1.644_844_807_072_889E-5,
89    -5.754_195_010_082_103E-5,
90    1.885_028_850_958_416E-4,
91    -5.763_755_745_385_823E-4,
92    1.639_475_616_941_335E-3,
93    -4.324_309_995_050_575E-3,
94    1.054_646_039_459_499E-2,
95    -2.373_741_480_589_946E-2,
96    4.930_528_423_967_07E-2,
97    -9.490_109_704_804_764E-2,
98    1.716_209_015_222_087E-1,
99    -3.046_826_723_431_983E-1,
100    6.767_952_744_094_76E-1,
101];
102
103const VALUES_I0_B: &[f64] = &[
104    -7.233_180_487_874_753E-18,
105    -4.830_504_485_944_182E-18,
106    4.465_621_420_296_759E-17,
107    3.461_222_867_697_461E-17,
108    -2.827_623_980_516_583E-16,
109    -3.425_485_619_677_219E-16,
110    1.772_560_133_056_526E-15,
111    3.811_680_669_352_622E-15,
112    -9.554_846_698_828_307E-15,
113    -4.150_569_347_287_222E-14,
114    1.540_086_217_521_409E-14,
115    3.852_778_382_742_142E-13,
116    7.180_124_451_383_666E-13,
117    -1.794_178_531_506_806E-12,
118    -1.321_581_184_044_771E-11,
119    -3.149_916_527_963_241E-11,
120    1.188_914_710_784_643E-11,
121    4.940_602_388_224_969E-10,
122    3.396_232_025_708_386E-9,
123    2.266_668_990_498_178E-8,
124    2.048_918_589_469_063E-7,
125    2.891_370_520_834_756E-6,
126    6.889_758_346_916_823E-5,
127    3.369_116_478_255_694E-3,
128    8.044_904_110_141_088E-1,
129];
130
131fn chbevl(x: f64, vals: &[f64]) -> f64 {
132    let mut b0 = vals[0];
133    let mut b1 = 0.;
134
135    let mut b2 = 0.;
136    for val in vals.iter().skip(1) {
137        b2 = b1;
138        b1 = b0;
139        b0 = x.mul_add(b1, -b2) + val;
140    }
141
142    0.5 * (b0 - b2)
143}
144
145fn i0(x: f64) -> f64 {
146    if x <= 8.0 { f64::exp(x) * chbevl(x / 2. - 2., VALUES_I0_A) }
147    else { f64::exp(x) * chbevl(32. / x - 2., VALUES_I0_B) / f64::sqrt(x) }
148}
arr_rs/math/operations/special.rs

arr_rs/math/operations/
special.rs
