mathhook_core/functions/special/
digamma.rs

1use crate::{Expression, MathConstant, Number};
2
3/// Digamma function ψ(z) = Γ'(z)/Γ(z)
4///
5/// The digamma function is the logarithmic derivative of the Gamma function:
6/// ψ(z) = d/dz ln(Γ(z))
7///
8/// # Mathematical Properties
9///
10/// - ψ(1) = -γ (Euler-Mascheroni constant ≈ -0.5772156649)
11/// - ψ(n+1) = ψ(n) + 1/n for positive integers n
12/// - ψ(z+1) = ψ(z) + 1/z (recurrence relation)
13/// - ψ(1/2) = -γ - ln(4)
14///
15/// # Implementation
16///
17/// Uses reflection formula for z < 0.5 and series expansion for z ≥ 0.5.
18/// Special values are computed exactly when possible.
19///
20/// # Arguments
21///
22/// * `z` - Expression to evaluate digamma function at
23///
24/// # Examples
25///
26/// ```rust
27/// use mathhook_core::{Expression, Number};
28/// use mathhook_core::functions::special::digamma;
29///
30/// let z = Expression::Number(Number::Integer(1));
31/// let result = digamma(&z);
32/// ```
33pub fn digamma(z: &Expression) -> Expression {
34    match z {
35        Expression::Number(Number::Integer(1)) => Expression::mul(vec![
36            Expression::integer(-1),
37            Expression::constant(MathConstant::EulerGamma),
38        ]),
39        Expression::Number(Number::Integer(n)) if *n > 1 => {
40            let mut sum = Expression::mul(vec![
41                Expression::integer(-1),
42                Expression::constant(MathConstant::EulerGamma),
43            ]);
44            for k in 1..*n {
45                sum = Expression::add(vec![sum, Expression::rational(1, k)]);
46            }
47            sum
48        }
49        Expression::Number(Number::Float(x)) => {
50            let result = digamma_numerical(*x);
51            Expression::Number(Number::Float(result))
52        }
53        _ => Expression::function("digamma", vec![z.clone()]),
54    }
55}
56
57/// Numerically evaluates the digamma function using series expansion
58///
59/// Uses reflection formula for z < 0.5 and asymptotic series for z ≥ 0.5.
60///
61/// # Arguments
62///
63/// * `z` - Value to evaluate digamma at
64///
65/// # Returns
66///
67/// ψ(z) value
68pub fn digamma_numerical(mut z: f64) -> f64 {
69    if z.is_nan() || z.is_infinite() {
70        return f64::NAN;
71    }
72
73    if z <= 0.0 && (z - z.round()).abs() < 1e-10 {
74        return f64::NAN;
75    }
76
77    let mut result = 0.0;
78    while z < 10.0 {
79        result -= 1.0 / z;
80        z += 1.0;
81    }
82
83    let z_inv = 1.0 / z;
84    let z2_inv = z_inv * z_inv;
85
86    result += f64::ln(z) - 0.5 * z_inv;
87    result -= z2_inv
88        * (1.0 / 12.0 - z2_inv * (1.0 / 120.0 - z2_inv * (1.0 / 252.0 - z2_inv * (1.0 / 240.0))));
89
90    result
91}
92
93#[cfg(test)]
94mod tests {
95    use super::*;
96
97    #[test]
98    fn test_digamma_special_values() {
99        let result = digamma(&Expression::Number(Number::Integer(1)));
100        match result {
101            Expression::Mul(ref factors) if factors.len() == 2 => {
102                assert!(matches!(
103                    factors[0],
104                    Expression::Number(Number::Integer(-1))
105                ));
106                assert!(matches!(
107                    factors[1],
108                    Expression::Constant(MathConstant::EulerGamma)
109                ));
110            }
111            _ => panic!("ψ(1) should be -γ"),
112        }
113    }
114
115    #[test]
116    fn test_digamma_numerical() {
117        let result = digamma(&Expression::Number(Number::Float(1.0)));
118        match result {
119            Expression::Number(Number::Float(val)) => {
120                const EULER_GAMMA: f64 = 0.577_215_664_901_532_9;
121                assert!((val + EULER_GAMMA).abs() < 1e-10, "ψ(1) ≈ -γ, got {}", val);
122            }
123            _ => panic!("ψ(1.0) should return numerical float"),
124        }
125    }
126
127    #[test]
128    fn test_digamma_recurrence() {
129        const EULER_GAMMA: f64 = 0.577_215_664_901_532_9;
130        let psi_2 = digamma_numerical(2.0);
131        let expected = -EULER_GAMMA + 1.0;
132        assert!(
133            (psi_2 - expected).abs() < 1e-10,
134            "ψ(2) = ψ(1) + 1, got {}",
135            psi_2
136        );
137    }
138
139    #[test]
140    fn test_digamma_integer_formula() {
141        let result = digamma(&Expression::Number(Number::Integer(3)));
142        let simplified = result.to_string();
143        assert!(simplified.contains("EulerGamma"));
144    }
145}
mathhook_core/functions/special/digamma.rs

mathhook_core/functions/special/
digamma.rs
