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

use num::complex::Complex;
use crate::cln::CLn;

/// Provides the 6-th order polylogarithm function `li6()` of a
/// number of type `T`.
pub trait Li6<T> {
    fn li6(&self) -> T;
}

impl Li6<Complex<f64>> for Complex<f64> {
    /// Returns the sixths order polylogarithm of a complex number of type
    /// `Complex<f64>`.
    ///
    /// # Example:
    /// ```
    /// use num::complex::Complex;
    /// use polylog::Li6;
    ///
    /// let z = Complex::new(1.0, 1.0);
    /// println!("Li6({}) = {}", z, z.li6());
    /// ```
    fn li6(&self) -> Complex<f64> {
        let pi  = std::f64::consts::PI;
        let pi2 = pi*pi;
        let z6  = 1.0173430619844491; // zeta(6)
        let bf  = [
            1.                    , -31./64.               ,
            1.5241340877914952e-01, -3.4365555877057613e-02,
            5.7174797239368999e-03, -6.8180453746570645e-04,
            4.9960361948734493e-05, -4.9166051196039048e-07,
           -3.0632975161302164e-07,  1.4414599270849095e-08,
            3.7272438230924107e-09, -3.7300867345487607e-10,
           -5.1246526816085832e-11,  9.0541930956636683e-12,
            6.7381882615512517e-13, -2.1215831150303135e-13,
           -6.8408811719011698e-15,  4.8691178462005581e-15
        ];

        if self.im == 0.0 {
            if self.re == 0.0 {
                return Complex::new(0., 0.);
            }
            if self.re == 1.0 {
                return Complex::new(z6, 0.);
            }
            if self.re == -1.0 {
                return Complex::new(-31./32.*z6, 0.0);
            }
        }

        let nz  = self.norm();
        let pz  = self.arg();
        let lnz = nz.ln();

        if lnz*lnz + pz*pz < 1.0 { // |log(z)| < 1
            let u  = Complex::new(lnz, pz);
            let u2 = u*u;
            let c1 = 1.0369277551433699; // zeta(5)
            let c2 = 0.54116161685556910;
            let c3 = 0.20034281719326571;
            let c4 = 0.068538919452009435;
            let c5 = (137./60. - (-u).cln())/120.;
            let c6 = -1./1440.;

            let cs = [
                -1.6534391534391534e-05, 2.2964432686654909e-08,
                -9.9413128513657614e-11, 6.6912682653423394e-13,
                -5.7933058574392549e-15
            ];

            return z6 + u * c1 +
                u2 * (c2 + u * c3 +
                u2 * (c4 + u * c5 +
                u2 * (c6 +
                u * (cs[0] +
                u2 * (cs[1] +
                u2 * (cs[2] +
                u2 * (cs[3] +
                u2 * (cs[4]))))))));
        }

        let (u, rest, sgn) = if nz <= 1.0 {
            (-(1.0 - self).cln(), Complex::new(0.0, 0.0), 1.)
        } else { // nz > 1.0
            let pi4 = pi2*pi2;
            let pi6 = pi2*pi4;
            let arg = if pz > 0.0 { pz - pi } else { pz + pi };
            let lmz = Complex::new(lnz, arg); // (-self).cln()
            let lmz2 = lmz*lmz;
            (-(1. - 1./self).cln(), -31.*pi6/15120. + lmz2*(-7./720.*pi4 + lmz2*(-1./144.*pi2 - 1./720.*lmz2)), -1.)
        };

        let u2 = u*u;
        let u4 = u2*u2;
        let u8 = u4*u4;

        rest + sgn * (
           u*bf[0] +
           u2*(bf[1] + u*bf[2]) +
           u4*(bf[3] + u*bf[4] + u2*(bf[5] + u*bf[6])) +
           u8*(bf[7] + u*bf[8] + u2*(bf[9] + u*bf[10]) +
               u4*(bf[11] + u*bf[12] + u2*(bf[13] + u*bf[14]))) +
           u8*u8*(bf[15] + u*bf[16] + u2*bf[17])
        )
    }
}