polylog/
li2.rs

1use num::complex::Complex;
2use crate::cln::CLn;
3
4trait Li2Approx<T> {
5    fn approx(&self) -> T;
6}
7
8impl Li2Approx<f32> for f32 {
9    /// rational function approximation of Re[Li2(x)] for x in [0, 1/2]
10    fn approx(&self) -> f32 {
11        let cp = [ 1.00000020_f32, -0.780790946_f32, 0.0648256871_f32 ];
12        let cq = [ 1.00000000_f32, -1.03077545_f32, 0.211216710_f32 ];
13
14        let x = *self;
15        let p = cp[0] + x*(cp[1] + x*cp[2]);
16        let q = cq[0] + x*(cq[1] + x*cq[2]);
17
18        x*p/q
19    }
20}
21
22impl Li2Approx<f64> for f64 {
23    /// rational function approximation of Re[Li2(x)] for x in [0, 1/2]
24    fn approx(&self) -> f64 {
25        let cp = [
26            0.9999999999999999502e+0_f64,
27           -2.6883926818565423430e+0_f64,
28            2.6477222699473109692e+0_f64,
29           -1.1538559607887416355e+0_f64,
30            2.0886077795020607837e-1_f64,
31           -1.0859777134152463084e-2_f64
32        ];
33        let cq = [
34            1.0000000000000000000e+0_f64,
35           -2.9383926818565635485e+0_f64,
36            3.2712093293018635389e+0_f64,
37           -1.7076702173954289421e+0_f64,
38            4.1596017228400603836e-1_f64,
39           -3.9801343754084482956e-2_f64,
40            8.2743668974466659035e-4_f64
41        ];
42
43        let x = *self;
44        let x2 = x*x;
45        let x4 = x2*x2;
46        let p = cp[0] + x*cp[1] + x2*(cp[2] + x*cp[3]) +
47            x4*(cp[4] + x*cp[5]);
48        let q = cq[0] + x*cq[1] + x2*(cq[2] + x*cq[3]) +
49            x4*(cq[4] + x*cq[5] + x2*cq[6]);
50
51        x*p/q
52    }
53}
54
55impl Li2Approx<Complex<f32>> for Complex<f32> {
56    /// series approximation of Li2(z) for Re(z) <= 1/2 and |z| <= 1
57    /// in terms of self = -ln(1 - z)
58    fn approx(&self) -> Complex<f32> {
59        // bf[1..N-1] are the even Bernoulli numbers / (2 n + 1)!
60        // generated by: Table[BernoulliB[2 n]/(2 n + 1)!, {n, 1, 19}]
61        let bf = [
62            -1.0_f32/4.0_f32,
63             1.0_f32/36.0_f32,
64            -1.0_f32/3600.0_f32,
65             1.0_f32/211680.0_f32,
66        ];
67        let x = *self;
68        let x2 = x*x;
69
70        x + x2*(bf[0] + x*(bf[1] + x2*(bf[2] + x2*bf[3])))
71    }
72}
73
74impl Li2Approx<Complex<f64>> for Complex<f64> {
75    /// series approximation of Li2(z) for Re(z) <= 1/2 and |z| <= 1
76    /// in terms of self = -ln(1 - z)
77    fn approx(&self) -> Complex<f64> {
78        // bf[1..N-1] are the even Bernoulli numbers / (2 n + 1)!
79        // generated by: Table[BernoulliB[2 n]/(2 n + 1)!, {n, 1, 19}]
80        let bf = [
81            -1.0_f64/4.0_f64,
82             1.0_f64/36.0_f64,
83            -1.0_f64/3600.0_f64,
84             1.0_f64/211680.0_f64,
85            -1.0_f64/10886400.0_f64,
86             1.0_f64/526901760.0_f64,
87            -4.0647616451442255e-11_f64,
88             8.9216910204564526e-13_f64,
89            -1.9939295860721076e-14_f64,
90             4.5189800296199182e-16_f64,
91        ];
92        let x = *self;
93        let x2 = x*x;
94        let x4 = x2*x2;
95
96        x + x2*(bf[0] +
97                x*(bf[1] +
98                   x2*(
99                       bf[2] +
100                       x2*bf[3] +
101                       x4*(bf[4] + x2*bf[5]) +
102                       x4*x4*(bf[6] + x2*bf[7] + x4*(bf[8] + x2*bf[9]))
103                   )
104                )
105        )
106    }
107}
108
109/// Provides the 2nd order polylogarithm (dilogarithm) function
110/// `li2()` of a number of type `T`.
111pub trait Li2<T> {
112    fn li2(&self) -> T;
113}
114
115impl Li2<f32> for f32 {
116    /// Returns the real dilogarithm of a real number of type `f32`.
117    ///
118    /// Implemented as rational function approximation.
119    ///
120    /// # Example:
121    /// ```
122    /// use polylog::Li2;
123    ///
124    /// assert!((1.0_f32.li2() - 1.64493407_f32).abs() < 2.0_f32*std::f32::EPSILON);
125    /// ```
126    fn li2(&self) -> f32 {
127        let z2 = std::f32::consts::PI*std::f32::consts::PI/6.0_f32;
128        let x = *self;
129
130        // transform to [0, 1/2]
131        if x < -1.0_f32 {
132            let l = (1.0_f32 - x).ln();
133            (1.0_f32/(1.0_f32 - x)).approx() - z2 + l*(0.5_f32*l - (-x).ln())
134        } else if x == -1.0_f32 {
135            -0.5_f32*z2
136        } else if x < 0.0_f32 {
137            let l = (-x).ln_1p();
138            -(x/(x - 1.0_f32)).approx() - 0.5_f32*l*l
139        } else if x == 0.0_f32 {
140            x
141        } else if x < 0.5_f32 {
142            x.approx()
143        } else if x < 1.0_f32 {
144            -(1.0_f32 - x).approx() + z2 - x.ln()*(-x).ln_1p()
145        } else if x == 1.0_f32 {
146            z2
147        } else if x < 2.0_f32 {
148            let l = x.ln();
149            (1.0_f32 - 1.0_f32/x).approx() + z2 - l*((1.0_f32 - 1.0_f32/x).ln() + 0.5_f32*l)
150        } else {
151            let l = x.ln();
152            -(1.0_f32/x).approx() + 2.0_f32*z2 - 0.5_f32*l*l
153        }
154    }
155}
156
157impl Li2<f64> for f64 {
158    /// Returns the real dilogarithm of a real number of type `f64`.
159    ///
160    /// Implemented as rational function approximation with a maximum
161    /// error of 5e-17 [[arXiv:2201.01678]].
162    ///
163    /// [arXiv:2201.01678]: https://arxiv.org/abs/2201.01678
164    ///
165    /// # Example:
166    /// ```
167    /// use polylog::Li2;
168    ///
169    /// assert!((1.0_f64.li2() - 1.6449340668482264_f64).abs() < 2.0_f64*std::f64::EPSILON);
170    /// ```
171    fn li2(&self) -> f64 {
172        let z2 = std::f64::consts::PI*std::f64::consts::PI/6.0_f64;
173        let x = *self;
174
175        // transform to [0, 1/2]
176        if x < -1.0_f64 {
177            let l = (1.0_f64 - x).ln();
178            (1.0_f64/(1.0_f64 - x)).approx() - z2 + l*(0.5_f64*l - (-x).ln())
179        } else if x == -1.0_f64 {
180            -0.5_f64*z2
181        } else if x < 0.0_f64 {
182            let l = (-x).ln_1p();
183            -(x/(x - 1.0_f64)).approx() - 0.5_f64*l*l
184        } else if x == 0.0_f64 {
185            x
186        } else if x < 0.5_f64 {
187            x.approx()
188        } else if x < 1.0_f64 {
189            -(1.0_f64 - x).approx() + z2 - x.ln()*(-x).ln_1p()
190        } else if x == 1.0_f64 {
191            z2
192        } else if x < 2.0_f64 {
193            let l = x.ln();
194            (1.0_f64 - 1.0_f64/x).approx() + z2 - l*((1.0_f64 - 1.0_f64/x).ln() + 0.5_f64*l)
195        } else {
196            let l = x.ln();
197            -(1.0_f64/x).approx() + 2.0_f64*z2 - 0.5_f64*l*l
198        }
199    }
200}
201
202impl Li2<Complex<f32>> for Complex<f32> {
203    /// Returns the dilogarithm of a complex number of type
204    /// `Complex<f32>`.
205    ///
206    /// # Example:
207    /// ```
208    /// use num::complex::Complex;
209    /// use polylog::Li2;
210    ///
211    /// assert!((Complex::new(1.0_f32, 1.0_f32).li2() - Complex::new(0.61685028_f32, 1.46036212_f32)).norm() < std::f32::EPSILON);
212    /// ```
213    fn li2(&self) -> Complex<f32> {
214        let pi = std::f32::consts::PI;
215        let rz = self.re;
216        let iz = self.im;
217
218        if iz == 0.0_f32 {
219            if rz <= 1.0_f32 {
220                Complex::new(rz.li2(), iz)
221            } else { // rz > 1
222                Complex::new(rz.li2(), -pi*rz.ln())
223            }
224        } else {
225            let nz = self.norm_sqr();
226
227            if nz < std::f32::EPSILON {
228                self*(1.0_f32 + 0.25_f32*self)
229            } else if rz <= 0.5_f32 {
230                if nz > 1.0_f32 {
231                    let l = (-self).cln();
232                    -(-(1.0_f32 - 1.0_f32/self).cln()).approx() - 0.5_f32*l*l - pi*pi/6.0_f32
233                } else { // nz <= 1
234                    (-(1.0_f32 - self).cln()).approx()
235                }
236            } else { // rz > 0.5
237                if nz <= 2.0_f32*rz {
238                    let l = -(self).cln();
239                    -l.approx() + l*(1.0_f32 - self).cln() + pi*pi/6.0_f32
240                } else { // nz > 2*rz
241                    let l = (-self).cln();
242                    -(-(1.0_f32 - 1.0_f32/self).cln()).approx() - 0.5_f32*l*l - pi*pi/6.0_f32
243                }
244            }
245        }
246    }
247}
248
249impl Li2<Complex<f64>> for Complex<f64> {
250    /// Returns the dilogarithm of a complex number of type
251    /// `Complex<f64>`.
252    ///
253    /// This function has been translated from the
254    /// [SPheno](https://spheno.hepforge.org/) package.
255    ///
256    /// # Example:
257    /// ```
258    /// use num::complex::Complex;
259    /// use polylog::Li2;
260    ///
261    /// assert!((Complex::new(1.0_f64, 1.0_f64).li2() - Complex::new(0.6168502750680849_f64, 1.4603621167531195_f64)).norm() < 2.0_f64*std::f64::EPSILON);
262    /// ```
263    fn li2(&self) -> Complex<f64> {
264        let pi = std::f64::consts::PI;
265        let rz = self.re;
266        let iz = self.im;
267
268        if iz == 0.0_f64 {
269            if rz <= 1.0_f64 {
270                Complex::new(rz.li2(), iz)
271            } else { // rz > 1
272                Complex::new(rz.li2(), -pi*rz.ln())
273            }
274        } else {
275            let nz = self.norm_sqr();
276
277            if nz < std::f64::EPSILON {
278                self*(1.0_f64 + 0.25_f64*self)
279            } else if rz <= 0.5_f64 {
280                if nz > 1.0_f64 {
281                    let l = (-self).cln();
282                    -(-(1.0_f64 - 1.0_f64/self).cln()).approx() - 0.5_f64*l*l - pi*pi/6.0_f64
283                } else { // nz <= 1
284                    (-(1.0_f64 - self).cln()).approx()
285                }
286            } else { // rz > 0.5
287                if nz <= 2.0_f64*rz {
288                    let l = -(self).cln();
289                    -l.approx() + l*(1.0_f64 - self).cln() + pi*pi/6.0_f64
290                } else { // nz > 2*rz
291                    let l = (-self).cln();
292                    -(-(1.0_f64 - 1.0_f64/self).cln()).approx() - 0.5_f64*l*l - pi*pi/6.0_f64
293                }
294            }
295        }
296    }
297}
polylog/li2.rs

polylog/
li2.rs
