1 /*
2 Copyright (C) 2015 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 "arb.h"
13 #include "bernoulli.h"
14 #include "arb_poly.h"
15
16 /* todo: don't use exact bernoulli numbers for large len */
17 /* todo: output exact integers when precise enough */
18 void
arb_power_sum_vec(arb_ptr res,const arb_t a,const arb_t b,slong len,slong prec)19 arb_power_sum_vec(arb_ptr res, const arb_t a, const arb_t b, slong len, slong prec)
20 {
21 arb_ptr t, u, v;
22 slong k;
23
24 if (len < 1)
25 return;
26
27 t = _arb_vec_init(len + 1);
28 u = _arb_vec_init(len + 1);
29 v = _arb_vec_init(len + 1);
30
31 /* exp(ax), exp(bx) */
32 arb_set(t + 1, a);
33 arb_set(u + 1, b);
34 _arb_poly_exp_series(t, t, 2, len + 1, prec);
35 _arb_poly_exp_series(u, u, 2, len + 1, prec);
36 _arb_vec_sub(t, u, t, len + 1, prec);
37
38 /* x/(exp(x)-1) */
39 BERNOULLI_ENSURE_CACHED(len + 1);
40 for (k = 0; k <= len; k++)
41 arb_set_fmpq(u + k, bernoulli_cache + k, prec);
42 _arb_poly_borel_transform(u, u, len + 1, prec);
43
44 _arb_poly_mullow(v, t, len + 1, u, len + 1, len + 1, prec);
45 _arb_poly_inv_borel_transform(v, v, len + 1, prec);
46
47 for (k = 0; k < len; k++)
48 arb_div_ui(res + k, v + k + 1, k + 1, prec);
49
50 _arb_vec_clear(t, len + 1);
51 _arb_vec_clear(u, len + 1);
52 _arb_vec_clear(v, len + 1);
53 }
54
55