adic 0.5.1

Arithmetic and rootfinding for p-adic numbers
Documentation 
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
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137

use crate::{
    error::{AdicError, AdicResult},
    normed::{Normed, UltraNormed},
    traits::{AdicPrimitive, CanTruncate},
    QAdic, ZAdic,
};
use super::iwasawa_log;


/// Adic digamma for x in Z_p
///
/// The [Digamma Function](<https://en.wikipedia.org/wiki/Digamma_function>) is the derivative of the
///  log of the Gamma function.
/// We calculate it using the [Volkenborn Integral](<https://en.wikipedia.org/wiki/Volkenborn_integral>)
///  of the [`iwasawa_log`](crate::function::special::iwasawa_log) function (Robert, p. 379).
///
/// The [Euler-Mascheroni constant](<https://en.wikipedia.org/wiki/Euler%27s_constant>) is equal to `-digamma(1)`.
/// This will not be the standard constant, it will be in Z_p, different for different p-adic numbers.
///
/// ```
/// # use adic::{function::special::adic_digamma, UAdic, ZAdic};
/// # fn main() -> Result<(), Box<dyn std::error::Error>> {
/// let digamma_five = adic_digamma(UAdic::new(5, vec![0, 1]), 4).unwrap();
/// let digamma_six = adic_digamma(UAdic::new(5, vec![1, 1]), 4).unwrap();
/// assert_eq!(digamma_six, digamma_five);
/// let euler_gamma_5 = ZAdic::new_approx(5, 4, vec![0, 1, 0, 3]);
/// assert_eq!(euler_gamma_5, -adic_digamma(UAdic::new(5, vec![1]), 4)?);
/// # Ok(()) }
/// ```
pub fn adic_digamma<Z>(x: Z, stopping_index : usize) -> AdicResult<ZAdic>
where Z: Into<ZAdic> {

    let x = QAdic::<ZAdic>::from_integer(x.into());
    let p = x.p();
    let zero = QAdic::zero(p);
    let one = QAdic::one(p);
    let mut sum = zero.clone();
    let mut t = x;

    //See if we can extend this to Q_p

    let integration_stop = usize::try_from(u32::from(p))?.pow(u32::try_from(stopping_index)?);
    for _ in 0..integration_stop {
        //We're not sure how quickly this converges but 2 * stopping_index seems like a good starting point
        if t.is_unit(){
            sum = sum + iwasawa_log(t.clone(), 2 * stopping_index)?;
        }
        t = t + one.clone();
    }
    let p_neg_n = QAdic::new(ZAdic::one(p), -isize::try_from(stopping_index)?);
    sum = sum * p_neg_n;
    if (sum.valuation() < 0.into()){
        Err(AdicError::IllDefined("Returns a fractional component, not in Z_p".to_string()))
    } else {
        Ok(sum.into_split(0).1)
    }
}



#[cfg(test)]
mod test {

    use num::Rational32;
    use crate::{
        traits::{AdicPrimitive, CanApproximate, TryPrimedFrom},
        RAdic, ZAdic,
    };
    use super::adic_digamma;

    #[test]
    fn euler_mascheroni() {

        let eg_constant_3 = zadic_approx!(3, 6, [0, 2, 2, 1, 2, 1]);
        let eg_constant_5 = zadic_approx!(5, 4, [0, 1, 0, 3]);
        let eg_constant_7 = zadic_approx!(7, 4, [5, 2, 4, 6]);
        let eg_constant_11 = zadic_approx!(11, 3, [1, 10, 2]);
        let eg_constant_13 = zadic_approx!(13, 3, [0, 4, 0]);
        assert_eq!(Ok(-eg_constant_3), adic_digamma(uadic!(3, [1]), 6));
        assert_eq!(Ok(-eg_constant_5), adic_digamma(uadic!(5, [1]), 4));
        assert_eq!(Ok(-eg_constant_7), adic_digamma(uadic!(7, [1]), 4));
        assert_eq!(Ok(-eg_constant_11), adic_digamma(uadic!(11, [1]), 3));
        assert_eq!(Ok(-eg_constant_13), adic_digamma(uadic!(13, [1]), 3));

    }

    #[test]
    fn simple() {
        let digamma_one = adic_digamma(uadic!(5, [1]), 4).unwrap();
        let digamma_two = adic_digamma(uadic!(5, [2]), 4).unwrap();
        let digamma_three = adic_digamma(uadic!(5, [3]), 4).unwrap();
        let digamma_four = adic_digamma(uadic!(5, [4]), 4).unwrap();
        let digamma_five = adic_digamma(uadic!(5, [0, 1]), 4).unwrap();
        let digamma_six = adic_digamma(uadic!(5, [1, 1]), 4).unwrap();
        let digamma_seven = adic_digamma(uadic!(5, [2, 1]), 4).unwrap();
        let digamma_eight = adic_digamma(uadic!(5, [3, 1]), 4).unwrap();
        let digamma_nine = adic_digamma(uadic!(5, [4, 1]), 4).unwrap();
        let digamma_ten = adic_digamma(uadic!(5, [0, 2]), 4).unwrap();
        assert_eq!(digamma_two, zadic_approx!(5, 4, [1, 4, 4, 1]));
        assert_eq!(digamma_two.clone() - digamma_one, zadic_approx!(5, 4, [1, 0, 0, 0]));
        assert_eq!(
            digamma_three.clone() - digamma_two,
            RAdic::try_primed_from(5, Rational32::new(1, 2)).unwrap().into_approximation(4)
        );
        assert_eq!(
            digamma_four.clone() - digamma_three,
            RAdic::try_primed_from(5, Rational32::new(1, 3)).unwrap().into_approximation(4)
        );
        assert_eq!(
            digamma_five.clone() - digamma_four,
            RAdic::try_primed_from(5, Rational32::new(1, 4)).unwrap().into_approximation(4)
        );
        assert_eq!(digamma_six, digamma_five);
        assert_eq!(
            digamma_seven.clone() - digamma_six,
            RAdic::try_primed_from(5, Rational32::new(1, 6)).unwrap().into_approximation(4)
        );
        assert_eq!(
            digamma_eight.clone() - digamma_seven,
            RAdic::try_primed_from(5, Rational32::new(1, 7)).unwrap().into_approximation(4)
        );
        assert_eq!(
            digamma_nine.clone() - digamma_eight,
            RAdic::try_primed_from(5, Rational32::new(1, 8)).unwrap().into_approximation(4)
        );
        assert_eq!(
            digamma_ten.clone() - digamma_nine,
            RAdic::try_primed_from(5, Rational32::new(1, 9)).unwrap().into_approximation(4)
        );
        let di_ten_minus_five = digamma_ten - digamma_five;
        let fraction_sums = (6..10).map(
            |n| RAdic::try_primed_from(5, Rational32::new(1, n)).unwrap().into_approximation(4)
        ).fold(ZAdic::zero(5), |acc , r| acc + r);
        assert_eq!(di_ten_minus_five, fraction_sums);
    }

}