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