1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
use crate::special::{Gamma, Probability};
use crate::Functions;
/// Provides methods for calculating polygamma functions.
pub struct Polygamma;
impl Polygamma {
/// Calculates the digamma function.
///
/// The digamma function, denoted as ψ(z), is the derivative of the natural logarithm of the gamma function.
///
/// # Parameters
///
/// - `z`: The value at which to calculate the digamma function.
///
/// # Returns
///
/// The value of the digamma function at the given parameter `z`.
///
/// # Example
///
/// ```rust
/// use numerilib::special::Polygamma;
///
/// let z = 2.0;
/// let digamma = Polygamma::digamma(z);
///
/// println!("Digamma({}) = {}", z, digamma);
/// ```
/// <hr/>
pub fn digamma(z: f64) -> f64 {
let p1 = |t: f64| Gamma::lanczosln(t).exp();
let dt = Functions::derivative(p1, z);
dt / Gamma::lanczosln(z).exp()
}
/// Calculates the polygamma function of a given degree.
///
/// The polygamma function, denoted as ψ^(n)(z), is the n-th derivative of the digamma function.
///
/// # Parameters
///
/// - `degree`: The degree of the polygamma function. Use -1 for the digamma function and 0 for the first derivative (digamma function).
/// - `z`: The value at which to calculate the polygamma function.
///
/// # Returns
///
/// The value of the polygamma function of the specified degree at the given parameter `z`.
///
/// # Example 1 (Digamma):
///
/// ```rust
/// use numerilib::special::Polygamma;
///
/// let degree = 0; // 0 corresponds to the digamma function
/// let z = 5.0;
///
/// let digamma = Polygamma::polygamma(degree, z);
///
/// println!("Digamma({}) = {}", z, digamma);
/// ```
///
/// # Example 2:
///
/// ```rust
/// use numerilib::special::Polygamma;
///
/// let degree = 1;
/// let z = 5.0;
///
/// let polygamma = Polygamma::polygamma(degree, z);
///
/// println!("Polygamma({},{}) = {}", degree, z, polygamma);
/// ```
/// <hr/>
pub fn polygamma(degree: i32, z: f64) -> f64 {
if degree == -1 {
Gamma::lanczosln(z)
} else if degree == 0 {
return Polygamma::digamma(z);
} else {
let func = |x: f64| {
(-1_f64).powi(degree + 1)
* Probability::factorial(degree as f64)
* (1_f64 / (z + x).powi(degree + 1))
};
Functions::summation(0_f64, 99_f64, func)
}
}
}