polylog/
li4.rs

1use num::complex::Complex;
2use crate::cln::CLn;
3
4/// Provides the 4-th order polylogarithm function `li4()` of a
5/// number of type `T`.
6pub trait Li4<T> {
7    fn li4(&self) -> T;
8}
9
10impl Li4<f64> for f64 {
11    /// Returns the fourth order polylogarithm of a real number of type `f64`.
12    ///
13    /// # Example:
14    /// ```
15    /// use polylog::Li4;
16    ///
17    /// assert!((1.0_f64.li4() - 1.0823232337111382_f64).abs() < std::f64::EPSILON);
18    /// ```
19    fn li4(&self) -> f64 {
20        let z2 = 1.6449340668482264;
21        let z4 = 1.0823232337111382;
22        let x = *self;
23
24        // transform x to y in [-1,1]
25        let (y, rest, sgn) = if x < -1.0 {
26            let l = (-x).ln();
27            let l2 = l*l;
28            (1.0/x, -7.0/4.0*z4 + l2*(-0.5*z2 - 1.0/24.0*l2), -1.0)
29        } else if x == -1.0 {
30            return -7.0/8.0*z4
31        } else if x == 0.0 {
32            return x
33        } else if x < 1.0 {
34            (x, 0.0, 1.0)
35        } else if x == 1.0 {
36            return z4
37        } else { // x > 1.0
38            let l = x.ln();
39            let l2 = l*l;
40            (1.0/x, 2.0*z4 + l2*(z2 - 1.0/24.0*l2), -1.0)
41        };
42
43        if y < 0.0 {
44            sgn*li4_neg(y) + rest
45        } else if y < 0.5 {
46            sgn*li4_half(y) + rest
47        } else if y < 0.8 {
48            sgn*li4_mid(y) + rest
49        } else { // y <= 1.0
50            sgn*li4_one(y) + rest
51        }
52    }
53}
54
55// Li_4(x) for x in [-1,0]
56fn li4_neg(x: f64) -> f64 {
57    let cp = [
58        0.9999999999999999952e+0, -1.8532099956062184217e+0,
59        1.1937642574034898249e+0, -3.1817912243893560382e-1,
60        3.2268284189261624841e-2, -8.3773570305913850724e-4
61    ];
62    let cq = [
63        1.0000000000000000000e+0, -1.9157099956062165688e+0,
64        1.3011504531166486419e+0, -3.7975653506939627186e-1,
65        4.5822723996558783670e-2, -1.8023912938765272341e-3,
66        1.0199621542882314929e-5
67    ];
68
69    let x2 = x*x;
70    let x4 = x2*x2;
71    let p = cp[0] + x * cp[1] + x2 * (cp[2] + x * cp[3]) +
72            x4 * (cp[4] + x * cp[5]);
73    let q = cq[0] + x * cq[1] + x2 * (cq[2] + x * cq[3]) +
74            x4 * (cq[4] + x * cq[5] + x2 * cq[6]);
75
76    x*p/q
77}
78
79// Li_4(x) for x in [0,1/2]
80fn li4_half(x: f64) -> f64 {
81    let cp = [
82        1.0000000000000000414e+0, -2.0588072418045364525e+0,
83        1.4713328756794826579e+0, -4.2608608613069811474e-1,
84        4.2975084278851543150e-2, -6.8314031819918920802e-4
85    ];
86    let cq = [
87        1.0000000000000000000e+0, -2.1213072418045207223e+0,
88        1.5915688992789175941e+0, -5.0327641401677265813e-1,
89        6.1467217495127095177e-2, -1.9061294280193280330e-3
90    ];
91
92    let x2 = x*x;
93    let x4 = x2*x2;
94    let p = cp[0] + x * cp[1] + x2 * (cp[2] + x * cp[3]) +
95            x4 * (cp[4] + x * cp[5]);
96    let q = cq[0] + x * cq[1] + x2 * (cq[2] + x * cq[3]) +
97            x4 * (cq[4] + x * cq[5]);
98
99    x*p/q
100}
101
102// Li_4(x) for x in [1/2,8/10]
103fn li4_mid(x: f64) -> f64 {
104    let cp = [
105        3.2009826406098890447e-9, 9.9999994634837574160e-1,
106       -2.9144851228299341318e+0, 3.1891031447462342009e+0,
107       -1.6009125158511117090e+0, 3.5397747039432351193e-1,
108       -2.5230024124741454735e-2
109    ];
110    let cq = [
111        1.0000000000000000000e+0, -2.9769855248411488460e+0,
112        3.3628208295110572579e+0, -1.7782471949702788393e+0,
113        4.3364007973198649921e-1, -3.9535592340362510549e-2,
114        5.7373431535336755591e-4
115    ];
116
117    let x2 = x*x;
118    let x4 = x2*x2;
119    let p = cp[0] + x * cp[1] + x2 * (cp[2] + x * cp[3]) +
120            x4 * (cp[4] + x * cp[5] + x2 * cp[6]);
121    let q = cq[0] + x * cq[1] + x2 * (cq[2] + x * cq[3]) +
122            x4 * (cq[4] + x * cq[5] + x2 * cq[6]);
123
124    p/q
125}
126
127// Li_4(x) for x in [8/10,1]
128fn li4_one(x: f64) -> f64 {
129    let z2 = 1.6449340668482264;
130    let z3 = 1.2020569031595943;
131    let z4 = 1.0823232337111382;
132    let l = x.ln();
133    let l2 = l*l;
134
135    z4 +
136    l*(z3 +
137    l*(0.5*z2 +
138    l*(11.0/36.0 - 1.0/6.0*(-l).ln() +
139    l*(-1.0/48.0 +
140    l*(-1.0/1440.0 +
141    l2*(1.0/604800.0 - 1.0/91445760.0*l2))))))
142}
143
144impl Li4<Complex<f64>> for Complex<f64> {
145    /// Returns the fourth order polylogarithm of a complex number of type
146    /// `Complex<f64>`.
147    ///
148    /// # Example:
149    /// ```
150    /// use num::complex::Complex;
151    /// use polylog::Li4;
152    ///
153    /// assert!((Complex::new(1.0_f64, 1.0_f64).li4() - Complex::new(0.9593189135784193_f64, 1.1380391966769828_f64)).norm() < 2.0_f64*std::f64::EPSILON);
154    /// ```
155    fn li4(&self) -> Complex<f64> {
156        let pi  = std::f64::consts::PI;
157        let pi2 = pi*pi;
158        let z4  = 1.0823232337111382;
159
160        if self.im == 0.0 {
161            if self.re <= 1.0 {
162                Complex::new(self.re.li4(), self.im)
163            } else { // rz > 1.0
164                let l = self.re.ln();
165                Complex::new(self.re.li4(), -pi/6.0*l*l*l)
166            }
167        } else {
168            let nz  = self.norm();
169            let pz  = self.arg();
170            let lnz = nz.ln();
171
172            if lnz*lnz + pz*pz < 1.0 { // |log(z)| < 1
173                let u  = Complex::new(lnz, pz);
174                let u2 = u*u;
175                let u4 = u2*u2;
176                let u8 = u4*u4;
177                let c1 = 1.2020569031595943; // zeta(3)
178                let c2 = 0.82246703342411322;
179                let c3 = (11.0/6.0 - (-u).cln())/6.0;
180                let c4 = -1.0/48.0;
181
182                let cs = [
183                    -6.9444444444444444e-04, 1.6534391534391534e-06,
184                    -1.0935444136502338e-08, 1.0438378493934049e-10,
185                    -1.2165942300622435e-12, 1.6130006528350101e-14,
186                    -2.3428810452879340e-16
187                ];
188
189                z4 + u2*(c2 + u2*c4) +
190                u*(c1 +
191                   c3*u2 +
192                   u4*(cs[0] + u2*cs[1]) +
193                   u8*(cs[2] + u2*cs[3] + u4*(cs[4] + u2*cs[5])) +
194                   u8*u8*cs[6]
195                )
196            } else if nz <= 1.0 {
197                cli4_unit_circle(-(1.0 - self).cln())
198            } else { // nz > 1.0
199                let pi4  = pi2*pi2;
200                let arg = if pz > 0.0 { pz - pi } else { pz + pi };
201                let lmz = Complex::new(lnz, arg); // (-self).cln()
202                let lmz2 = lmz*lmz;
203                -cli4_unit_circle(-(1.0 - 1.0/self).cln()) + 1.0/360.0*(-7.0*pi4 + lmz2*(-30.0*pi2 - 15.0*lmz2))
204            }
205        }
206    }
207}
208
209/// series approximation of Li4(z) for |z| <= 1
210/// in terms of x = -ln(1 - z)
211fn cli4_unit_circle(x: Complex<f64>) -> Complex<f64> {
212    let bf  = [
213        1.0                   , -7.0/16.0              ,
214        1.1651234567901235e-01, -1.9820601851851852e-02,
215        1.9279320987654321e-03, -3.1057098765432099e-05,
216       -1.5624009114857835e-05,  8.4851235467732066e-07,
217        2.2909616603189711e-07, -2.1832614218526917e-08,
218       -3.8828248791720156e-09,  5.4462921032203321e-10,
219        6.9608052106827254e-11, -1.3375737686445215e-11,
220       -1.2784852685266572e-12,  3.2605628580248922e-13,
221        2.3647571168618257e-14, -7.9231351220311617e-15,
222    ];
223
224    let x2 = x*x;
225    let x4 = x2*x2;
226    let x8 = x4*x4;
227
228    x*bf[0] +
229    x2*(bf[1] + x*bf[2]) +
230    x4*(bf[3] + x*bf[4] + x2*(bf[5] + x*bf[6])) +
231    x8*(bf[7] + x*bf[8] + x2*(bf[9] + x*bf[10]) +
232        x4*(bf[11] + x*bf[12] + x2*(bf[13] + x*bf[14]))) +
233    x8*x8*(bf[15] + x*bf[16] + x2*bf[17])
234}
polylog/li4.rs

polylog/
li4.rs
