number_diff/functions/utils/
useful_functions.rs

1use std::f64::{consts::E, NAN};
2
3use crate::{Elementary::*, Integrate, EULER_MASCHERONI};
4
5/// An infinit number in the complex plane with an unknown or undefined complex argument.
6///
7/// For instance will approach complex infinity as x approaches 0. The limit could be either plus
8/// or minus infinity wich makes its magnitude infinit and its complex argument undefined.
9///
10/// See [this article](https://mathworld.wolfram.com/ComplexInfinity.html) for further information.
11pub const COMPLEX_INFINITY: f64 = NAN;
12
13/// returns n! for numbers n ∈ ℕ
14fn factorial_integer(numb: u128) -> u128 {
15    if numb == 0 {
16        1
17    } else {
18        let mut res: u128 = 1;
19        for i in 1..=numb {
20            res *= i;
21        }
22
23        res
24    }
25}
26
27// TODO: make the gamma function work for values < 1
28/// Returns the value of 𝜞(z) as defined by ∫t^(z-1)e^(-t)dt evaluated from 0 to ∞.
29pub fn gamma_function(z: f64) -> f64 {
30    let inner_funciton = Mul(
31        Pow(X.into(), Sub(Con(z).into(), Con(1.).into()).into()).into(),
32        Pow(Con(E).into(), Mul(X.into(), Con(-1.).into()).into()).into(),
33    );
34
35    // for whole numbers
36    if z.fract() == 0.0 {
37        // fraction part of the number is zero, meaning that the number is an integer
38        inner_funciton
39            .integrate()
40            .set_lower_bound(0.)
41            .set_upper_bound(100.)
42            .set_precision(1000)
43            .evaluate()
44            .unwrap()
45            .round_to(0)
46    } else {
47        inner_funciton
48            .integrate()
49            .set_lower_bound(0.)
50            .set_upper_bound(100.)
51            .set_precision(100000)
52            .evaluate()
53            .unwrap()
54    }
55}
56
57/// the polygamma function 𝛙m(z) describes the relationship between 𝜞(z) and its derivatives. For instance 𝛙0(z) =
58/// 𝜞'(z)/𝜞(z). [See article](https://en.wikipedia.org/wiki/Polygamma_function)
59pub fn polygamma_function(z: f64, m: usize) -> f64 {
60    if m == 0 {
61        digamma_function(z)
62    } else {
63        let inner_funciton = Pow(Log(Con(E).into(), X.into()).into(), Con(m as f64).into())
64            * Pow(X.into(), Con(z - 1.).into())
65            / (Con(1.) - X);
66
67        -inner_funciton
68            .integrate()
69            .set_lower_bound(1e-10)
70            .set_upper_bound(1. - 1e-10)
71            .set_precision(10)
72            .evaluate()
73            .unwrap()
74    }
75}
76
77/// Special case of the polygamma function 𝛙m(z) where m=0, The integral definition of the function
78/// then changes. [See article](https://en.wikipedia.org/wiki/Digamma_function#Integral_representations)
79pub fn digamma_function(z: f64) -> f64 {
80    let inner_funciton = (Con(1.) - Pow(X.into(), Con(z - 1.).into())) / (Con(1.) - X);
81
82    let integral_value = inner_funciton
83        .integrate()
84        .set_lower_bound(0.)
85        .set_upper_bound(1. - 1e-10)
86        .set_precision(1000)
87        .evaluate()
88        .unwrap();
89
90    integral_value - EULER_MASCHERONI
91}
92
93/// Allows the usage of factorials i.e. `x!`
94/// usually defined as:
95/// * x! = x*(x-1)!
96/// * 0! = 1
97///
98/// for x ∈ ℕ.
99///
100/// The definition of the factorial function can be expanded to the domain of all real numbers
101/// using the [gamma function](https://en.wikipedia.org/wiki/Gamma_function):
102///
103/// x! = 𝜞(x+1)
104///
105pub trait Factorial {
106    type Output;
107    fn factorial(&self) -> Self::Output;
108}
109
110/// implement factorial method for natural number types
111macro_rules! impl_factorial_natural {
112    (for $($t:ty), +) => {
113        $(impl Factorial for $t {
114            type Output = u128;
115            fn factorial(&self) -> Self::Output {
116                factorial_integer(self.clone() as u128)
117            }
118        })*
119    };
120}
121impl_factorial_natural!(for u8, u16, u32, u64, u128, usize);
122
123/// implement factorial for integer number types
124macro_rules! impl_factorial_integer {
125    (for $($t:ty), +) => {
126        $(impl Factorial for $t {
127            type Output = f64;
128            fn factorial(&self) -> Self::Output {
129                if self.is_negative() {
130                    // the continuos factorial funciton approaches ±∞ (complex infinity) for any
131                    // negative integer
132                    NAN
133                } else {
134                    factorial_integer(self.clone() as u128) as f64
135                }
136            }
137        })*
138    };
139}
140impl_factorial_integer!(for i8, i16, i32, i64, i128, isize);
141
142/// implement factorial method for float types
143macro_rules! impl_factorial_float {
144    (for $($t:ty), +) => {
145        $(impl Factorial for $t {
146            type Output = Self;
147            fn factorial(&self) -> Self::Output{
148                gamma_function(*self as f64 + 1.) as Self
149            }
150        })*
151    }
152}
153impl_factorial_float!(for f32, f64);
154
155/// Allows the usage of rounding methods that are more specific than rust std's round() method.
156pub trait Round {
157    /// Rounds self (a number) to the given number of decimal places. This method is
158    /// mainly made for the f32 and f64 types since integer types already have no decimal places.
159    ///
160    /// Example:
161    /// ```rust
162    /// assert_eq!(23.3274.round_to(2), 23.33);
163    ///
164    /// assert_eq!((1. / 3.).round_to(5), 0.33333);
165    ///
166    /// // For integer types, rounding to a decimal point is the same as casting it to f64
167    /// assert_eq!(100_u8.round_to(10), 100.);
168    /// ```
169    fn round_to(&mut self, decimal_places: u32) -> f64;
170
171    /// Rounds self (a number) to the given number of significant figures.
172    /// Example:
173    /// ```rust
174    /// assert_eq!(14912387964_u128.with_significant_figures(5), 14912000000);
175    ///
176    /// assert_eq!(-4095_i32.with_significant_figures(1), -4000);
177    ///
178    /// assert_eq!(1234.5678_f64.with_significant_figures(6), 1234.57)
179    /// ```
180    fn with_significant_figures(&mut self, digits: u64) -> Self;
181}
182
183macro_rules! impl_round_float {
184    (for $($t:ty), +) => {
185        $(impl Round for $t {
186            fn round_to(&mut self, decimal_places: u32) -> f64 {
187                let mut value = *self as f64;
188                // move the decimal point
189                value *= 10_usize.pow(decimal_places) as f64;
190                // round to an integer
191                value = value.round();
192                // move the decimal point back
193                value /= 10_usize.pow(decimal_places) as f64;
194                // give self the value of the rounded number
195                *self = value as Self;
196                value
197            }
198
199            fn with_significant_figures(&mut self, digits: u64) -> Self {
200                let value = if *self >= 0. {
201
202                    let order = (*self).log10().trunc() as i32;
203                        if digits as i32 <= order {
204                            ((*self) as isize).with_significant_figures(digits) as Self
205                        } else {
206                            if *self >= 1. {
207                                (*self * (10 as Self).powi((digits as i32 - order -1) as i32)).round() / (10 as Self).powi((digits as i32 - order -1) as i32)
208                            } else {
209                                (*self * (10 as Self).powi((digits as i32 - order) as i32)).round() / (10 as Self).powi((digits as i32 - order) as i32)
210                            }
211                        }
212                } else {
213                    -1. *(*self *-1.).with_significant_figures(digits)
214                };
215
216                *self = value as Self;
217                value
218            }
219        })*
220    };
221}
222impl_round_float!(for f32, f64);
223
224macro_rules! impl_round_int {
225    (for $($t:ty), +) => {
226        $(impl Round for $t {
227            #[allow(unused_variables)]
228            fn round_to(&mut self, decimal_places: u32) -> f64 {
229                *self as f64
230            }
231            fn with_significant_figures(&mut self, digits: u64) -> Self {
232                // move the decimal point to the appropriate spot so that we can round and then
233                // move it back
234                let order = (*self).ilog10() as u64;
235                let new_value = ((*self) as f64 / 10_f64.powi((order - digits +1) as i32)).round() * 10_f64.powi((order - digits +1) as i32);
236                // set new value
237                *self = new_value as Self;
238                *self
239            }
240        })*
241
242    };
243}
244impl_round_int!(for u8, u16, u32, u64, u128, usize, i8, i16, i32, i64, i128, isize);
number_diff/functions/utils/useful_functions.rs

number_diff/functions/utils/
useful_functions.rs
