polylog/li/
rli.rs

1use num::complex::Complex;
2use crate::cln::CLn;
3use crate::{Li0, Li1, Li2, Li3, Li4};
4use super::eta::neg_eta;
5use super::fac::{fac, inv_fac};
6use super::harmonic::harmonic;
7use super::zeta::zeta;
8
9/// returns real n-th order polylogarithm Re[Li(n,x)] for real x
10pub fn rli(n: i32, x: f64) -> f64 {
11    let odd_sgn = |n| if is_even(n) { -1.0 } else { 1.0 };
12
13    if x == 0.0 {
14        x
15    } else if x == 1.0 {
16        zeta(n)
17    } else if x == -1.0 {
18        neg_eta(n)
19    } else if x.is_nan() {
20        std::f64::NAN
21    } else if n < -1 {
22        // arXiv:2010.09860
23        let c = 4.0*std::f64::consts::PI*std::f64::consts::PI;
24        let l2 = ln_sqr(x);
25        if c*x*x < l2 {
26            li_series(n, x)
27        } else if l2 < 0.512*0.512*c {
28            li_unity_neg(n, Complex::new(x, 0.0)).re
29        } else {
30            odd_sgn(n)*li_series(n, x.recip())
31        }
32    } else if n == -1 {
33        x/((1.0 - x)*(1.0 - x))
34    } else if n == 0 {
35        x.li0()
36    } else if n == 1 {
37        x.li1()
38    } else if n == 2 {
39        x.li2()
40    } else if n == 3 {
41        x.li3()
42    } else if n == 4 {
43        x.li4()
44    } else {
45        // transform x to y in [-1,1]
46        let (y, rest, sgn) = if x < -1.0 {
47            (x.recip(), li_neg_rest(n, x), odd_sgn(n))
48        } else if x < 1.0 {
49            (x, 0.0, 1.0)
50        } else { // x > 1.0
51            (x.recip(), li_pos_rest(n, x), odd_sgn(n))
52        };
53
54        if n < 20 && y > 0.75 {
55            sgn*li_unity_pos(n, y) + rest
56        } else {
57            sgn*li_series(n, y) + rest
58        }
59    }
60}
61
62/// returns true if x is even, false otherwise
63fn is_even(x: i32) -> bool {
64    x & 1 == 0
65}
66
67/// returns |ln(x)|^2 for all x
68fn ln_sqr(x: f64) -> f64 {
69    if x < 0.0 {
70        let l = (-x).ln();
71        l*l + std::f64::consts::PI*std::f64::consts::PI
72    } else if x == 0.0 {
73        std::f64::NAN
74    } else {
75        let l = x.ln();
76        l*l
77    }
78}
79
80#[test]
81fn test_ln_sqr() {
82    assert!(ln_sqr(2.0) == 2.0_f64.ln()*2.0_f64.ln());
83    assert!(ln_sqr(0.0).is_nan());
84    assert!(ln_sqr(-2.0) == Complex::new(-2.0, 0.0).ln().norm_sqr());
85}
86
87/// returns r.h.s. of inversion formula for x < -1:
88///
89/// Li(n,-x) + (-1)^n Li(n,-1/x)
90///    = -ln(n,x)^n/n! + 2 sum(r=1:(n÷2), ln(x)^(n-2r)/(n-2r)! Li(2r,-1))
91fn li_neg_rest(n: i32, x: f64) -> f64 {
92    let l = (-x).ln();
93    let l2 = l*l;
94
95    if is_even(n) {
96        let mut sum = 0.0;
97        let mut p = 1.0; // collects l^(2u)
98        for u in 0..n/2 {
99            let old_sum = sum;
100            sum += p*inv_fac(2*u)*neg_eta(n - 2*u);
101            p *= l2;
102            if sum == old_sum { break; }
103        }
104        2.0*sum - p*inv_fac(n)
105    } else {
106        let mut sum = 0.0;
107        let mut p = l; // collects l^(2u + 1)
108        for u in 0..(n - 1)/2 {
109            let old_sum = sum;
110            sum += p*inv_fac(2*u + 1)*neg_eta(n - 1 - 2*u);
111            p *= l2;
112            if sum == old_sum { break; }
113        }
114        2.0*sum - p*inv_fac(n)
115    }
116}
117
118/// returns (sin((2n+1)x), cos((2n+1)x)), given
119/// (sn, cn) = (sin(2nx), cos(2nx))   (previous value)
120/// (s2, c2) = (sin(2x), sin(2x))     (initial value)
121fn next_cosi((sn, cn): (f64, f64), (s2, c2): (f64, f64)) -> (f64, f64) {
122    (sn*c2 + cn*s2, cn*c2 - sn*s2)
123}
124
125/// returns r.h.s. of inversion formula for x > 1;
126/// same expression as in li_neg_rest(n,x), but with
127/// complex logarithm ln(-x)
128fn li_pos_rest(n: i32, x: f64) -> f64 {
129    let pi = std::f64::consts::PI;
130    let l = x.ln();
131    let mag = l.hypot(pi); // |ln(-x)|
132    let arg = pi.atan2(l); // arg(ln(-x))
133    let l2 = mag*mag;      // |ln(-x)|^2
134
135    if is_even(n) {
136        let mut sum = 0.0;
137        let mut p = 1.0; // collects mag^(2u)
138        let mut cosi = (0.0, 1.0); // collects (sin(2*u*arg), cos(2*u*arg))
139        let cosi2 = (2.0*arg).sin_cos();
140        for u in 0..n/2 {
141            let old_sum = sum;
142            sum += p*cosi.1*inv_fac(2*u)*neg_eta(n - 2*u);
143            if sum == old_sum { break; }
144            p *= l2;
145            cosi = next_cosi(cosi, cosi2);
146        }
147        2.0*sum - p*cosi.1*inv_fac(n)
148    } else {
149        let mut sum = 0.0;
150        let mut p = mag; // collects mag^(2u + 1)
151        let (s, c) = arg.sin_cos();
152        let mut cosi = (s, c); // collects (sin((2*u + 1)*arg), cos((2*u + 1)*arg))
153        let cosi2 = (2.0*s*c, 2.0*c*c - 1.0); // (2.0*arg).sin_cos()
154        for u in 0..(n - 1)/2 {
155            let old_sum = sum;
156            sum += p*cosi.1*inv_fac(2*u + 1)*neg_eta(n - 1 - 2*u);
157            if sum == old_sum { break; }
158            p *= l2;
159            cosi = next_cosi(cosi, cosi2);
160        }
161        2.0*sum - p*cosi.1*inv_fac(n)
162    }
163}
164
165/// returns Li(n,x) using the series expansion of Li(n,x) for n > 0
166/// and real x ~ 1 where 0 < x < 1:
167///
168/// Li(n,x) = sum(j=0:Inf, zeta(n-j) ln(x)^j/j!)
169///
170/// where
171///
172/// zeta(1) = -ln(-ln(x)) + harmonic(n - 1)
173///
174/// harmonic(n) = sum(k=1:n, 1/k)
175fn li_unity_pos(n: i32, x: f64) -> f64 {
176    let l = x.ln();
177    let mut sum = zeta(n);
178    let mut p = 1.0; // collects l^j/j!
179
180    for j in 1..(n - 1) {
181        p *= l/(j as f64);
182        sum += zeta(n - j)*p;
183    }
184
185    p *= l/((n - 1) as f64);
186    sum += (harmonic(n - 1) - (-l).ln())*p;
187
188    p *= l/(n as f64);
189    sum += zeta(0)*p;
190
191    p *= l/((n + 1) as f64);
192    sum += zeta(-1)*p;
193
194    let l2 = l*l;
195
196    for j in ((n + 3)..i32::MAX).step_by(2) {
197        p *= l2/(((j - 1)*j) as f64);
198        let old_sum = sum;
199        sum += zeta(n - j)*p;
200        if sum == old_sum { break; }
201    }
202
203    sum
204}
205
206/// returns Li(n,x) using the series expansion for n < 0 and x ~ 1
207///
208/// Li(n,x) = gamma(1-n) (-ln(x))^(n-1)
209///           + sum(k=0:Inf, zeta(n-k) ln(x)^k/k!)
210fn li_unity_neg(n: i32, z: Complex<f64>) -> Complex<f64> {
211    let lnz = z.cln();
212    let lnz2 = lnz*lnz;
213    let mut sum = fac(-n)*(-lnz).powi(n - 1);
214    let (mut k, mut lnzk) = if is_even(n) {
215        (1, lnz)
216    } else {
217        sum += zeta(n);
218        (2, lnz2)
219    };
220
221    loop {
222        let term = zeta(n - k)*inv_fac(k)*lnzk;
223        if !term.is_finite() { break; }
224        let sum_old = sum;
225        sum += term;
226        if sum == sum_old || k >= i32::MAX - 2 { break; }
227        lnzk *= lnz2;
228        k += 2;
229    }
230
231    sum
232}
233
234/// returns Li(n,x) using the naive series expansion of Li(n,x)
235/// for |x| < 1:
236///
237/// Li(n,x) = sum(k=1:Inf, x^k/k^n)
238fn li_series(n: i32, x: f64) -> f64
239{
240    let mut sum = x;
241    let mut xn = x*x;
242
243    for k in 2..i32::MAX {
244        let term = xn/(k as f64).powi(n);
245        if !term.is_finite() { break; }
246        let old_sum = sum;
247        sum += term;
248        if sum == old_sum { break; }
249        xn *= x;
250    }
251
252    sum
253}
polylog/li/rli.rs

polylog/li/
rli.rs
