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
_acb_poly_acb_invpow_cpx(acb_ptr res,const acb_t N,const acb_t c,slong trunc,slong prec)14 void _acb_poly_acb_invpow_cpx(acb_ptr res, const acb_t N, const acb_t c, slong trunc, slong prec)
15 {
16 slong i;
17 acb_t logN;
18
19 acb_init(logN);
20 acb_log(logN, N, prec);
21 acb_mul(res + 0, logN, c, prec);
22 acb_neg(res + 0, res + 0);
23 acb_exp(res + 0, res + 0, prec);
24
25 for (i = 1; i < trunc; i++)
26 {
27 acb_mul(res + i, res + i - 1, logN, prec);
28 acb_div_si(res + i, res + i, -i, prec);
29 }
30
31 acb_clear(logN);
32 }
33
34 void
_acb_poly_powsum_series_naive(acb_ptr z,const acb_t s,const acb_t a,const acb_t q,slong n,slong len,slong prec)35 _acb_poly_powsum_series_naive(acb_ptr z,
36 const acb_t s, const acb_t a, const acb_t q, slong n, slong len, slong prec)
37 {
38 slong k, i;
39 int q_one, s_int;
40 acb_t ak, logak, t, qpow, negs;
41
42 acb_init(ak);
43 acb_init(logak);
44 acb_init(t);
45 acb_init(qpow);
46 acb_init(negs);
47
48 _acb_vec_zero(z, len);
49 acb_one(qpow);
50 acb_neg(negs, s);
51
52 q_one = acb_is_one(q);
53 s_int = arb_is_int(acb_realref(s)) && arb_is_zero(acb_imagref(s));
54
55 for (k = 0; k < n; k++)
56 {
57 acb_add_ui(ak, a, k, prec);
58
59 if (len == 1)
60 {
61 acb_pow(t, ak, negs, prec);
62 }
63 else
64 {
65 acb_log(logak, ak, prec);
66
67 if (s_int)
68 {
69 acb_pow(t, ak, negs, prec);
70 }
71 else
72 {
73 acb_mul(t, logak, negs, prec);
74 acb_exp(t, t, prec);
75 }
76 }
77
78 if (!q_one)
79 {
80 acb_mul(t, t, qpow, prec);
81 if (k < n - 1)
82 acb_mul(qpow, qpow, q, prec);
83 }
84
85 acb_add(z, z, t, prec);
86
87 for (i = 1; i < len; i++)
88 {
89 acb_mul(t, t, logak, prec);
90 acb_div_si(t, t, -i, prec);
91 acb_add(z + i, z + i, t, prec);
92 }
93 }
94
95 acb_clear(ak);
96 acb_clear(logak);
97 acb_clear(t);
98 acb_clear(qpow);
99 acb_clear(negs);
100 }
101
102