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_bound_deriv(mag_t res,const acb_t z,const acb_t ez1,const fmpz_t k)15 acb_lambertw_bound_deriv(mag_t res, const acb_t z, const acb_t ez1, const fmpz_t k)
16 {
17     mag_t t, u, v;
18 
19     mag_init(t);
20     mag_init(u);
21     mag_init(v);
22 
23     if (fmpz_is_zero(k))
24     {
25         acb_get_mag(t, z);
26 
27         /* |z| <= 64 */
28         if (mag_cmp_2exp_si(t, 6) < 0)
29         {
30             /* 2.25 / sqrt(|t|) / sqrt(1+|t|), t = |ez+1| */
31             acb_get_mag_lower(t, ez1);
32             mag_one(u);
33             mag_add_lower(u, u, t);
34             mag_mul_lower(t, t, u);
35             mag_rsqrt(t, t);
36 
37             if (arb_is_positive(acb_realref(ez1)))
38             {
39                 mag_mul_ui(t, t, 135);      /* x0.9375 for small improvement */
40                 mag_mul_2exp_si(t, t, -6);
41             }
42             else
43             {
44                 mag_mul_ui(t, t, 9);
45                 mag_mul_2exp_si(t, t, -2);
46             }
47 
48             mag_set(res, t);
49         }
50         else
51         {
52             acb_get_mag_lower(t, z);
53 
54             if (mag_cmp_2exp_si(t, 2) >= 0)
55             {
56                 mag_one(u);
57                 mag_div(res, u, t);
58             }
59             else   /* unlikely */
60             {
61                 acb_get_mag_lower(u, ez1);
62                 mag_rsqrt(u, u);
63                 mag_mul_2exp_si(u, u, -1);
64                 mag_add_ui(u, u, 1);
65                 mag_mul_ui(u, u, 3);
66                 mag_div(res, u, t);
67             }
68         }
69     }
70     else if (fmpz_is_pm1(k))
71     {
72         if (arb_is_nonnegative(acb_realref(z)) ||
73             (fmpz_is_one(k) && arb_is_nonnegative(acb_imagref(z))) ||
74             (fmpz_equal_si(k, -1) && arb_is_negative(acb_imagref(z))))
75         {
76             /* (1 + 1/(4+|z|^2))/|z| */
77             acb_get_mag_lower(t, z);
78             mag_mul_lower(u, t, t);
79             mag_set_ui_lower(v, 4);
80             mag_add_lower(u, u, v);
81             mag_one(v);
82             mag_div(u, v, u);
83             mag_add(u, u, v);
84             mag_div(res, u, t);
85         }
86         else
87         {
88             /* (1 + 0.71875/sqrt(|e*z+1|)) / |z| */
89             acb_get_mag_lower(t, ez1);
90             mag_rsqrt(t, t);
91             mag_mul_ui(t, t, 23);
92             mag_mul_2exp_si(t, t, -5);
93             mag_one(u);
94             mag_add(t, t, u);
95             acb_get_mag_lower(u, z);
96             mag_div(res, t, u);
97         }
98 
99         mag_clear(t);
100         mag_clear(u);
101         mag_clear(v);
102         return;
103     }
104     else
105     {
106         /* |W'(z)|  = |1/z| |W(z)/(1+W(z))|
107                    <= |1/z| (2pi) / (2pi-1) */
108         mag_set_ui_2exp_si(t, 77, -6);
109         acb_get_mag_lower(res, z);
110         mag_div(res, t, res);
111     }
112 
113     mag_clear(t);
114     mag_clear(u);
115     mag_clear(v);
116 }
117 
118