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.h"
13
14 void
acb_lambertw_asymp(acb_t res,const acb_t z,const fmpz_t k,slong L,slong M,slong prec)15 acb_lambertw_asymp(acb_t res, const acb_t z, const fmpz_t k, slong L, slong M, slong prec)
16 {
17 acb_t L1, L2, sigma, tau, s, c, u;
18 slong l, m;
19 fmpz_t t;
20 fmpz * sc;
21
22 /* For k = 0, the asymptotic expansion is not valid near 0. */
23 /* (It is sufficient to look at the midpoint as a test here.) */
24 if (fmpz_is_zero(k) && arf_cmpabs_2exp_si(arb_midref(acb_realref(z)), 0) < 0
25 && arf_cmpabs_2exp_si(arb_midref(acb_imagref(z)), 0) < 0)
26 {
27 acb_indeterminate(res);
28 return;
29 }
30
31 acb_init(L1);
32 acb_init(L2);
33 acb_init(sigma);
34 acb_init(tau);
35 acb_init(s);
36 acb_init(c);
37 acb_init(u);
38 fmpz_init(t);
39
40 acb_const_pi(L2, prec);
41 acb_mul_2exp_si(L2, L2, 1);
42 acb_mul_fmpz(L2, L2, k, prec);
43 acb_mul_onei(L2, L2);
44 acb_log(L1, z, prec);
45 acb_add(L1, L1, L2, prec);
46 acb_log(L2, L1, prec);
47
48 acb_inv(sigma, L1, prec);
49 acb_mul(tau, L2, sigma, prec);
50
51 acb_zero(s);
52
53 /* Stirling numbers */
54 sc = _fmpz_vec_init(L);
55
56 acb_one(u);
57
58 for (m = 1; m < M; m++)
59 {
60 if (m == 1)
61 {
62 for (l = 0; l < L; l++)
63 fmpz_one(sc + l);
64 }
65 else
66 {
67 for (l = 0; l < L; l++)
68 {
69 fmpz_mul_ui(sc + l, sc + l, m + l - 1);
70 if (l > 0)
71 fmpz_add(sc + l, sc + l, sc + l - 1);
72 }
73 }
74
75 acb_zero(c);
76
77 /* todo: precompute powers instead of horner */
78 for (l = L - 1; l >= 0; l--)
79 {
80 acb_mul(c, c, sigma, prec);
81 if (l % 2)
82 acb_sub_fmpz(c, c, sc + l, prec);
83 else
84 acb_add_fmpz(c, c, sc + l, prec);
85 }
86
87 acb_mul(u, u, tau, prec);
88 acb_div_ui(u, u, m, prec);
89 acb_addmul(s, c, u, prec);
90 }
91
92 _fmpz_vec_clear(sc, L);
93
94 acb_sub(s, s, L2, prec);
95 acb_add(s, s, L1, prec);
96
97 {
98 mag_t m4s, m4t, one, q, r;
99
100 mag_init(m4s);
101 mag_init(m4t);
102 mag_init(one);
103 mag_init(q);
104 mag_init(r);
105
106 acb_get_mag(m4s, sigma);
107 mag_mul_2exp_si(m4s, m4s, 2);
108 acb_get_mag(m4t, tau);
109 mag_mul_2exp_si(m4t, m4t, 2);
110
111 mag_one(one);
112
113 mag_sub_lower(q, one, m4s);
114 mag_sub_lower(r, one, m4t);
115 mag_mul(q, q, r);
116
117 mag_pow_ui(r, m4s, L);
118 mag_mul(r, r, m4t);
119 mag_pow_ui(m4t, m4t, M);
120 mag_add(r, r, m4t);
121
122 mag_div(q, r, q);
123
124 acb_add_error_mag(s, q);
125
126 mag_clear(m4s);
127 mag_clear(m4t);
128 mag_clear(one);
129 mag_clear(q);
130 mag_clear(r);
131 }
132
133 acb_set(res, s);
134
135 acb_clear(sigma);
136 acb_clear(tau);
137 acb_clear(s);
138 acb_clear(c);
139 acb_clear(L1);
140 acb_clear(L2);
141 acb_clear(u);
142 fmpz_clear(t);
143 }
144
145