1 /*
2 Copyright (C) 2014 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_poly.h"
13
14
15 static void
bound_I(arb_ptr I,const arb_t A,const arb_t B,const arb_t C,slong len,slong wp)16 bound_I(arb_ptr I, const arb_t A, const arb_t B, const arb_t C, slong len, slong wp)
17 {
18 slong k;
19
20 arb_t D, Dk, L, T, Bm1;
21
22 arb_init(D);
23 arb_init(Dk);
24 arb_init(Bm1);
25 arb_init(T);
26 arb_init(L);
27
28 arb_sub_ui(Bm1, B, 1, wp);
29 arb_one(L);
30
31 /* T = 1 / (A^Bm1 * Bm1) */
32 arb_inv(T, A, wp);
33 arb_pow(T, T, Bm1, wp);
34 arb_div(T, T, Bm1, wp);
35
36 if (len > 1)
37 {
38 arb_log(D, A, wp);
39 arb_add(D, D, C, wp);
40 arb_mul(D, D, Bm1, wp);
41 arb_set(Dk, D);
42 }
43
44 for (k = 0; k < len; k++)
45 {
46 if (k > 0)
47 {
48 arb_mul_ui(L, L, k, wp);
49 arb_add(L, L, Dk, wp);
50 arb_mul(Dk, Dk, D, wp);
51 }
52
53 arb_mul(I + k, L, T, wp);
54 arb_div(T, T, Bm1, wp);
55 }
56
57 arb_clear(D);
58 arb_clear(Dk);
59 arb_clear(Bm1);
60 arb_clear(T);
61 arb_clear(L);
62 }
63
64 /* 0.5*(B/AN)^2 + |B|/AN */
65 static void
bound_C(arb_t C,const arb_t AN,const arb_t B,slong wp)66 bound_C(arb_t C, const arb_t AN, const arb_t B, slong wp)
67 {
68 arb_t t;
69 arb_init(t);
70 arb_abs(t, B);
71 arb_div(t, t, AN, wp);
72 arb_mul_2exp_si(C, t, -1);
73 arb_add_ui(C, C, 1, wp);
74 arb_mul(C, C, t, wp);
75 arb_clear(t);
76 }
77
78 static void
bound_K(arb_t C,const arb_t AN,const arb_t B,const arb_t T,slong wp)79 bound_K(arb_t C, const arb_t AN, const arb_t B, const arb_t T, slong wp)
80 {
81 if (arb_is_zero(B) || arb_is_zero(T))
82 {
83 arb_one(C);
84 }
85 else
86 {
87 arb_div(C, B, AN, wp);
88 /* TODO: atan is dumb, should also bound by pi/2 */
89 arb_atan(C, C, wp);
90 arb_mul(C, C, T, wp);
91 if (arb_is_nonpositive(C))
92 arb_one(C);
93 else
94 arb_exp(C, C, wp);
95 }
96 }
97
98 static void
bound_rfac(arb_ptr F,const acb_t s,ulong n,slong len,slong wp)99 bound_rfac(arb_ptr F, const acb_t s, ulong n, slong len, slong wp)
100 {
101 if (len == 1)
102 {
103 acb_rising_ui_get_mag(arb_radref(F), s, n);
104 arf_set_mag(arb_midref(F), arb_radref(F));
105 mag_zero(arb_radref(F + 0));
106 }
107 else
108 {
109 arb_struct sx[2];
110 arb_init(sx + 0);
111 arb_init(sx + 1);
112 acb_abs(sx + 0, s, wp);
113 arb_one(sx + 1);
114 _arb_vec_zero(F, len);
115 _arb_poly_rising_ui_series(F, sx, 2, n, len, wp);
116 arb_clear(sx + 0);
117 arb_clear(sx + 1);
118 }
119 }
120
121 void
_acb_poly_zeta_em_bound(arb_ptr bound,const acb_t s,const acb_t a,ulong N,ulong M,slong len,slong wp)122 _acb_poly_zeta_em_bound(arb_ptr bound, const acb_t s, const acb_t a, ulong N, ulong M, slong len, slong wp)
123 {
124 arb_t K, C, AN, S2M;
125 arb_ptr F, R;
126 slong k;
127
128 arb_srcptr alpha = acb_realref(a);
129 arb_srcptr beta = acb_imagref(a);
130 arb_srcptr sigma = acb_realref(s);
131 arb_srcptr tau = acb_imagref(s);
132
133 arb_init(AN);
134 arb_init(S2M);
135
136 /* require alpha + N > 1, sigma + 2M > 1 */
137 arb_add_ui(AN, alpha, N - 1, wp);
138 arb_add_ui(S2M, sigma, 2*M - 1, wp);
139
140 if (!arb_is_positive(AN) || !arb_is_positive(S2M) || N < 1 || M < 1)
141 {
142 arb_clear(AN);
143 arb_clear(S2M);
144
145 for (k = 0; k < len; k++)
146 arb_pos_inf(bound + k);
147
148 return;
149 }
150
151 /* alpha + N, sigma + 2M */
152 arb_add_ui(AN, AN, 1, wp);
153 arb_add_ui(S2M, S2M, 1, wp);
154
155 R = _arb_vec_init(len);
156 F = _arb_vec_init(len);
157
158 arb_init(K);
159 arb_init(C);
160
161 /* bound for power integral */
162 bound_C(C, AN, beta, wp);
163 bound_K(K, AN, beta, tau, wp);
164 bound_I(R, AN, S2M, C, len, wp);
165
166 for (k = 0; k < len; k++)
167 {
168 arb_mul(R + k, R + k, K, wp);
169 arb_div_ui(K, K, k + 1, wp);
170 }
171
172 /* bound for rising factorial */
173 bound_rfac(F, s, 2*M, len, wp);
174
175 /* product (TODO: only need upper bound; write a function for this) */
176 _arb_poly_mullow(bound, F, len, R, len, len, wp);
177
178 /* bound for bernoulli polynomials, 4 / (2pi)^(2M) */
179 arb_const_pi(C, wp);
180 arb_mul_2exp_si(C, C, 1);
181 arb_pow_ui(C, C, 2 * M, wp);
182 arb_ui_div(C, 4, C, wp);
183 _arb_vec_scalar_mul(bound, bound, len, C, wp);
184
185 arb_clear(K);
186 arb_clear(C);
187 arb_clear(AN);
188 arb_clear(S2M);
189
190 _arb_vec_clear(R, len);
191 _arb_vec_clear(F, len);
192 }
193
194 void
_acb_poly_zeta_em_bound1(mag_t bound,const acb_t s,const acb_t a,slong N,slong M,slong len,slong wp)195 _acb_poly_zeta_em_bound1(mag_t bound,
196 const acb_t s, const acb_t a, slong N, slong M, slong len, slong wp)
197 {
198 arb_ptr vec = _arb_vec_init(len);
199 _acb_poly_zeta_em_bound(vec, s, a, N, M, len, wp);
200 _arb_vec_get_mag(bound, vec, len);
201 _arb_vec_clear(vec, len);
202 }
203
204