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