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 int
_acb_lambertw_check_branch(const acb_t w,const fmpz_t k,slong prec)15 _acb_lambertw_check_branch(const acb_t w, const fmpz_t k, slong prec)
16 {
17 arb_t t, u, v, a, b;
18 int res = 0;
19
20 arb_init(t);
21 arb_init(u);
22 arb_init(v);
23 arb_init(a);
24 arb_init(b);
25
26 /* t = x sinc(y), v = -cos(y) */
27 if (arb_is_exact(acb_imagref(w)))
28 {
29 if (arb_is_zero(acb_imagref(w)))
30 {
31 arb_one(t);
32 arb_one(v);
33 }
34 else
35 {
36 arb_sin_cos(t, v, acb_imagref(w), prec);
37 arb_div(t, t, acb_imagref(w), prec);
38 }
39 }
40 else
41 {
42 arb_sinc(t, acb_imagref(w), prec);
43 arb_cos(v, acb_imagref(w), prec);
44 }
45 arb_mul(t, t, acb_realref(w), prec);
46 arb_neg(v, v);
47
48 /* u = y / pi, with conjugate relation for k */
49 arb_const_pi(u, prec);
50 arb_div(u, acb_imagref(w), u, prec);
51 if (fmpz_sgn(k) < 0)
52 arb_neg(u, u);
53
54 if (fmpz_is_zero(k))
55 {
56 /* -1 < u < 1 and t > v */
57 arb_set_si(a, -1);
58 arb_set_si(b, 1);
59
60 if (arb_gt(u, a) && arb_lt(u, b) && arb_gt(t, v))
61 {
62 res = 1;
63 }
64 }
65 else
66 {
67 arb_set_fmpz(a, k);
68 arb_abs(a, a);
69 arb_mul_2exp_si(a, a, 1);
70 arb_add_ui(b, a, 1, prec);
71 arb_sub_ui(a, a, 2, prec);
72
73 /* if u > 2|k|-2 and u < 2|k|+1 */
74 if (arb_gt(u, a) && arb_lt(u, b))
75 {
76 arb_add_ui(a, a, 1, prec);
77 arb_sub_ui(b, b, 1, prec);
78
79 /* if u > 2|k|-1 and u < 2|k| */
80 if (arb_gt(u, a) && arb_lt(u, b))
81 {
82 res = 1;
83 }
84 else if (arb_lt(u, b) && arb_lt(t, v)) /* u < 2|k| and t < v */
85 {
86 res = 1;
87 }
88 else if (arb_gt(u, a) && arb_gt(t, v)) /* u > 2|k|-1 and t > v */
89 {
90 res = 1;
91 }
92 }
93 }
94
95 arb_clear(t);
96 arb_clear(u);
97 arb_clear(v);
98 arb_clear(a);
99 arb_clear(b);
100
101 return res;
102 }
103
104 int
acb_lambertw_check_branch(const acb_t w,const fmpz_t k,slong prec)105 acb_lambertw_check_branch(const acb_t w, const fmpz_t k, slong prec)
106 {
107 if (prec > 64 && _acb_lambertw_check_branch(w, k, 32))
108 return 1;
109
110 return _acb_lambertw_check_branch(w, k, prec);
111 }
112
113