1 /*
2 Copyright (C) 2017 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_lambertw_series(acb_ptr res,acb_srcptr z,slong zlen,const fmpz_t k,int flags,slong len,slong prec)14 void _acb_poly_lambertw_series(acb_ptr res,
15 acb_srcptr z, slong zlen, const fmpz_t k, int flags, slong len, slong prec)
16 {
17 acb_ptr w, ew, t, u;
18 acb_t ew0;
19
20 zlen = FLINT_MIN(zlen, len);
21
22 if (zlen == 1)
23 {
24 acb_lambertw(res, z, k, flags, prec);
25 _acb_vec_zero(res + 1, len - 1);
26 return;
27 }
28
29 w = _acb_vec_init(len);
30 ew = _acb_vec_init(len);
31 t = _acb_vec_init(len);
32 u = _acb_vec_init(len);
33 acb_init(ew0);
34
35 acb_lambertw(w, z, k, flags, prec);
36
37 if (acb_contains_zero(w))
38 acb_exp(ew0, w, prec);
39 else
40 acb_div(ew0, z, w, prec);
41
42 acb_add(t, ew0, z, prec);
43 acb_div(w + 1, z + 1, t, prec);
44
45 NEWTON_INIT(2, len)
46 NEWTON_LOOP(m, n)
47
48 /* _acb_poly_exp_series(ew, w, m, n, prec); */
49 acb_zero(t);
50 _acb_vec_set(t + 1, w + 1, m - 1);
51 _acb_poly_exp_series(ew, t, m, n, prec);
52 _acb_vec_scalar_mul(ew, ew, n, ew0, prec);
53
54 _acb_poly_mullow(t, ew, n, w, m, n, prec);
55 _acb_poly_sub(u, t, n, z, FLINT_MIN(zlen, n), prec);
56 _acb_vec_add(t, t, ew, n, prec);
57 _acb_poly_div_series(ew, u, n, t, n, n, prec);
58 _acb_vec_neg(w + m, ew + m, n - m);
59
60 NEWTON_END_LOOP
61 NEWTON_END
62
63 _acb_vec_set(res, w, len);
64
65 _acb_vec_clear(w, len);
66 _acb_vec_clear(ew, len);
67 _acb_vec_clear(t, len);
68 _acb_vec_clear(u, len);
69 acb_clear(ew0);
70 }
71
72 void
acb_poly_lambertw_series(acb_poly_t res,const acb_poly_t z,const fmpz_t k,int flags,slong len,slong prec)73 acb_poly_lambertw_series(acb_poly_t res,
74 const acb_poly_t z, const fmpz_t k, int flags, slong len, slong prec)
75 {
76 if (len == 0 || (fmpz_is_zero(k) && z->length == 0))
77 {
78 acb_poly_zero(res);
79 return;
80 }
81
82 if (z->length == 0)
83 {
84 acb_poly_fit_length(res, len);
85 _acb_vec_indeterminate(res->coeffs, len);
86 _acb_poly_set_length(res, len);
87 return;
88 }
89
90 acb_poly_fit_length(res, len);
91 _acb_poly_lambertw_series(res->coeffs, z->coeffs, z->length, k, flags, len, prec);
92 _acb_poly_set_length(res, len);
93 _acb_poly_normalise(res);
94 }
95
96