ekore 0.15.3

EKO expressions
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
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162

//! Cache harmonic sums for given Mellin N.

use num::{complex::Complex, Zero};
use std::collections::HashMap;

use crate::harmonics::{g_functions, w1, w2, w3, w4, w5};

/// List of available elements.
#[derive(Debug, PartialEq, Eq, Hash)]
pub enum K {
    /// $S_1(N)$
    S1,
    /// $S_2(N)$
    S2,
    /// $S_3(N)$
    S3,
    /// $S_4(N)$
    S4,
    /// $S_5(N)$
    S5,
    /// $S_1(N/2)$
    S1h,
    /// $S_2(N/2)$
    S2h,
    /// $S_3(N/2)$
    S3h,
    /// $S_1((N-1)/2)$
    S1mh,
    /// $S_2((N-1)/2)$
    S2mh,
    /// $S_3((N-1)/2)$
    S3mh,
    /// $g_3(N)$
    G3,
    /// $S_{-1}(N)$ even moments
    Sm1e,
    /// $S_{-1}(N)$ odd moments
    Sm1o,
    /// $S_{-2}(N)$ even moments
    Sm2e,
    /// $S_{-2}(N)$ odd moments
    Sm2o,
    /// $S_{-3}(N)$ even moments
    Sm3e,
    /// $S_{-3}(N)$ odd moments
    Sm3o,
    /// $S_{-2,1}(N)$ even moments
    Sm21e,
    /// $S_{-2,1}(N)$ odd moments
    Sm21o,
}

/// Hold all cached values.
pub struct Cache {
    /// Mellin N
    n: Complex<f64>,
    /// Mapping
    m: HashMap<K, Complex<f64>>,
}

impl Cache {
    /// Initialize new, empty Cache at given Mellin N.
    pub fn new(n: Complex<f64>) -> Self {
        Self {
            n,
            m: HashMap::new(),
        }
    }

    /// Get Mellin N.
    pub fn n(&self) -> Complex<f64> {
        Complex::new(self.n.re, self.n.im)
    }

    /// Retrieve an element.
    pub fn get(&mut self, k: K) -> Complex<f64> {
        let val = self.m.get(&k);
        // already there?
        if let Some(value) = val {
            return *value;
        }
        // compute new
        let val = match k {
            K::S1 => w1::S1(self.n),
            K::S2 => w2::S2(self.n),
            K::S3 => w3::S3(self.n),
            K::S4 => w4::S4(self.n),
            K::S5 => w5::S5(self.n),
            K::S1h => w1::S1(self.n / 2.),
            K::S2h => w2::S2(self.n / 2.),
            K::S3h => w3::S3(self.n / 2.),
            K::S1mh => w1::S1((self.n - 1.) / 2.),
            K::S2mh => w2::S2((self.n - 1.) / 2.),
            K::S3mh => w3::S3((self.n - 1.) / 2.),
            K::G3 => g_functions::g3(self),
            K::Sm1e => w1::Sm1e(self),
            K::Sm1o => w1::Sm1o(self),
            K::Sm2e => w2::Sm2e(self),
            K::Sm2o => w2::Sm2o(self),
            K::Sm3e => w3::Sm3e(self),
            K::Sm3o => w3::Sm3o(self),
            K::Sm21e => w3::Sm21e(self),
            K::Sm21o => w3::Sm21o(self),
        };
        // insert
        self.m.insert(k, val);
        val
    }
}

/// Recursive computation of harmonic sums.
///
/// Compute the harmonic sum $S_{w}(N+k)$ stating from the value $S_{w}(N)$ via the recurrence relations.
pub fn recursive_harmonic_sum(
    base_value: Complex<f64>,
    n: Complex<f64>,
    iterations: usize,
    weight: u32,
) -> Complex<f64> {
    let mut fact = Complex::zero();
    for i in 1..iterations + 1 {
        fact += (1.0 / (n + (i as f64))).powu(weight);
    }
    base_value + fact
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::harmonics::{w1, w2, w3, w4};
    use crate::{assert_approx_eq_cmplx, cmplx};
    use num::complex::Complex;

    #[test]
    fn n() {
        let n = cmplx!(1., 0.);
        let c = Cache::new(n);
        let mut m = c.n();
        m += cmplx!(1., 0.);
        assert_approx_eq_cmplx!(f64, c.n(), n);
        assert_approx_eq_cmplx!(f64, m, cmplx!(2., 0.));
    }

    #[test]
    fn test_recursive_harmonic_sum() {
        const SX: [fn(Complex<f64>) -> Complex<f64>; 5] = [w1::S1, w2::S2, w3::S3, w4::S4, w5::S5];
        const NS: [Complex<f64>; 2] = [cmplx!(1.0, 0.0), cmplx!(2.34, 3.45)];
        const ITERS: [usize; 2] = [1, 2];
        for sit in SX.iter().enumerate() {
            for nit in NS.iter().enumerate() {
                let n = *nit.1;
                for iit in ITERS.iter().enumerate() {
                    let iterations = *iit.1;
                    let s_base = sit.1(n);
                    let s_test = sit.1(n + (iterations as f64));
                    let s_ref = recursive_harmonic_sum(s_base, n, iterations, 1 + (sit.0 as u32));
                    assert_approx_eq_cmplx!(f64, s_test, s_ref);
                }
            }
        }
    }
}