1 /*
2     Copyright (C) 2019 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_dirichlet.h"
13 
14 /*
15 Claim: the error is bounded by 1/64 if n <= 1 and (1/64) (log(n)/n) if n >= 2.
16 
17 A crude lower bound for g_n is 2 pi exp(W(n)), or 8*n/log(n) for n >= 8.
18 
19 We want to solve pi n = -t/2 log(2 pi/t) - t/2 - pi/8 + epsilon for t (= g_n).
20 Using (47) in Brent [https://arxiv.org/abs/1609.03682], |epsilon| <= 1/(8 t) for for t >= 2.
21 Also, for x >= 3, |f'(x)| < 0.5 where f(x) = exp(W(x)).
22 
23 Assume n >= 9, so that (n+1/8)/e >= 3.35. Then inverting gives
24 
25 t = 2 pi exp[W( [pi n - epsilon + pi/8] / (pi e) ) + 1]
26   = 2 pi e exp[W((n+1/8)/e  - epsilon / (pi e))]
27   = 2 pi e exp[W((n+1/8)/e)] + epsilon2, |epsilon2| <= 1/(8 t) <= (1/64) (log(n)/n)
28 
29 One can check 0 <= n <= 8 separately.
30 */
31 static void
gram_point_initial(arb_t x,const fmpz_t n,slong prec)32 gram_point_initial(arb_t x, const fmpz_t n, slong prec)
33 {
34     arb_t pi, e;
35     mag_t b;
36 
37     arb_init(pi);
38     arb_init(e);
39     mag_init(b);
40 
41     arb_const_pi(pi, prec);
42     arb_const_e(e, prec);
43 
44     /* x = 2*pi*exp(1 + W((n+1/8)/e)) */
45     arb_one(x);
46     arb_mul_2exp_si(x, x, -3);
47     arb_add_fmpz(x, x, n, prec);
48     arb_div(x, x, e, prec);
49     arb_lambertw(x, x, 0, prec);
50     arb_add_ui(x, x, 1, prec);
51     arb_exp(x, x, prec);
52     arb_mul(x, x, pi, prec);
53     arb_mul_2exp_si(x, x, 1);
54 
55     if (fmpz_cmp_ui(n, 1) <= 0)
56     {
57         mag_set_ui_2exp_si(b, 1, -6);
58     }
59     else
60     {
61         mag_set_fmpz(b, n);
62         mag_log(b, b);
63         mag_div_fmpz(b, b, n);
64         mag_mul_2exp_si(b, b, -6);
65     }
66 
67     arb_add_error_mag(x, b);
68 
69     arb_clear(pi);
70     arb_clear(e);
71     mag_clear(b);
72 }
73 
74 void
acb_dirichlet_gram_point(arb_t res,const fmpz_t n,const dirichlet_group_t G,const dirichlet_char_t chi,slong prec)75 acb_dirichlet_gram_point(arb_t res, const fmpz_t n, const dirichlet_group_t G, const dirichlet_char_t chi, slong prec)
76 {
77     slong asymp_accuracy;
78 
79     /* Only implemented for n >= -1 and Riemann zeta. */
80     if (fmpz_cmp_si(n, -1) < 0 || G != NULL || chi != NULL)
81     {
82         arb_indeterminate(res);
83         return;
84     }
85 
86     asymp_accuracy = 2 * fmpz_bits(n);
87     asymp_accuracy = FLINT_MIN(asymp_accuracy, prec);
88 
89     gram_point_initial(res, n, asymp_accuracy + 20);
90     asymp_accuracy = arb_rel_accuracy_bits(res);
91 
92     if (asymp_accuracy < prec)
93     {
94         acb_struct tmp[2];
95         arb_t f, fprime, root;
96         mag_t C, r;
97         slong * steps;
98         slong wp, step;
99 
100         acb_init(tmp);
101         acb_init(tmp + 1);
102         arb_init(f);
103         arb_init(fprime);
104         arb_init(root);
105         mag_init(C);
106         mag_init(r);
107         steps = flint_malloc(sizeof(slong) * FLINT_BITS);
108 
109         step = 0;
110         steps[step] = prec * 1.05 + 10;
111 
112         while (steps[step] / 2 > asymp_accuracy)
113         {
114             steps[step + 1] = steps[step] / 2;
115             step++;
116         }
117 
118         arb_set(root, res);
119 
120         /* theta''(x) <= C = 1/x, x >= 1 */
121         arb_get_mag_lower(C, root);
122         if (mag_cmp_2exp_si(C, 0) >= 0)
123             mag_inv(C, C);
124         else
125             mag_inf(C);
126 
127         arb_set(root, res);
128 
129         for ( ; step >= 0; step--)
130         {
131             wp = steps[step] + 10;
132             wp = FLINT_MAX(wp, arb_rel_accuracy_bits(root) + 10);
133 
134             /* store radius, set root to the midpoint */
135             mag_set(r, arb_radref(root));
136             mag_zero(arb_radref(root));
137 
138             acb_set_arb(tmp, root);
139             acb_dirichlet_hardy_theta(tmp, tmp, NULL, NULL, 2, wp);
140             arb_set(f, acb_realref(tmp));
141             arb_const_pi(acb_imagref(tmp), wp);
142             arb_submul_fmpz(f, acb_imagref(tmp), n, wp);
143 
144             arb_set(fprime, acb_realref(tmp + 1));
145 
146             /* f'([m+/-r]) = f'(m) +/- f''([m +/- r]) * r */
147             mag_mul(r, C, r);
148             arb_add_error_mag(fprime, r);
149             arb_div(f, f, fprime, wp);
150             arb_sub(root, root, f, wp);
151 
152             /* Verify inclusion so that C is still valid. */
153             if (!arb_contains(res, root))
154             {
155                 flint_printf("unexpected: no containment computing Gram point\n");
156                 arb_set(root, res);
157                 break;
158             }
159         }
160 
161         arb_set(res, root);
162 
163         acb_clear(tmp);
164         acb_clear(tmp + 1);
165         arb_clear(f);
166         arb_clear(fprime);
167         arb_clear(root);
168         mag_clear(C);
169         mag_clear(r);
170         flint_free(steps);
171     }
172 
173     arb_set_round(res, res, prec);
174 }
175