vml/special/func/func.rs
1use crate::special::Gamma;
2
3pub struct Functions;
4
5impl Functions {
6 /// A simple implementation of the Error Function that is used to Calculate The NormalCDF Function. <br>
7 /// Learn more about the Error Function at: <https://wikipedia.org/wiki/Error_function>
8 /// <hr/>
9 ///
10 /// # Example:
11 /// ```
12 /// use vml::special::Functions;
13 ///
14 /// let bound = 4_f64;
15 /// let low_erf = Functions::lower_erf(bound);
16 ///
17 /// assert_eq!(low_erf, -0.9999999845827403);
18 /// ```
19 /// <hr/>
20 ///
21 pub fn lower_erf(bound: f64) -> f64 {
22 let f = |x: f64| (std::f64::consts::E).powf(-1_f64 * (x.powi(2)));
23
24 let error_function: f64 = (2.0 / std::f64::consts::PI.sqrt()) * Functions::integral(bound, 0_f64, f).unwrap();
25 error_function
26 }
27 /// A traditional implementation of the standard Error Function. <br>
28 /// Learn more about the Error Function at: <https://wikipedia.org/wiki/Error_function>
29 /// <hr/>
30 ///
31 /// # Example:
32 /// ```
33 /// use vml::special::Functions;
34 ///
35 /// let bound = 4_f64;
36 /// let low_erf = Functions::erf(bound);
37 ///
38 /// assert_eq!(low_erf, 0.9999999845827289);
39 /// ```
40 /// <hr/>
41 ///
42 pub fn erf(bound: f64) -> f64 {
43 let f = |x: f64| (std::f64::consts::E).powf(-1_f64 * (x.powi(2)));
44
45 let error_function: f64 = (2.0 / std::f64::consts::PI.sqrt()) * Functions::integral(0_f64, bound, f).unwrap();
46 error_function
47 }
48 /// A Rust implementation of the this approximation by Alijah Ahmed at: <https://scistatcalc.blogspot.com/2013/09/numerical-estimate-of-inverse-error.html>
49 ///
50 /// # Example
51 ///
52 /// ```
53 /// use vml::special::Functions;
54 ///
55 /// let x = 0.975;
56 /// let inverf = Functions::inverf(x);
57 ///
58 /// assert_eq!(inverf, 1.584911068059619);
59 /// ```
60 /// <hr/>
61 ///
62 pub fn inverf(x: f64) -> f64 {
63 let w = -((1.0 - x) * (1.0 + x)).ln();
64 let mut p;
65
66 if w < 5.0 {
67 let w_minus_2_5 = w - 2.5;
68 p = 2.81022636e-08;
69 p = 3.43273939e-07 + p * w_minus_2_5;
70 p = -3.5233877e-06 + p * w_minus_2_5;
71 p = -4.39150654e-06 + p * w_minus_2_5;
72 p = 0.00021858087 + p * w_minus_2_5;
73 p = -0.00125372503 + p * w_minus_2_5;
74 p = -0.00417768164 + p * w_minus_2_5;
75 p = 0.246640727 + p * w_minus_2_5;
76 p = 1.50140941 + p * w_minus_2_5;
77 } else {
78 let sqrt_w_minus_3 = (w).sqrt() - 3.0;
79 p = -0.000200214257;
80 p = 0.000100950558 + p * sqrt_w_minus_3;
81 p = 0.00134934322 + p * sqrt_w_minus_3;
82 p = -0.00367342844 + p * sqrt_w_minus_3;
83 p = 0.00573950773 + p * sqrt_w_minus_3;
84 p = -0.0076224613 + p * sqrt_w_minus_3;
85 p = 0.00943887047 + p * sqrt_w_minus_3;
86 p = 1.00167406 + p * sqrt_w_minus_3;
87 p = 2.83297682 + p * sqrt_w_minus_3;
88 }
89
90 let res_ra = p * x;
91
92 let fx = Functions::erf(res_ra) - x;
93 let df = (2.0 / (std::f64::consts::PI).sqrt()) * (-res_ra * res_ra).exp();
94 let d2f = -2.0 * res_ra * df;
95
96 res_ra - (2.0 * fx * df) / ((2.0 * df * df) - (fx * d2f))
97 }
98 /// Uses the definition of a derivative to calculate the derivative of a function at a specific point of a given function. <br>
99 /// Learn more about Derivatives and Differentiation at: <https://wikipedia.org/wiki/Derivative>
100 /// <hr/>
101 ///
102 /// # Example:
103 /// ```
104 /// use vml::special::Functions;
105 ///
106 /// let function = |x: f64| x.powi(2);
107 /// let x = 2;
108 ///
109 /// let derivative = Functions::derivative(function, x).unwrap();
110 ///
111 /// assert_eq!(derivative, 4);
112 /// ```
113 /// <hr/>
114 ///
115 pub fn derivative(f: fn(f64) -> f64, x: impl Into<f64> + Copy) -> Result<i64, f64> {
116 let h = 1e-7;
117 let definition = ((f(x.into() + h)) - f(x.into())) / h;
118
119 if definition.fract().abs() < 1e-7 {
120 Ok(definition.round() as i64)
121 } else {
122 Err(definition)
123 }
124 }
125 /// Uses Simpson's 1/3rd Rule to calculate the integral.<br>
126 /// Learn more about Integrals at: <https://wikipedia.org/wiki/Integral> <br>
127 /// Learn more about Simpson's 1/3rd Rule at: <https://wikipedia.org/wiki/Simpson's_rule> <br>
128 /// <hr/>
129 ///
130 /// # Example
131 /// ```
132 /// use vml::special::Functions;
133 ///
134 /// let lower_bound = 0_f64;
135 /// let upper_bound = 6_f64;
136 /// let function = |x: f64| x.powi(2);
137 ///
138 /// let integral = Functions::integral(lower_bound, upper_bound, function).unwrap();
139 ///
140 /// assert_eq!(integral, 72_f64)
141 /// ```
142 /// <hr/>
143 ///
144 pub fn integral<F: Fn(f64) -> f64>(lower_bound: f64, upper_bound: f64, f: F) -> Result<f64, String> {
145 let lower = lower_bound;
146 let upper = upper_bound;
147
148 let n = 100000;
149 let h = (upper - lower) / n as f64;
150
151 let mut sum = f(lower) + f(upper);
152
153 for i in 1..n {
154 let x = lower + i as f64 * h;
155 sum += if i % 2 == 0 {
156 2.0 * f(x)
157 } else {
158 4.0 * f(x)
159 };
160 }
161
162 let r1 = (h / 3.0) * sum;
163
164 if r1.is_infinite() || r1.is_nan() {
165 return Err("Function is Divergent".to_string());
166 }
167
168 let r2 = (r1 * 1e+6).round();
169
170 if ((r1 * 1e+6).floor() + 1e-6) > r2 || ((r1 * 1e+6).ceil() - 1e-6) < r2 {
171 Ok(r2 / 1e+6)
172 } else {
173 Ok(r1)
174 }
175 }
176 /// Summations in Rust. <br>
177 /// Learn more at: <https://wikipedia.org/wiki/Summation>
178 /// <hr/>
179 ///
180 /// # Example #1: Constant
181 ///
182 /// ```
183 /// use vml::special::Functions;
184 ///
185 /// let start = 0;
186 /// let limit = 9;
187 /// let function = |x: f64| 3_f64;
188 ///
189 /// let summation = Functions::summation(start, limit, function);
190 ///
191 /// assert_eq!(summation, 30_f64);
192 /// ```
193 /// <hr/>
194 ///
195 /// # Example #2: Function
196 /// ```
197 /// use vml::special::Functions;
198 ///
199 /// let start = 4.5;
200 /// let limit = 100;
201 /// let function = |x: f64| 1_f64 / x;
202 ///
203 /// let summation = Functions::summation(start, limit, function);
204 ///
205 /// assert_eq!(summation, 3.1040441843062854);
206 /// ```
207 /// <hr/>
208 ///
209 pub fn summation(start: impl Into<f64> + Copy,limit: impl Into<f64> + Copy, f: fn(f64) -> f64) -> f64 {
210 let start = start.into().round() as i64;
211 let limit = limit.into().round() as i64;
212 let mut result = 0.0;
213
214 for i in start..=limit {
215 result += f(i as f64);
216 }
217
218 result
219 }
220 /// Calculates a Factorial by using Lanczos's Gamma Function Approximation. <br>
221 /// Learn more at: <https://wikipedia.org/wiki/Factorial>
222 /// <hr/>
223 ///
224 /// # Example:
225 /// ```
226 /// use vml::special::Functions;
227 ///
228 /// let n = 6_f64;
229 /// let factorial = Functions::factorial(n);
230 /// assert_eq!(factorial, 720_f64);
231 /// ```
232 /// <hr/>
233 ///
234 pub fn factorial(n: f64) -> f64 {
235 Gamma::lanczos(n + 1_f64)
236 }
237 /// Calculates the product of a function. a.k.a capital Pi Notation. <br>
238 /// Learn more at: <https://wikipedia.org/wiki/Product_(mathematics)#Product_of_a_sequence>
239 /// <hr/>
240 ///
241 /// # Example #1: Constant
242 ///
243 /// ```
244 /// use vml::special::Functions;
245 ///
246 /// let start = 2_f64;
247 /// let limit = 7_f64;
248 /// let f = |x: f64| 3_f64;
249 ///
250 /// let product_series = Functions::product(start, limit, f);
251 /// assert_eq!(product_series, 729_f64);
252 /// ```
253 /// <hr/>
254 ///
255 /// # Example #2: Function
256 ///
257 /// ```
258 /// use vml::special::Functions;
259 ///
260 /// let start = 1_f64;
261 /// let limit = 6_f64;
262 /// let f = |x: f64| x.powi(2);
263 ///
264 /// let product_series = Functions::product(start, limit, f);
265 /// assert_eq!(product_series, 518400_f64);
266 /// ```
267 /// <hr/>
268 ///
269 pub fn product<F: Fn(f64) -> f64>(start: f64, limit: f64, f: F) -> f64 {
270 let mut result = 1.0;
271 let start = start.round() as i64;
272 let limit = limit.round() as i64;
273
274 for i in start..=limit {
275 result *= f(i as f64);
276 }
277 result
278 }
279}