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