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