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_modular.h"
13
14 static const int pentagonal_best_m[] = {
15 2, 5, 7, 11, 13, 17, 19, 23, 25, 35,
16 55, 65, 77, 91, 119, 133, 143, 175, 275, 325,
17 385, 455, 595, 665, 715, 935, 1001, 1309, 1463, 1547,
18 1729, 1925, 2275, 2975, 3325, 3575, 4675, 5005, 6545, 7315,
19 7735, 8645, 10465, 11305, 12155, 13585, 16445, 17017, 19019, 23023,
20 24871, 25025, 32725, 36575, 38675, 43225, 52325, 56525, 60775, 67925,
21 82225, 85085, 95095, 115115, 124355, 145145, 146965, 168245, 177905, 198835,
22 224315, 230945, 279565, 312455, 323323, 391391, 425425, 475475, 575575, 621775,
23 725725, 734825, 841225, 889525, 994175, 0
24 };
25
26 static const int pentagonal_best_m_residues[] = {
27 2, 3, 4, 6, 7, 9, 10, 12, 11, 12,
28 18, 21, 24, 28, 36, 40, 42, 44, 66, 77,
29 72, 84, 108, 120, 126, 162, 168, 216, 240, 252,
30 280, 264, 308, 396, 440, 462, 594, 504, 648, 720,
31 756, 840, 1008, 1080, 1134, 1260, 1512, 1512, 1680, 2016,
32 2160, 1848, 2376, 2640, 2772, 3080, 3696, 3960, 4158, 4620,
33 5544, 4536, 5040, 6048, 6480, 7560, 7560, 8640, 9072, 10080,
34 11340, 11340, 13608, 15120, 15120, 18144, 16632, 18480, 22176, 23760,
35 27720, 27720, 31680, 33264, 36960, 0
36 };
37
38 slong acb_modular_rs_optimal_m(const int * best_ms, const int * num_residues, slong N);
39
40 #define PENTAGONAL(N) ((((N)+2)/2) * ((3*(N)+5)/2)/2)
41
42 void
_acb_modular_mul(acb_t z,acb_t tmp1,acb_t tmp2,const acb_t x,const acb_t y,slong wprec,slong prec)43 _acb_modular_mul(acb_t z, acb_t tmp1, acb_t tmp2, const acb_t x, const acb_t y, slong wprec, slong prec)
44 {
45 if (prec <= 1024)
46 {
47 acb_mul(z, x, y, wprec);
48 }
49 else if (x == y)
50 {
51 acb_set_round(tmp1, x, wprec);
52 acb_mul(z, tmp1, tmp1, wprec);
53 }
54 else
55 {
56 acb_set_round(tmp1, x, wprec);
57 acb_set_round(tmp2, y, wprec);
58 acb_mul(z, tmp1, tmp2, wprec);
59 }
60 }
61
62 void
_acb_modular_eta_sum_basecase(acb_t eta,const acb_t q,double log2q_approx,slong N,slong prec)63 _acb_modular_eta_sum_basecase(acb_t eta, const acb_t q, double log2q_approx, slong N, slong prec)
64 {
65 slong e, e1, e2, k, k1, k2, num, term_prec;
66 slong *exponents, *aindex, *bindex;
67 acb_ptr qpow;
68 acb_t tmp1, tmp2;
69 double log2term_approx;
70
71 if (N <= 5)
72 {
73 if (N <= 1)
74 {
75 acb_set_ui(eta, N != 0);
76 }
77 else if (N == 2)
78 {
79 acb_sub_ui(eta, q, 1, prec);
80 acb_neg(eta, eta);
81 }
82 else
83 {
84 acb_mul(eta, q, q, prec);
85 acb_add(eta, eta, q, prec);
86 acb_neg(eta, eta);
87 acb_add_ui(eta, eta, 1, prec);
88 }
89 return;
90 }
91
92 num = 1;
93 while (PENTAGONAL(num) < N)
94 num++;
95
96 acb_init(tmp1);
97 acb_init(tmp2);
98
99 exponents = flint_malloc(sizeof(slong) * 3 * num);
100 aindex = exponents + num;
101 bindex = aindex + num;
102
103 qpow = _acb_vec_init(num);
104
105 acb_modular_addseq_eta(exponents, aindex, bindex, num);
106 acb_set_round(qpow + 0, q, prec);
107
108 acb_zero(eta);
109
110 for (k = 0; k < num; k++)
111 {
112 e = exponents[k];
113
114 log2term_approx = e * log2q_approx;
115 term_prec = FLINT_MIN(FLINT_MAX(prec + log2term_approx + 16.0, 16.0), prec);
116
117 if (k > 0)
118 {
119 k1 = aindex[k];
120 k2 = bindex[k];
121 e1 = exponents[k1];
122 e2 = exponents[k2];
123
124 if (e == e1 + e2)
125 {
126 _acb_modular_mul(qpow + k, tmp1, tmp2, qpow + k1, qpow + k2, term_prec, prec);
127 }
128 else if (e == 2 * e1 + e2)
129 {
130 _acb_modular_mul(qpow + k, tmp1, tmp2, qpow + k1, qpow + k1, term_prec, prec);
131 _acb_modular_mul(qpow + k, tmp1, tmp2, qpow + k, qpow + k2, term_prec, prec);
132 }
133 else
134 {
135 flint_printf("exponent not in addition sequence!\n");
136 flint_abort();
137 }
138 }
139
140 if (k % 4 <= 1)
141 acb_sub(eta, eta, qpow + k, prec);
142 else
143 acb_add(eta, eta, qpow + k, prec);
144 }
145
146 acb_add_ui(eta, eta, 1, prec);
147
148 flint_free(exponents);
149 _acb_vec_clear(qpow, num);
150 acb_clear(tmp1);
151 acb_clear(tmp2);
152 }
153
154 void
_acb_modular_eta_sum_rs(acb_t eta,const acb_t q,double log2q_approx,slong N,slong prec)155 _acb_modular_eta_sum_rs(acb_t eta, const acb_t q, double log2q_approx, slong N, slong prec)
156 {
157 slong * tab;
158 slong k, term_prec, i, e, eprev;
159 slong m, num_pentagonal;
160 double log2term_approx;
161 acb_ptr qpow;
162 acb_t tmp1, tmp2;
163
164 acb_init(tmp1);
165 acb_init(tmp2);
166
167 /* choose rectangular splitting parameters */
168 m = acb_modular_rs_optimal_m(pentagonal_best_m, pentagonal_best_m_residues, N);
169
170 /* build addition sequence */
171 tab = flint_calloc(m + 1, sizeof(slong));
172
173 for (k = 0; PENTAGONAL(k) < N; k++)
174 tab[PENTAGONAL(k) % m] = -1;
175 num_pentagonal = k;
176 tab[m] = -1;
177
178 /* compute powers in addition sequence */
179 qpow = _acb_vec_init(m + 1);
180 acb_modular_fill_addseq(tab, m + 1);
181
182 for (k = 0; k < m + 1; k++)
183 {
184 if (k == 0)
185 {
186 acb_one(qpow + k);
187 }
188 else if (k == 1)
189 {
190 acb_set_round(qpow + k, q, prec);
191 }
192 else if (tab[k] != 0)
193 {
194 log2term_approx = k * log2q_approx;
195 term_prec = FLINT_MIN(FLINT_MAX(prec + log2term_approx + 16.0, 16.0), prec);
196 _acb_modular_mul(qpow + k, tmp1, tmp2, qpow + tab[k], qpow + k - tab[k], term_prec, prec);
197 }
198 }
199
200 /* compute eta */
201 acb_zero(eta);
202 term_prec = prec;
203
204 for (k = num_pentagonal - 1; k >= 0; k--)
205 {
206 e = PENTAGONAL(k); /* exponent */
207 eprev = PENTAGONAL(k+1);
208
209 log2term_approx = e * log2q_approx;
210 term_prec = FLINT_MIN(FLINT_MAX(prec + log2term_approx + 16.0, 16.0), prec);
211
212 /* giant steps */
213 for (i = e / m; i < eprev / m; i++)
214 {
215 if (!acb_is_zero(eta))
216 _acb_modular_mul(eta, tmp1, tmp2, eta, qpow + m, term_prec, prec);
217 }
218
219 if (k % 4 <= 1)
220 acb_sub(eta, eta, qpow + (e % m), prec);
221 else
222 acb_add(eta, eta, qpow + (e % m), prec);
223 }
224
225 acb_add_ui(eta, eta, 1, prec);
226
227 acb_clear(tmp1);
228 acb_clear(tmp2);
229
230 _acb_vec_clear(qpow, m + 1);
231 flint_free(tab);
232 }
233
234 void
acb_modular_eta_sum(acb_t eta,const acb_t q,slong prec)235 acb_modular_eta_sum(acb_t eta, const acb_t q, slong prec)
236 {
237 mag_t err, qmag;
238 double log2q_approx;
239 int q_is_real;
240 slong N;
241
242 mag_init(err);
243 mag_init(qmag);
244
245 q_is_real = arb_is_zero(acb_imagref(q));
246
247 acb_get_mag(qmag, q);
248 log2q_approx = mag_get_d_log2_approx(qmag);
249
250 if (log2q_approx >= 0.0)
251 {
252 N = 1;
253 mag_inf(err);
254 }
255 else /* Pick N and compute error bound */
256 {
257 N = 0;
258 while (0.05 * N * N < prec)
259 {
260 if (log2q_approx * PENTAGONAL(N) < -prec - 2)
261 break;
262 N++;
263 }
264 N = PENTAGONAL(N);
265
266 mag_geom_series(err, qmag, N);
267 if (mag_is_inf(err))
268 N = 1;
269 }
270
271 if (N < 400)
272 _acb_modular_eta_sum_basecase(eta, q, log2q_approx, N, prec);
273 else
274 _acb_modular_eta_sum_rs(eta, q, log2q_approx, N, prec);
275
276 if (q_is_real)
277 arb_add_error_mag(acb_realref(eta), err);
278 else
279 acb_add_error_mag(eta, err);
280
281 mag_clear(err);
282 mag_clear(qmag);
283 }
284
285