polylog/
li5.rs

1use num::complex::Complex;
2use crate::cln::CLn;
3
4/// Provides the 5-th order polylogarithm function `li5()` of a
5/// number of type `T`.
6pub trait Li5<T> {
7    fn li5(&self) -> T;
8}
9
10impl Li5<Complex<f64>> for Complex<f64> {
11    /// Returns the fifth order polylogarithm of a complex number of type
12    /// `Complex<f64>`.
13    ///
14    /// # Example:
15    /// ```
16    /// use num::complex::Complex;
17    /// use polylog::Li5;
18    ///
19    /// assert!((Complex::new(1.0_f64, 1.0_f64).li5() - Complex::new(0.9874666591701124_f64, 1.0684416071074221_f64)).norm() < 2.0_f64*std::f64::EPSILON);
20    /// ```
21    fn li5(&self) -> Complex<f64> {
22        let pi  = std::f64::consts::PI;
23        let pi2 = pi*pi;
24        let z5  = 1.0369277551433699; // zeta(5)
25
26        if self.im == 0.0 && self.re == 0.0 {
27            *self
28        } else if self.im == 0.0 && self.re == 1.0 {
29            Complex::new(z5, self.im)
30        } else if self.im == 0.0 && self.re == -1.0 {
31            Complex::new(-15.0/16.0*z5, self.im)
32        } else {
33            let nz  = self.norm();
34            let pz  = self.arg();
35            let lnz = nz.ln();
36
37            if lnz*lnz + pz*pz < 1.0 { // |log(z)| < 1
38                let u  = Complex::new(lnz, pz);
39                let u2 = u*u;
40                let c1 = 1.0823232337111382; // zeta(4)
41                let c2 = 0.60102845157979714; // zeta(3)/2
42                let c3 = 0.27415567780803774;
43                let c4 = (25.0/12.0 - (-u).cln())/24.0;
44                let c5 = -1.0/240.0;
45                let cs = [
46                    -1.1574074074074074e-04, 2.0667989417989418e-07,
47                    -1.0935444136502338e-09, 8.6986487449450412e-12,
48                    -8.6899587861588824e-14, 1.0081254080218813e-15
49                ];
50
51                z5 + u * c1 +
52                u2 * (c2 + u * c3 +
53                u2 * (c4 + u * c5 +
54                u2 * (cs[0] +
55                u2 * (cs[1] +
56                u2 * (cs[2] +
57                u2 * (cs[3] +
58                u2 * (cs[4] +
59                u2 * (cs[5]))))))))
60            } else if nz <= 1.0 {
61                cli5_unit_circle(-(1.0 - self).cln())
62            } else { // nz > 1.0
63                let pi4  = pi2*pi2;
64                let arg = if pz > 0.0 { pz - pi } else { pz + pi };
65                let lmz = Complex::new(lnz, arg); // (-self).cln()
66                let lmz2 = lmz*lmz;
67                cli5_unit_circle(-(1.0 - 1.0/self).cln()) - 1.0/360.0*lmz*(7.0*pi4 + lmz2*(10.0*pi2 + 3.0*lmz2))
68            }
69        }
70    }
71}
72
73/// series approximation of Li5(z) for |z| <= 1
74/// in terms of x = -ln(1 - z)
75fn cli5_unit_circle(x: Complex<f64>) -> Complex<f64> {
76    let bf  = [
77        1.0                   , -15.0/32.0             ,
78        1.3953189300411523e-01, -2.8633777006172840e-02,
79        4.0317412551440329e-03, -3.3985018004115226e-04,
80        4.5445184621617666e-06,  2.3916808048569012e-06,
81       -1.2762692600122747e-07, -3.1628984306505932e-08,
82        3.2848118445335192e-09,  4.7613713995660579e-10,
83       -8.0846898171909830e-11, -7.2387648587737207e-12,
84        1.9439760115173968e-12,  1.0256978405977236e-13,
85       -4.6180551009884830e-14, -1.1535857196470580e-15,
86        1.0903545401333394e-15
87    ];
88
89    let x2 = x*x;
90    let x4 = x2*x2;
91    let x8 = x4*x4;
92
93    x*bf[0] +
94    x2*(bf[1] + x*bf[2]) +
95    x4*(bf[3] + x*bf[4] + x2*(bf[5] + x*bf[6])) +
96    x8*(bf[7] + x*bf[8] + x2*(bf[9] + x*bf[10]) +
97        x4*(bf[11] + x*bf[12] + x2*(bf[13] + x*bf[14]))) +
98    x8*x8*(bf[15] + x*bf[16] + x2*(bf[17] + x*bf[18]))
99}
polylog/li5.rs

polylog/
li5.rs
