1 /*
2     Copyright (C) 2014 Fredrik Johansson
3 
4     This file is part of Arb.
5 
6     Arb is free software: you can redistribute it and/or modify it under
7     the terms of the GNU Lesser General Public License (LGPL) as published
8     by the Free Software Foundation; either version 2.1 of the License, or
9     (at your option) any later version.  See <http://www.gnu.org/licenses/>.
10 */
11 
12 #include "acb_poly.h"
13 #include "acb_hypgeom.h"
14 
15 /* note: will not return a wrong value, as arf_get_si aborts on overflow */
16 slong
arb_get_si_lower(const arb_t x)17 arb_get_si_lower(const arb_t x)
18 {
19     arf_t t;
20     slong v;
21 
22     arf_init(t);
23     arf_set_mag(t, arb_radref(x));
24     arf_sub(t, arb_midref(x), t, 2 * FLINT_BITS, ARF_RND_FLOOR);
25 
26     v = arf_get_si(t, ARF_RND_FLOOR);
27 
28     arf_clear(t);
29 
30     return v;
31 }
32 
33 slong
polylog_choose_terms(mag_t err,slong sigma,const mag_t z,slong d,slong prec)34 polylog_choose_terms(mag_t err, slong sigma, const mag_t z, slong d, slong prec)
35 {
36     slong N;
37 
38     for (N = 3; ; N = FLINT_MAX(N+3, N*1.1))
39     {
40         mag_polylog_tail(err, z, sigma, d, N);
41 
42         /* TODO: do something else when |Li_s(z)| is very small/very large? */
43         if (mag_cmp_2exp_si(err, -prec) < 0)
44             break;
45 
46         if (N > 100 * prec)
47         {
48             N = 3;
49             mag_inf(err);
50             break;
51         }
52     }
53 
54     return N;
55 }
56 
57 int
polylog_is_real(const acb_t s,const acb_t z)58 polylog_is_real(const acb_t s, const acb_t z)
59 {
60     if (!arb_is_zero(acb_imagref(s)))
61         return 0;
62     else if (!arb_is_zero(acb_imagref(z)))
63         return 0;
64     else if (arb_contains_si(acb_realref(z), 1))
65         return 0;
66     else if (acb_is_int(s) && arb_is_nonpositive(acb_realref(s)))
67         return 1;
68     else
69         return (arf_cmp_2exp_si(arb_midref(acb_realref(z)), 0) < 0);
70 }
71 
72 void
_acb_poly_polylog_cpx_zeta(acb_ptr w,const acb_t s,const acb_t z,slong len,slong prec)73 _acb_poly_polylog_cpx_zeta(acb_ptr w, const acb_t s, const acb_t z, slong len, slong prec)
74 {
75     acb_ptr e1, e2, z1, z2, e1z1, e2z2;
76     acb_t t, u, v;
77     slong k, len2;
78     int deflate_zeta, deflate_gamma, is_real;
79 
80     if (!acb_is_finite(s) || !acb_is_finite(z))
81     {
82         _acb_vec_indeterminate(w, len);
83         return;
84     }
85 
86     if (acb_is_one(z))
87     {
88         if (arb_gt(acb_realref(s), acb_realref(z))) /* Re(s) > 1 */
89         {
90             acb_zeta(w, s, prec);
91             _acb_vec_indeterminate(w + 1, len - 1);
92         }
93         else
94         {
95             _acb_vec_indeterminate(w, len);
96         }
97 
98         return;
99     }
100 
101     is_real = polylog_is_real(s, z);
102 
103     acb_init(t);
104     acb_init(u);
105     acb_init(v);
106 
107     /* v = 1-s */
108     acb_one(v);
109     acb_sub(v, v, s, prec);
110 
111     /* pole of zeta */
112     deflate_zeta = acb_is_one(v);
113 
114     /* poles of gamma at nonpositive integer v */
115     deflate_gamma = (arb_is_zero(acb_imagref(v)) &&
116             arb_is_int(acb_realref(v)) &&
117             arf_sgn(arb_midref(acb_realref(v))) <= 0);
118 
119     len2 = len + deflate_gamma;
120 
121     e1 = _acb_vec_init(len + 1);
122     e2 = _acb_vec_init(len + 1);
123     z1 = _acb_vec_init(len + 1);
124     z2 = _acb_vec_init(len + 1);
125     e1z1 = _acb_vec_init(len + 1);
126     e2z2 = _acb_vec_init(len + 1);
127 
128     /* u = log(-z)/(pi*i) */
129     acb_neg(t, z);
130     acb_log(t, t, prec);
131     acb_const_pi(u, prec);
132     acb_mul_onei(u, u);
133     acb_div(u, t, u, prec);
134 
135     /* z1 = zeta(v+x, 1/2 + log(-z)/(2*pi*i)) */
136     acb_one(t);
137     acb_add(t, t, u, prec);
138     acb_mul_2exp_si(t, t, -1);
139     _acb_poly_zeta_cpx_series(z1, v, t, deflate_zeta, len2, prec);
140 
141     /* z2 = zeta(v+x, 1/2 - log(-z)/(2*pi*i)) */
142     acb_one(t);
143     acb_sub(t, t, u, prec);
144     acb_mul_2exp_si(t, t, -1);
145     _acb_poly_zeta_cpx_series(z2, v, t, deflate_zeta, len2, prec);
146 
147     /* e1 = (i/(2pi))^(v+x) */
148     acb_onei(t);
149     acb_const_pi(u, prec);
150     acb_div(t, t, u, prec);
151     acb_mul_2exp_si(t, t, -1);
152     _acb_poly_acb_pow_cpx(e1, t, v, len + (deflate_zeta || deflate_gamma), prec);
153 
154     /* e2 = (1/(2 pi i))^(v+x) */
155     acb_conj(t, t);
156     _acb_poly_acb_pow_cpx(e2, t, v, len + (deflate_zeta || deflate_gamma), prec);
157 
158     _acb_poly_mullow(e1z1, e1, len2, z1, len2, len2, prec);
159     _acb_poly_mullow(e2z2, e2, len2, z2, len2, len2, prec);
160     _acb_vec_add(z1, e1z1, e2z2, len2, prec);
161 
162     if (deflate_gamma)
163     {
164         /* gamma(v+x) = pi/sin(pi(v+x)) * 1/gamma(1-v-x) */
165 
166         /* TODO: write a csc function? */
167         acb_zero(e1);
168         acb_const_pi(e1 + 1, prec);
169         acb_mul_2exp_si(e2, v, -1);
170         if (!arb_is_int(acb_realref(e2)))
171             acb_neg(e1 + 1, e1 + 1);
172         _acb_poly_sin_series(e2, e1, 2, len2, prec);
173         _acb_poly_inv_series(e1, e2 + 1, len, len, prec);
174         acb_const_pi(e2, prec);
175         _acb_vec_scalar_mul(e1, e1, len, e2, prec);
176 
177         acb_set(z2, s);
178         acb_set_si(z2 + 1, -1);
179         _acb_poly_rgamma_series(e2, z2, 2, len, prec);
180         _acb_poly_mullow(z2, e1, len, e2, len, len, prec);
181 
182         _acb_poly_mullow(w, z1 + 1, len, z2, len, len, prec);
183     }
184     else
185     {
186         if (deflate_zeta)
187         {
188             for (k = 0; k < len; k++)
189             {
190                 arb_mul_2exp_si(acb_realref(e1 + k + 1), acb_realref(e1 + k + 1), 1);
191                 arb_add(acb_realref(z1 + k), acb_realref(z1 + k), acb_realref(e1 + k + 1), prec);
192             }
193 
194         }
195 
196         /* gamma(v+x) */
197         acb_set(e1, v);
198         if (len > 1)
199             acb_one(e1 + 1);
200         _acb_poly_gamma_series(z2, e1, FLINT_MIN(len, 2), len, prec);
201 
202         _acb_poly_mullow(w, z2, len, z1, len, len, prec);
203     }
204 
205     /* correct signs (from s -> 1-s) */
206     for (k = 1; k < len; k += 2)
207         acb_neg(w + k, w + k);
208 
209     if (is_real)
210         if (acb_is_finite(w))
211             arb_zero(acb_imagref(w));
212 
213     _acb_vec_clear(e1, len + 1);
214     _acb_vec_clear(e2, len + 1);
215     _acb_vec_clear(z1, len + 1);
216     _acb_vec_clear(z2, len + 1);
217     _acb_vec_clear(e1z1, len + 1);
218     _acb_vec_clear(e2z2, len + 1);
219 
220     acb_clear(t);
221     acb_clear(u);
222     acb_clear(v);
223 }
224 
225 void
_acb_poly_polylog_cpx_small(acb_ptr w,const acb_t s,const acb_t z,slong len,slong prec)226 _acb_poly_polylog_cpx_small(acb_ptr w, const acb_t s, const acb_t z, slong len, slong prec)
227 {
228     slong k, N, sigma;
229     int is_real;
230     mag_t zmag, err, errf;
231     acb_t a;
232 
233     acb_init(a);
234     mag_init(zmag);
235     mag_init(err);
236     mag_init(errf);
237 
238     is_real = polylog_is_real(s, z);
239     acb_get_mag(zmag, z);
240     sigma = arb_get_si_lower(acb_realref(s));
241 
242     N = polylog_choose_terms(err, sigma, zmag, len - 1, prec);
243 
244     /* TODO: allow threading */
245     acb_one(a);
246     _acb_poly_powsum_series_naive(w, s, a, z, N - 1, len, prec);
247     _acb_vec_scalar_mul(w, w, len, z, prec);
248 
249     for (k = 0; k < len; k++)
250     {
251         mag_polylog_tail(err, zmag, sigma, k, N);
252         mag_rfac_ui(errf, k);
253         mag_mul(err, err, errf);
254 
255         if (is_real && mag_is_finite(err))
256             arb_add_error_mag(acb_realref(w + k), err);
257         else
258             acb_add_error_mag(w + k, err);
259     }
260 
261     acb_clear(a);
262     mag_clear(zmag);
263     mag_clear(err);
264     mag_clear(errf);
265 }
266 
267 void
_acb_poly_polylog_cpx(acb_ptr w,const acb_t s,const acb_t z,slong len,slong prec)268 _acb_poly_polylog_cpx(acb_ptr w, const acb_t s, const acb_t z, slong len, slong prec)
269 {
270     mag_t zmag;
271 
272     if (len == 1 && acb_equal_si(s, 2))
273     {
274         acb_hypgeom_dilog(w, z, prec);
275         return;
276     }
277 
278     mag_init(zmag);
279     acb_get_mag(zmag, z);
280 
281     if (mag_cmp_2exp_si(zmag, -1) < 0)
282         _acb_poly_polylog_cpx_small(w, s, z, len, prec);
283     else
284         _acb_poly_polylog_cpx_zeta(w, s, z, len, prec);
285 
286     mag_clear(zmag);
287 }
288 
289 void
_acb_poly_polylog_series(acb_ptr res,acb_srcptr s,slong slen,const acb_t z,slong len,slong prec)290 _acb_poly_polylog_series(acb_ptr res, acb_srcptr s, slong slen, const acb_t z, slong len, slong prec)
291 {
292     acb_ptr t, u;
293 
294     slen = FLINT_MIN(slen, len);
295 
296     t = _acb_vec_init(len);
297     u = _acb_vec_init(len);
298 
299     _acb_poly_polylog_cpx(t, s, z, len, prec);
300 
301     /* compose with nonconstant part */
302     acb_zero(u);
303     _acb_vec_set(u + 1, s + 1, slen - 1);
304     _acb_poly_compose_series(res, t, len, u, slen, len, prec);
305 
306     _acb_vec_clear(t, len);
307     _acb_vec_clear(u, len);
308 }
309 
310 void
acb_poly_polylog_series(acb_poly_t res,const acb_poly_t s,const acb_t z,slong n,slong prec)311 acb_poly_polylog_series(acb_poly_t res, const acb_poly_t s, const acb_t z, slong n, slong prec)
312 {
313     if (n == 0)
314     {
315         acb_poly_zero(res);
316         return;
317     }
318 
319     acb_poly_fit_length(res, n);
320 
321     if (s->length == 0)
322     {
323         acb_t t;
324         acb_init(t);
325         _acb_poly_polylog_series(res->coeffs, t, 1, z, n, prec);
326         acb_clear(t);
327     }
328     else
329     {
330         _acb_poly_polylog_series(res->coeffs, s->coeffs, s->length, z, n, prec);
331     }
332 
333     _acb_poly_set_length(res, n);
334     _acb_poly_normalise(res);
335 }
336