1 /*
2 Copyright (C) 2015 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 square_best_m[] = {
15 2, 3, 4, 8, 12, 16, 32, 48, 80, 96,
16 112, 144, 240, 288, 336, 480, 560, 576, 720, 1008,
17 1440, 1680, 2016, 2640, 2880, 3600, 4032, 5040, 7920, 9360,
18 10080, 15840, 18480, 20160, 25200, 31680, 37440, 39600, 44352, 50400,
19 55440, 65520, 85680, 95760, 102960, 110880, 131040, 171360, 191520, 205920,
20 221760, 262080, 277200, 327600, 383040, 411840, 514800, 554400, 655200, 720720,
21 942480, 0
22 };
23
24 static const int square_best_m_residues[] = {
25 2, 2, 2, 3, 4, 4, 7, 8, 12, 14,
26 16, 16, 24, 28, 32, 42, 48, 48, 48, 64,
27 84, 96, 112, 144, 144, 176, 192, 192, 288, 336,
28 336, 504, 576, 576, 704, 864, 1008, 1056, 1152, 1232,
29 1152, 1344, 1728, 1920, 2016, 2016, 2352, 3024, 3360, 3528,
30 3456, 4032, 4224, 4928, 5760, 6048, 7392, 7392, 8624, 8064,
31 10368, 0
32 };
33
34 static const int trigonal_best_m[] = {
35 2, 6, 10, 14, 18, 30, 42, 66, 70, 90,
36 126, 198, 210, 330, 390, 450, 630, 990, 1170, 1386,
37 1638, 2142, 2310, 2730, 3150, 4950, 5850, 6930, 8190, 10710,
38 11970, 12870, 16830, 18018, 23562, 26334, 27846, 30030, 34650, 40950,
39 53550, 59850, 64350, 84150, 90090, 117810, 131670, 139230, 155610, 188370,
40 203490, 218790, 244530, 270270, 296010, 306306, 342342, 414414, 447678, 450450,
41 589050, 658350, 696150, 778050, 941850, 0
42 };
43
44 static const int trigonal_best_m_residues[] = {
45 1, 2, 3, 4, 4, 6, 8, 12, 12, 12,
46 16, 24, 24, 36, 42, 44, 48, 72, 84, 96,
47 112, 144, 144, 168, 176, 264, 308, 288, 336, 432,
48 480, 504, 648, 672, 864, 960, 1008, 1008, 1056, 1232,
49 1584, 1760, 1848, 2376, 2016, 2592, 2880, 3024, 3360, 4032,
50 4320, 4536, 5040, 5544, 6048, 6048, 6720, 8064, 8640, 7392,
51 9504, 10560, 11088, 12320, 14784, 0,
52 };
53
54 slong
acb_modular_rs_optimal_m(const int * best_ms,const int * num_residues,slong N)55 acb_modular_rs_optimal_m(const int * best_ms, const int * num_residues, slong N)
56 {
57 slong i, m, cost, best_i, best_m, best_cost;
58
59 best_i = 0;
60 best_m = best_ms[0];
61 best_cost = WORD_MAX;
62
63 for (i = 0; (m = best_ms[i]) != 0; i++)
64 {
65 cost = N / m + num_residues[i];
66
67 if (i == 0 || cost < best_cost)
68 {
69 best_i = i;
70 best_cost = cost;
71 best_m = m;
72 }
73 }
74
75 /* flint_printf("N = %wd, best_m = %wd, best_cost = %wd, s(m) = %d\n",
76 N, best_m, best_cost, num_residues[best_i]); */
77 i = best_i;
78
79 return best_m;
80 }
81
82 void
acb_modular_theta_const_sum_rs(acb_t theta2,acb_t theta3,acb_t theta4,const acb_t q,slong N,slong prec)83 acb_modular_theta_const_sum_rs(acb_t theta2, acb_t theta3, acb_t theta4,
84 const acb_t q, slong N, slong prec)
85 {
86 slong * tab;
87 slong k, term_prec, i, e, eprev;
88 slong M, m2, m3, num_square, num_trigonal;
89 double log2q_approx, log2term_approx;
90 acb_ptr qpow;
91 acb_t tmp1, tmp2;
92 mag_t qmag;
93
94 mag_init(qmag);
95 acb_get_mag(qmag, q);
96 log2q_approx = mag_get_d_log2_approx(qmag);
97 mag_clear(qmag);
98
99 acb_init(tmp1);
100 acb_init(tmp2);
101
102 /* choose rectangular splitting parameters */
103 m2 = acb_modular_rs_optimal_m(trigonal_best_m, trigonal_best_m_residues, N);
104 m3 = acb_modular_rs_optimal_m(square_best_m, square_best_m_residues, N);
105 M = FLINT_MAX(m2, m3) + 1;
106
107 /* build addition sequence */
108 tab = flint_calloc(M, sizeof(slong));
109
110 for (k = 0; k*(k+1) < N; k++)
111 tab[(k*(k+1)) % m2] = -1;
112 num_trigonal = k;
113
114 for (k = 0; k*k < N; k++)
115 tab[(k*k) % m3] = -1;
116 num_square = k;
117
118 tab[m2] = -1;
119 tab[m3] = -1;
120
121 /* compute powers in addition sequence */
122 qpow = _acb_vec_init(M);
123 acb_modular_fill_addseq(tab, M);
124
125 for (k = 0; k < M; k++)
126 {
127 if (k == 0)
128 {
129 acb_one(qpow + k);
130 }
131 else if (k == 1)
132 {
133 acb_set_round(qpow + k, q, prec);
134 }
135 else if (tab[k] != 0)
136 {
137 log2term_approx = k * log2q_approx;
138 term_prec = FLINT_MIN(FLINT_MAX(prec + log2term_approx + 16.0, 16.0), prec);
139 _acb_modular_mul(qpow + k, tmp1, tmp2, qpow + tab[k], qpow + k - tab[k], term_prec, prec);
140 }
141 }
142
143 /* compute theta2 */
144 acb_zero(theta2);
145 term_prec = prec;
146
147 for (k = num_trigonal - 1; k >= 0; k--)
148 {
149 e = k * (k + 1); /* exponent */
150 eprev = (k + 1) * (k + 2);
151
152 log2term_approx = e * log2q_approx;
153 term_prec = FLINT_MIN(FLINT_MAX(prec + log2term_approx + 16.0, 16.0), prec);
154
155 /* giant steps */
156 for (i = e / m2; i < eprev / m2; i++)
157 {
158 if (!acb_is_zero(theta2))
159 _acb_modular_mul(theta2, tmp1, tmp2, theta2, qpow + m2, term_prec, prec);
160 }
161
162 acb_add(theta2, theta2, qpow + (e % m2), prec);
163 }
164
165 acb_mul_2exp_si(theta2, theta2, 1);
166
167 /* compute theta3, theta4 */
168 acb_zero(theta3);
169 acb_zero(theta4);
170 term_prec = prec;
171
172 for (k = num_square - 1; k >= 0; k--)
173 {
174 e = k * k; /* exponent */
175 eprev = (k + 1) * (k + 1);
176
177 log2term_approx = e * log2q_approx;
178 term_prec = FLINT_MIN(FLINT_MAX(prec + log2term_approx + 16.0, 16.0), prec);
179
180 /* giant steps */
181 for (i = e / m3; i < eprev / m3; i++)
182 {
183 if (!acb_is_zero(theta3))
184 _acb_modular_mul(theta3, tmp1, tmp2, theta3, qpow + m3, term_prec, prec);
185
186 if (!acb_is_zero(theta4))
187 _acb_modular_mul(theta4, tmp1, tmp2, theta4, qpow + m3, term_prec, prec);
188 }
189
190 if (k == 0)
191 {
192 acb_mul_2exp_si(theta3, theta3, 1);
193 acb_mul_2exp_si(theta4, theta4, 1);
194 }
195
196 acb_add(theta3, theta3, qpow + (e % m3), prec);
197
198 if (k % 2 == 0)
199 acb_add(theta4, theta4, qpow + (e % m3), prec);
200 else
201 acb_sub(theta4, theta4, qpow + (e % m3), prec);
202 }
203
204 acb_clear(tmp1);
205 acb_clear(tmp2);
206
207 _acb_vec_clear(qpow, M);
208 flint_free(tab);
209 }
210
211