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