1 /* Copyright (C) 2014 The PARI group.
2
3 This file is part of the PARI/GP package.
4
5 PARI/GP is free software; you can redistribute it and/or modify it under the
6 terms of the GNU General Public License as published by the Free Software
7 Foundation; either version 2 of the License, or (at your option) any later
8 version. It is distributed in the hope that it will be useful, but WITHOUT
9 ANY WARRANTY WHATSOEVER.
10
11 Check the License for details. You should have received a copy of it, along
12 with the package; see the file 'COPYING'. If not, write to the Free Software
13 Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
14
15 #include "pari.h"
16 #include "paripriv.h"
17
18 #define dbg_printf(lvl) if (DEBUGLEVEL >= (lvl) + 3) err_printf
19
20 /*
21 * SECTION: Functions dedicated to finding a j-invariant with a given
22 * trace.
23 */
24
25 /* TODO: This code is shared with
26 * torsion_compatible_with_characteristic() in 'torsion.c'. */
27 static void
hasse_bounds(long * low,long * high,long p)28 hasse_bounds(long *low, long *high, long p)
29 {
30 long two_sqrt_p = usqrt(4*p);
31 *low = p + 1 - two_sqrt_p;
32 *high = p + 1 + two_sqrt_p;
33 }
34
35 /* a / b : a and b are from factoru and b must divide a exactly */
36 INLINE GEN
famatsmall_divexact(GEN a,GEN b)37 famatsmall_divexact(GEN a, GEN b)
38 {
39 GEN a1 = gel(a,1), a2 = gel(a,2), c1, c2;
40 GEN b1 = gel(b,1), b2 = gel(b,2);
41 long i, j, k, la = lg(a1);
42 c1 = cgetg(la, t_VECSMALL);
43 c2 = cgetg(la, t_VECSMALL);
44 for (i = j = k = 1; j < la; j++)
45 {
46 c1[k] = a1[j];
47 c2[k] = a2[j];
48 if (a1[j] == b1[i]) { c2[k] -= b2[i++]; if (!c2[k]) continue; }
49 k++;
50 }
51 setlg(c1, k);
52 setlg(c2, k); return mkvec2(c1,c2);
53 }
54
55 /* This is Sutherland, 2009, TestCurveOrder.
56 *
57 * [a4, a6] and p specify an elliptic curve over FF_p. N0,N1 are the two
58 * possible curve orders, and n0,n1 their factoru */
59 static long
test_curve_order(norm_eqn_t ne,ulong a4,ulong a6,long N0,long N1,GEN n0,GEN n1,const long hasse[2])60 test_curve_order(norm_eqn_t ne, ulong a4, ulong a6,
61 long N0, long N1, GEN n0, GEN n1, const long hasse[2])
62 {
63 pari_sp ltop = avma, av;
64 ulong a4t, a6t, p = ne->p, pi = ne->pi, T = ne->T, swapped = 0;
65 long m0, m1, hasse_low, hasse_high;
66
67 if (p <= 11) {
68 long card = (long)p + 1 - Fl_elltrace(a4, a6, p);
69 return card == N0 || card == N1;
70 }
71 /* [a4, a6] is the given curve and [a4t, a6t] is its quadratic twist */
72 Fl_elltwist_disc(a4, a6, T, p, &a4t, &a6t);
73
74 m0 = m1 = 1;
75 if (N0 + N1 != 2 * (long)p + 2) pari_err_BUG("test_curve_order");
76
77 hasse_low = hasse[0];
78 hasse_high = hasse[1];
79 for (av = avma;;)
80 {
81 GEN pt, Q, fa0;
82 long a1, x, n_s;
83
84 pt = random_Flj_pre(a4, a6, p, pi);
85 Q = Flj_mulu_pre(pt, m0, a4, p, pi);
86 fa0 = m0 == 1? n0: famatsmall_divexact(n0, factoru(m0));
87 n_s = Flj_order_ufact(Q, N0 / m0, fa0, a4, p, pi);
88 if (n_s == 0) {
89 /* If m0 divides N1 and m1 divides N0 and N0 < N1, then swap */
90 if (!swapped && N1 % m0 == 0 && N0 % m1 == 0) {
91 swapspec(n0, n1, N0, N1);
92 swapped = 1; continue;
93 }
94 return gc_long(ltop,0);
95 }
96
97 m0 *= n_s;
98 a1 = (2 * p + 2) % m1;
99 x = (hasse_low + m0 - 1) / m0; /* using ceil(n/d) = (n + d - 1)/d */
100 x *= m0;
101 for ( ; x <= hasse_high; x += m0)
102 if ((x % m1) == a1 && x != N0 && x != N1) break;
103 /* every x in N was either N0 or N1, so we return true */
104 if (x > hasse_high) return gc_long(ltop,1);
105
106 lswap(a4, a4t);
107 lswap(a6, a6t);
108 lswap(m0, m1); set_avma(av);
109 }
110 }
111
112 static GEN
random_FleV(GEN x,GEN a6,ulong p,ulong pi)113 random_FleV(GEN x, GEN a6, ulong p, ulong pi)
114 { pari_APPLY_type(t_VEC, random_Fle_pre(uel(x,i), uel(a6,i), p, pi)) }
115
116 /**
117 * START Code from AVSs "torcosts.h"
118 */
119
120 struct torctab_rec {
121 int m;
122 int fix2, fix3;
123 int N;
124 int s2_flag;
125 int t3_flag;
126 double rating;
127 };
128
129 /*
130 These costs assume p=2 mod 3, 3 mod 4 and not 1 mod N
131 */
132
133 static struct torctab_rec torctab1[] = {
134 { 11, 1, 1, 11, 1, 1, 0.047250 },
135 { 33, 1, 0, 11, 1, 2, 0.047250 },
136 { 22, 1, 1, 11, 3, 1, 0.055125 },
137 { 66, 1, 0, 11, 3, 2, 0.055125 },
138 { 11, 1, 0, 11, 1, 0, 0.058000 },
139 { 13, 1, 1, 13, 1, 1, 0.058542 },
140 { 39, 1, 0, 13, 1, 2, 0.058542 },
141 { 22, 0, 1, 11, 2, 1, 0.061333 },
142 { 66, 0, 0, 11, 2, 2, 0.061333 },
143 { 22, 1, 0, 11, 3, 0, 0.061750 },
144 { 14, 1, 1, 14, 3, 1, 0.062500 },
145 { 42, 1, 0, 14, 3, 2, 0.062500 },
146 { 26, 1, 1, 13, 3, 1, 0.064583 },
147 { 78, 1, 0, 13, 3, 2, 0.064583 },
148 { 28, 0, 1, 14, 4, 1, 0.065625 },
149 { 84, 0, 0, 14, 4, 2, 0.065625 },
150 { 7, 1, 1, 7, 1, 1, 0.068750 },
151 { 13, 1, 0, 13, 1, 0, 0.068750 },
152 { 21, 1, 0, 7, 1, 2, 0.068750 },
153 { 26, 1, 0, 13, 3, 0, 0.069583 },
154 { 17, 1, 1, 17, 1, 1, 0.069687 },
155 { 51, 1, 0, 17, 1, 2, 0.069687 },
156 { 11, 0, 1, 11, 0, 1, 0.072500 },
157 { 33, 0, 0, 11, 0, 2, 0.072500 },
158 { 44, 1, 0, 11, 130, 0, 0.072667 },
159 { 52, 0, 1, 13, 4, 1, 0.073958 },
160 { 156, 0, 0, 13, 4, 2, 0.073958 },
161 { 34, 1, 1, 17, 3, 1, 0.075313 },
162 { 102, 1, 0, 17, 3, 2, 0.075313 },
163 { 15, 1, 0, 15, 1, 0, 0.075625 },
164 { 13, 0, 1, 13, 0, 1, 0.076667 },
165 { 39, 0, 0, 13, 0, 2, 0.076667 },
166 { 44, 0, 0, 11, 4, 0, 0.076667 },
167 { 30, 1, 0, 15, 3, 0, 0.077188 },
168 { 22, 0, 0, 11, 2, 0, 0.077333 },
169 { 34, 1, 0, 17, 3, 0, 0.077969 },
170 { 17, 1, 0, 17, 1, 0, 0.078750 },
171 { 14, 0, 1, 14, 0, 1, 0.080556 },
172 { 28, 0, 0, 14, 4, 0, 0.080556 },
173 { 42, 0, 0, 14, 0, 2, 0.080556 },
174 { 7, 1, 0, 7, 1, 0, 0.080833 },
175 { 9, 1, 0, 9, 1, 0, 0.080833 },
176 { 68, 0, 1, 17, 4, 1, 0.081380 },
177 { 204, 0, 0, 17, 4, 2, 0.081380 },
178 { 52, 0, 0, 13, 4, 0, 0.082292 },
179 { 10, 1, 1, 10, 3, 1, 0.084687 },
180 { 17, 0, 1, 17, 0, 1, 0.084687 },
181 { 51, 0, 0, 17, 0, 2, 0.084687 },
182 { 20, 0, 1, 10, 4, 1, 0.085938 },
183 { 60, 0, 0, 10, 4, 2, 0.085938 },
184 { 19, 1, 1, 19, 1, 1, 0.086111 },
185 { 57, 1, 0, 19, 1, 2, 0.086111 },
186 { 68, 0, 0, 17, 4, 0, 0.088281 },
187 { 38, 1, 1, 19, 3, 1, 0.089514 },
188 { 114, 1, 0, 19, 3, 2, 0.089514 },
189 { 20, 0, 0, 10, 4, 0, 0.090625 },
190 { 36, 0, 0, 18, 4, 0, 0.090972 },
191 { 26, 0, 0, 13, 2, 0, 0.091667 },
192 { 11, 0, 0, 11, 0, 0, 0.092000 },
193 { 19, 1, 0, 19, 1, 0, 0.092778 },
194 { 38, 1, 0, 19, 3, 0, 0.092778 },
195 { 14, 1, 0, 7, 3, 0, 0.092917 },
196 { 18, 1, 0, 9, 3, 0, 0.092917 },
197 { 76, 0, 1, 19, 4, 1, 0.095255 },
198 { 228, 0, 0, 19, 4, 2, 0.095255 },
199 { 10, 0, 1, 10, 0, 1, 0.096667 },
200 { 13, 0, 0, 13, 0, 0, 0.096667 },
201 { 30, 0, 0, 10, 0, 2, 0.096667 },
202 { 19, 0, 1, 19, 0, 1, 0.098333 },
203 { 57, 0, 0, 19, 0, 2, 0.098333 },
204 { 17, 0, 0, 17, 0, 0, 0.100000 },
205 { 23, 1, 1, 23, 1, 1, 0.100227 },
206 { 69, 1, 0, 23, 1, 2, 0.100227 },
207 { 7, 0, 1, 7, 0, 1, 0.100833 },
208 { 21, 0, 0, 7, 0, 2, 0.100833 },
209 { 76, 0, 0, 19, 4, 0, 0.102083 },
210 { 14, 0, 0, 14, 0, 0, 0.102222 },
211 { 18, 0, 0, 9, 2, 0, 0.102222 },
212 { 5, 1, 1, 5, 1, 1, 0.103125 },
213 { 46, 1, 1, 23, 3, 1, 0.104318 },
214 { 138, 1, 0, 23, 3, 2, 0.104318 },
215 { 23, 1, 0, 23, 1, 0, 0.105682 },
216 { 46, 1, 0, 23, 3, 0, 0.106705 },
217 { 92, 0, 1, 23, 4, 1, 0.109091 },
218 { 276, 0, 0, 23, 4, 2, 0.109091 },
219 { 19, 0, 0, 19, 0, 0, 0.110000 },
220 { 23, 0, 1, 23, 0, 1, 0.112273 },
221 { 69, 0, 0, 23, 0, 2, 0.112273 },
222 { 7, 0, 0, 7, 0, 0, 0.113333 },
223 { 9, 0, 0, 9, 0, 0, 0.113333 },
224 { 92, 0, 0, 23, 4, 0, 0.113826 },
225 { 16, 0, 1, 16, 0, 1, 0.118125 },
226 { 48, 0, 0, 16, 0, 2, 0.118125 },
227 { 5, 1, 0, 5, 1, 0, 0.121250 },
228 { 15, 0, 0, 15, 0, 0, 0.121250 },
229 { 10, 0, 0, 10, 0, 0, 0.121667 },
230 { 23, 0, 0, 23, 0, 0, 0.123182 },
231 { 12, 0, 0, 12, 0, 0, 0.141667 },
232 { 5, 0, 1, 5, 0, 1, 0.145000 },
233 { 16, 0, 0, 16, 0, 0, 0.145000 },
234 { 8, 0, 1, 8, 0, 1, 0.151250 },
235 { 29, 1, 1, 29, 1, 1, 0.153036 },
236 { 87, 1, 0, 29, 1, 2, 0.153036 },
237 { 25, 0, 0, 25, 0, 0, 0.155000 },
238 { 58, 1, 1, 29, 3, 1, 0.156116 },
239 { 174, 1, 0, 29, 3, 2, 0.156116 },
240 { 29, 1, 0, 29, 1, 0, 0.157500 },
241 { 58, 1, 0, 29, 3, 0, 0.157500 },
242 { 116, 0, 1, 29, 4, 1, 0.161086 },
243 { 29, 0, 1, 29, 0, 1, 0.163393 },
244 { 87, 0, 0, 29, 0, 2, 0.163393 },
245 { 116, 0, 0, 29, 4, 0, 0.163690 },
246 { 5, 0, 0, 5, 0, 0, 0.170000 },
247 { 8, 0, 0, 8, 0, 0, 0.170000 },
248 { 29, 0, 0, 29, 0, 0, 0.171071 },
249 { 31, 1, 1, 31, 1, 1, 0.186583 },
250 { 93, 1, 0, 31, 1, 2, 0.186583 },
251 { 62, 1, 1, 31, 3, 1, 0.189750 },
252 { 186, 1, 0, 31, 3, 2, 0.189750 },
253 { 31, 1, 0, 31, 1, 0, 0.191333 },
254 { 62, 1, 0, 31, 3, 0, 0.192167 },
255 { 124, 0, 1, 31, 4, 1, 0.193056 },
256 { 31, 0, 1, 31, 0, 1, 0.195333 },
257 { 93, 0, 0, 31, 0, 2, 0.195333 },
258 { 124, 0, 0, 31, 4, 0, 0.197917 },
259 { 2, 1, 1, 2, 3, 1, 0.200000 },
260 { 6, 1, 0, 2, 3, 2, 0.200000 },
261 { 31, 0, 0, 31, 0, 0, 0.206667 },
262 { 4, 1, 1, 4, 130, 1, 0.214167 },
263 { 6, 0, 0, 6, 0, 0, 0.226667 },
264 { 3, 1, 0, 3, 1, 0, 0.230000 },
265 { 4, 0, 1, 4, 0, 1, 0.241667 },
266 { 4, 1, 0, 2, 130, 0, 0.266667 },
267 { 4, 0, 0, 4, 0, 0, 0.283333 },
268 { 3, 0, 0, 3, 0, 0, 0.340000 },
269 { 1, 1, 1, 1, 1, 1, 0.362500 },
270 { 2, 0, 1, 2, 0, 1, 0.386667 },
271 { 1, 1, 0, 1, 1, 0, 0.410000 },
272 { 2, 0, 0, 2, 0, 0, 0.453333 },
273 };
274
275 static struct torctab_rec torctab2[] = {
276 { 11, 1, 1, 11, 1, 1, 0.047250 },
277 { 33, 1, 0, 11, 1, 2, 0.047250 },
278 { 22, 1, 1, 11, 3, 1, 0.055125 },
279 { 66, 1, 0, 11, 3, 2, 0.055125 },
280 { 13, 1, 1, 13, 1, 1, 0.057500 },
281 { 39, 1, 0, 13, 1, 2, 0.057500 },
282 { 11, 1, 0, 11, 1, 0, 0.058000 },
283 { 22, 0, 1, 11, 2, 1, 0.061333 },
284 { 66, 0, 0, 11, 2, 2, 0.061333 },
285 { 14, 1, 1, 14, 3, 1, 0.061458 },
286 { 42, 1, 0, 14, 3, 2, 0.061458 },
287 { 22, 1, 0, 11, 3, 0, 0.061750 },
288 { 26, 1, 1, 13, 3, 1, 0.064062 },
289 { 78, 1, 0, 13, 3, 2, 0.064062 },
290 { 28, 0, 1, 14, 4, 1, 0.065625 },
291 { 84, 0, 0, 14, 4, 2, 0.065625 },
292 { 13, 1, 0, 13, 1, 0, 0.066667 },
293 { 26, 1, 0, 13, 3, 0, 0.069583 },
294 { 17, 1, 1, 17, 1, 1, 0.069687 },
295 { 51, 1, 0, 17, 1, 2, 0.069687 },
296 { 11, 0, 1, 11, 0, 1, 0.070000 },
297 { 33, 0, 0, 11, 0, 2, 0.070000 },
298 { 7, 1, 1, 7, 1, 1, 0.070417 },
299 { 21, 1, 0, 7, 1, 2, 0.070417 },
300 { 15, 1, 0, 15, 1, 0, 0.072500 },
301 { 52, 0, 1, 13, 4, 1, 0.073090 },
302 { 156, 0, 0, 13, 4, 2, 0.073090 },
303 { 34, 1, 1, 17, 3, 1, 0.074219 },
304 { 102, 1, 0, 17, 3, 2, 0.074219 },
305 { 7, 1, 0, 7, 1, 0, 0.076667 },
306 { 13, 0, 1, 13, 0, 1, 0.076667 },
307 { 39, 0, 0, 13, 0, 2, 0.076667 },
308 { 44, 0, 0, 11, 4, 0, 0.076667 },
309 { 17, 1, 0, 17, 1, 0, 0.077188 },
310 { 22, 0, 0, 11, 2, 0, 0.077333 },
311 { 34, 1, 0, 17, 3, 0, 0.077969 },
312 { 30, 1, 0, 15, 3, 0, 0.080312 },
313 { 14, 0, 1, 14, 0, 1, 0.080556 },
314 { 28, 0, 0, 14, 4, 0, 0.080556 },
315 { 42, 0, 0, 14, 0, 2, 0.080556 },
316 { 9, 1, 0, 9, 1, 0, 0.080833 },
317 { 68, 0, 1, 17, 4, 1, 0.081380 },
318 { 204, 0, 0, 17, 4, 2, 0.081380 },
319 { 52, 0, 0, 13, 4, 0, 0.082292 },
320 { 10, 1, 1, 10, 3, 1, 0.083125 },
321 { 20, 0, 1, 10, 4, 1, 0.083333 },
322 { 60, 0, 0, 10, 4, 2, 0.083333 },
323 { 17, 0, 1, 17, 0, 1, 0.084687 },
324 { 51, 0, 0, 17, 0, 2, 0.084687 },
325 { 19, 1, 1, 19, 1, 1, 0.084722 },
326 { 57, 1, 0, 19, 1, 2, 0.084722 },
327 { 11, 0, 0, 11, 0, 0, 0.087000 },
328 { 68, 0, 0, 17, 4, 0, 0.088281 },
329 { 38, 1, 1, 19, 3, 1, 0.090139 },
330 { 114, 1, 0, 19, 3, 2, 0.090139 },
331 { 36, 0, 0, 18, 4, 0, 0.090972 },
332 { 19, 1, 0, 19, 1, 0, 0.091389 },
333 { 26, 0, 0, 13, 2, 0, 0.091667 },
334 { 13, 0, 0, 13, 0, 0, 0.092500 },
335 { 38, 1, 0, 19, 3, 0, 0.092778 },
336 { 14, 1, 0, 7, 3, 0, 0.092917 },
337 { 18, 1, 0, 9, 3, 0, 0.092917 },
338 { 20, 0, 0, 10, 4, 0, 0.095833 },
339 { 76, 0, 1, 19, 4, 1, 0.096412 },
340 { 228, 0, 0, 19, 4, 2, 0.096412 },
341 { 17, 0, 0, 17, 0, 0, 0.096875 },
342 { 19, 0, 1, 19, 0, 1, 0.098056 },
343 { 57, 0, 0, 19, 0, 2, 0.098056 },
344 { 23, 1, 1, 23, 1, 1, 0.100682 },
345 { 69, 1, 0, 23, 1, 2, 0.100682 },
346 { 7, 0, 1, 7, 0, 1, 0.100833 },
347 { 21, 0, 0, 7, 0, 2, 0.100833 },
348 { 30, 0, 0, 15, 2, 0, 0.100833 },
349 { 76, 0, 0, 19, 4, 0, 0.102083 },
350 { 14, 0, 0, 14, 0, 0, 0.102222 },
351 { 5, 1, 1, 5, 1, 1, 0.103125 },
352 { 46, 1, 1, 23, 3, 1, 0.104034 },
353 { 138, 1, 0, 23, 3, 2, 0.104034 },
354 { 23, 1, 0, 23, 1, 0, 0.104545 },
355 { 7, 0, 0, 7, 0, 0, 0.105000 },
356 { 10, 0, 1, 10, 0, 1, 0.105000 },
357 { 16, 0, 1, 16, 0, 1, 0.105417 },
358 { 48, 0, 0, 16, 0, 2, 0.105417 },
359 { 46, 1, 0, 23, 3, 0, 0.106705 },
360 { 18, 0, 0, 9, 2, 0, 0.107778 },
361 { 92, 0, 1, 23, 4, 1, 0.108239 },
362 { 276, 0, 0, 23, 4, 2, 0.108239 },
363 { 19, 0, 0, 19, 0, 0, 0.110000 },
364 { 23, 0, 1, 23, 0, 1, 0.111136 },
365 { 69, 0, 0, 23, 0, 2, 0.111136 },
366 { 9, 0, 0, 9, 0, 0, 0.113333 },
367 { 10, 0, 0, 10, 0, 0, 0.113333 },
368 { 92, 0, 0, 23, 4, 0, 0.113826 },
369 { 5, 1, 0, 5, 1, 0, 0.115000 },
370 { 15, 0, 0, 15, 0, 0, 0.115000 },
371 { 23, 0, 0, 23, 0, 0, 0.120909 },
372 { 8, 0, 1, 8, 0, 1, 0.126042 },
373 { 24, 0, 0, 8, 0, 2, 0.126042 },
374 { 16, 0, 0, 16, 0, 0, 0.127188 },
375 { 8, 0, 0, 8, 0, 0, 0.141667 },
376 { 25, 0, 1, 25, 0, 1, 0.144000 },
377 { 5, 0, 1, 5, 0, 1, 0.151250 },
378 { 12, 0, 0, 12, 0, 0, 0.152083 },
379 { 29, 1, 1, 29, 1, 1, 0.153929 },
380 { 87, 1, 0, 29, 1, 2, 0.153929 },
381 { 25, 0, 0, 25, 0, 0, 0.155000 },
382 { 58, 1, 1, 29, 3, 1, 0.155045 },
383 { 174, 1, 0, 29, 3, 2, 0.155045 },
384 { 29, 1, 0, 29, 1, 0, 0.156429 },
385 { 58, 1, 0, 29, 3, 0, 0.157857 },
386 { 116, 0, 1, 29, 4, 1, 0.158631 },
387 { 116, 0, 0, 29, 4, 0, 0.163542 },
388 { 29, 0, 1, 29, 0, 1, 0.164286 },
389 { 87, 0, 0, 29, 0, 2, 0.164286 },
390 { 29, 0, 0, 29, 0, 0, 0.169286 },
391 { 5, 0, 0, 5, 0, 0, 0.170000 },
392 { 31, 1, 1, 31, 1, 1, 0.187000 },
393 { 93, 1, 0, 31, 1, 2, 0.187000 },
394 { 62, 1, 1, 31, 3, 1, 0.188500 },
395 { 186, 1, 0, 31, 3, 2, 0.188500 },
396 { 31, 1, 0, 31, 1, 0, 0.191333 },
397 { 62, 1, 0, 31, 3, 0, 0.192083 },
398 { 124, 0, 1, 31, 4, 1, 0.193472 },
399 { 31, 0, 1, 31, 0, 1, 0.196167 },
400 { 93, 0, 0, 31, 0, 2, 0.196167 },
401 { 124, 0, 0, 31, 4, 0, 0.197083 },
402 { 2, 1, 1, 2, 3, 1, 0.200000 },
403 { 6, 1, 0, 2, 3, 2, 0.200000 },
404 { 31, 0, 0, 31, 0, 0, 0.205000 },
405 { 6, 0, 0, 6, 0, 0, 0.226667 },
406 { 3, 1, 0, 3, 1, 0, 0.230000 },
407 { 4, 0, 1, 4, 0, 1, 0.241667 },
408 { 4, 0, 0, 4, 0, 0, 0.283333 },
409 { 3, 0, 0, 3, 0, 0, 0.340000 },
410 { 1, 1, 1, 1, 1, 1, 0.362500 },
411 { 2, 0, 1, 2, 0, 1, 0.370000 },
412 { 1, 1, 0, 1, 1, 0, 0.385000 },
413 { 2, 0, 0, 2, 0, 0, 0.453333 },
414 };
415
416 static struct torctab_rec torctab3[] = {
417 { 66, 1, 0, 11, 3, 2, 0.040406 },
418 { 33, 1, 0, 11, 1, 2, 0.043688 },
419 { 78, 1, 0, 13, 3, 2, 0.045391 },
420 { 132, 1, 0, 11, 130, 2, 0.046938 },
421 { 39, 1, 0, 13, 1, 2, 0.047656 },
422 { 102, 1, 0, 17, 3, 2, 0.049922 },
423 { 42, 1, 0, 14, 3, 2, 0.050000 },
424 { 51, 1, 0, 17, 1, 2, 0.051680 },
425 { 132, 0, 0, 11, 4, 2, 0.052188 },
426 { 156, 1, 0, 13, 130, 2, 0.053958 },
427 { 156, 0, 0, 13, 4, 2, 0.054818 },
428 { 84, 1, 0, 14, 130, 2, 0.055000 },
429 { 15, 1, 0, 15, 1, 0, 0.056719 },
430 { 204, 0, 0, 17, 4, 2, 0.057227 },
431 { 114, 1, 0, 19, 3, 2, 0.057500 },
432 { 11, 1, 0, 11, 1, 0, 0.058000 },
433 { 66, 0, 0, 11, 2, 2, 0.058000 },
434 { 57, 1, 0, 19, 1, 2, 0.059062 },
435 { 30, 1, 0, 15, 3, 0, 0.059063 },
436 { 84, 0, 0, 14, 4, 2, 0.060677 },
437 { 22, 1, 0, 11, 3, 0, 0.061750 },
438 { 78, 0, 0, 13, 2, 2, 0.063542 },
439 { 228, 0, 0, 19, 4, 2, 0.063889 },
440 { 21, 1, 0, 7, 1, 2, 0.065000 },
441 { 138, 1, 0, 23, 3, 2, 0.065028 },
442 { 69, 1, 0, 23, 1, 2, 0.066903 },
443 { 13, 1, 0, 13, 1, 0, 0.068750 },
444 { 102, 0, 0, 17, 2, 2, 0.068906 },
445 { 26, 1, 0, 13, 3, 0, 0.069583 },
446 { 51, 0, 0, 17, 0, 2, 0.070312 },
447 { 60, 1, 0, 15, 130, 0, 0.071094 },
448 { 276, 0, 0, 23, 4, 2, 0.071236 },
449 { 39, 0, 0, 13, 0, 2, 0.071250 },
450 { 33, 0, 0, 11, 0, 2, 0.072750 },
451 { 44, 1, 0, 11, 130, 0, 0.073500 },
452 { 60, 0, 0, 15, 4, 0, 0.073828 },
453 { 9, 1, 0, 9, 1, 0, 0.074097 },
454 { 30, 0, 0, 15, 2, 0, 0.075625 },
455 { 57, 0, 0, 19, 0, 2, 0.075625 },
456 { 7, 1, 0, 7, 1, 0, 0.076667 },
457 { 44, 0, 0, 11, 4, 0, 0.076667 },
458 { 22, 0, 0, 11, 2, 0, 0.077333 },
459 { 17, 1, 0, 17, 1, 0, 0.078750 },
460 { 34, 1, 0, 17, 3, 0, 0.078750 },
461 { 69, 0, 0, 23, 0, 2, 0.079943 },
462 { 28, 0, 0, 14, 4, 0, 0.080556 },
463 { 42, 0, 0, 14, 0, 2, 0.080833 },
464 { 52, 0, 0, 13, 4, 0, 0.082292 },
465 { 14, 1, 1, 14, 3, 1, 0.083333 },
466 { 36, 0, 0, 18, 4, 0, 0.083391 },
467 { 18, 1, 0, 9, 3, 0, 0.085174 },
468 { 68, 0, 0, 17, 4, 0, 0.089583 },
469 { 15, 0, 0, 15, 0, 0, 0.090938 },
470 { 19, 1, 0, 19, 1, 0, 0.091389 },
471 { 26, 0, 0, 13, 2, 0, 0.091667 },
472 { 11, 0, 0, 11, 0, 0, 0.092000 },
473 { 13, 0, 0, 13, 0, 0, 0.092500 },
474 { 38, 1, 0, 19, 3, 0, 0.092778 },
475 { 14, 1, 0, 7, 3, 0, 0.092917 },
476 { 18, 0, 0, 9, 2, 0, 0.093704 },
477 { 174, 1, 0, 29, 3, 2, 0.095826 },
478 { 20, 0, 0, 10, 4, 0, 0.095833 },
479 { 96, 1, 0, 16, 133, 2, 0.096562 },
480 { 21, 0, 0, 21, 0, 0, 0.096875 },
481 { 87, 1, 0, 29, 1, 2, 0.096964 },
482 { 17, 0, 0, 17, 0, 0, 0.100000 },
483 { 348, 0, 0, 29, 4, 2, 0.100558 },
484 { 76, 0, 0, 19, 4, 0, 0.100926 },
485 { 14, 0, 0, 14, 0, 0, 0.102222 },
486 { 9, 0, 0, 9, 0, 0, 0.103889 },
487 { 46, 1, 0, 23, 3, 0, 0.105114 },
488 { 23, 1, 0, 23, 1, 0, 0.105682 },
489 { 48, 0, 0, 16, 0, 2, 0.106406 },
490 { 87, 0, 0, 29, 0, 2, 0.107545 },
491 { 19, 0, 0, 19, 0, 0, 0.107778 },
492 { 7, 0, 0, 7, 0, 0, 0.113333 },
493 { 10, 0, 0, 10, 0, 0, 0.113333 },
494 { 92, 0, 0, 23, 4, 0, 0.113636 },
495 { 12, 0, 0, 12, 0, 0, 0.114062 },
496 { 5, 1, 0, 5, 1, 0, 0.115000 },
497 { 186, 1, 0, 31, 3, 2, 0.115344 },
498 { 93, 1, 0, 31, 1, 2, 0.118125 },
499 { 23, 0, 0, 23, 0, 0, 0.120909 },
500 { 93, 0, 0, 31, 0, 2, 0.128250 },
501 { 16, 0, 0, 16, 0, 0, 0.138750 },
502 { 25, 0, 0, 25, 0, 0, 0.155000 },
503 { 58, 1, 0, 29, 3, 0, 0.155714 },
504 { 29, 1, 0, 29, 1, 0, 0.158214 },
505 { 3, 1, 0, 3, 1, 0, 0.163125 },
506 { 116, 0, 0, 29, 4, 0, 0.163690 },
507 { 5, 0, 0, 5, 0, 0, 0.170000 },
508 { 6, 0, 0, 6, 0, 0, 0.170000 },
509 { 8, 0, 0, 8, 0, 0, 0.170000 },
510 { 29, 0, 0, 29, 0, 0, 0.172857 },
511 { 31, 1, 0, 31, 1, 0, 0.191333 },
512 { 62, 1, 0, 31, 3, 0, 0.191750 },
513 { 124, 0, 0, 31, 4, 0, 0.197917 },
514 { 31, 0, 0, 31, 0, 0, 0.201667 },
515 { 3, 0, 0, 3, 0, 0, 0.236250 },
516 { 4, 0, 0, 4, 0, 0, 0.262500 },
517 { 2, 1, 1, 2, 3, 1, 0.317187 },
518 { 1, 1, 0, 1, 1, 0, 0.410000 },
519 { 2, 0, 0, 2, 0, 0, 0.453333 },
520 };
521
522 static struct torctab_rec torctab4[] = {
523 { 66, 1, 0, 11, 3, 2, 0.041344 },
524 { 33, 1, 0, 11, 1, 2, 0.042750 },
525 { 78, 1, 0, 13, 3, 2, 0.045781 },
526 { 39, 1, 0, 13, 1, 2, 0.046875 },
527 { 264, 1, 0, 11, 131, 2, 0.049043 },
528 { 42, 1, 0, 14, 3, 2, 0.050000 },
529 { 102, 1, 0, 17, 3, 2, 0.050508 },
530 { 51, 1, 0, 17, 1, 2, 0.051094 },
531 { 528, 1, 0, 11, 132, 2, 0.052891 },
532 { 132, 0, 0, 11, 4, 2, 0.052969 },
533 { 168, 1, 0, 14, 131, 2, 0.053965 },
534 { 156, 0, 0, 13, 4, 2, 0.054948 },
535 { 336, 1, 0, 14, 132, 2, 0.056120 },
536 { 15, 1, 0, 15, 1, 0, 0.056719 },
537 { 66, 0, 0, 11, 2, 2, 0.057000 },
538 { 114, 1, 0, 19, 3, 2, 0.057812 },
539 { 11, 1, 0, 11, 1, 0, 0.058000 },
540 { 204, 0, 0, 17, 4, 2, 0.058203 },
541 { 57, 1, 0, 19, 1, 2, 0.058542 },
542 { 84, 0, 0, 14, 4, 2, 0.059375 },
543 { 30, 1, 0, 15, 3, 0, 0.061406 },
544 { 22, 1, 0, 11, 3, 0, 0.063000 },
545 { 78, 0, 0, 13, 2, 2, 0.063542 },
546 { 138, 1, 0, 23, 3, 2, 0.064815 },
547 { 21, 1, 0, 7, 1, 2, 0.065000 },
548 { 228, 0, 0, 19, 4, 2, 0.065104 },
549 { 69, 1, 0, 23, 1, 2, 0.066477 },
550 { 13, 1, 0, 13, 1, 0, 0.068750 },
551 { 102, 0, 0, 17, 2, 2, 0.068906 },
552 { 51, 0, 0, 17, 0, 2, 0.069141 },
553 { 26, 1, 0, 13, 3, 0, 0.070625 },
554 { 276, 0, 0, 23, 4, 2, 0.071236 },
555 { 39, 0, 0, 13, 0, 2, 0.071250 },
556 { 33, 0, 0, 11, 0, 2, 0.072750 },
557 { 60, 0, 0, 15, 4, 0, 0.073828 },
558 { 9, 1, 0, 9, 1, 0, 0.074097 },
559 { 57, 0, 0, 19, 0, 2, 0.074583 },
560 { 30, 0, 0, 15, 2, 0, 0.075625 },
561 { 44, 0, 0, 11, 4, 0, 0.076667 },
562 { 17, 1, 0, 17, 1, 0, 0.077188 },
563 { 22, 0, 0, 11, 2, 0, 0.077333 },
564 { 69, 0, 0, 23, 0, 2, 0.080114 },
565 { 36, 0, 0, 18, 4, 0, 0.080208 },
566 { 34, 1, 0, 17, 3, 0, 0.080312 },
567 { 28, 0, 0, 14, 4, 0, 0.080556 },
568 { 7, 1, 0, 7, 1, 0, 0.080833 },
569 { 52, 0, 0, 13, 4, 0, 0.082292 },
570 { 42, 0, 0, 14, 0, 2, 0.082500 },
571 { 14, 1, 1, 14, 3, 1, 0.083333 },
572 { 15, 0, 0, 15, 0, 0, 0.086250 },
573 { 18, 1, 0, 9, 3, 0, 0.087083 },
574 { 26, 0, 0, 13, 2, 0, 0.088889 },
575 { 68, 0, 0, 17, 4, 0, 0.089583 },
576 { 48, 1, 0, 16, 132, 2, 0.089844 },
577 { 19, 1, 0, 19, 1, 0, 0.091389 },
578 { 11, 0, 0, 11, 0, 0, 0.092000 },
579 { 38, 1, 0, 19, 3, 0, 0.092917 },
580 { 18, 0, 0, 9, 2, 0, 0.093704 },
581 { 14, 1, 0, 7, 3, 0, 0.095000 },
582 { 96, 1, 0, 16, 133, 2, 0.095391 },
583 { 20, 0, 0, 10, 4, 0, 0.095833 },
584 { 174, 1, 0, 29, 3, 2, 0.095893 },
585 { 13, 0, 0, 13, 0, 0, 0.096667 },
586 { 17, 0, 0, 17, 0, 0, 0.096875 },
587 { 21, 0, 0, 21, 0, 0, 0.096875 },
588 { 87, 1, 0, 29, 1, 2, 0.097366 },
589 { 48, 0, 0, 16, 0, 2, 0.097969 },
590 { 24, 1, 0, 12, 131, 0, 0.098789 },
591 { 76, 0, 0, 19, 4, 0, 0.100926 },
592 { 348, 0, 0, 29, 4, 2, 0.101116 },
593 { 14, 0, 0, 14, 0, 0, 0.102222 },
594 { 9, 0, 0, 9, 0, 0, 0.103889 },
595 { 23, 1, 0, 23, 1, 0, 0.104545 },
596 { 46, 1, 0, 23, 3, 0, 0.105682 },
597 { 12, 0, 0, 12, 0, 0, 0.106250 },
598 { 87, 0, 0, 29, 0, 2, 0.108348 },
599 { 19, 0, 0, 19, 0, 0, 0.110000 },
600 { 7, 0, 0, 7, 0, 0, 0.113333 },
601 { 10, 0, 0, 10, 0, 0, 0.113333 },
602 { 92, 0, 0, 23, 4, 0, 0.113826 },
603 { 186, 1, 0, 31, 3, 2, 0.116094 },
604 { 93, 1, 0, 31, 1, 2, 0.116813 },
605 { 23, 0, 0, 23, 0, 0, 0.120909 },
606 { 5, 1, 0, 5, 1, 0, 0.121250 },
607 { 93, 0, 0, 31, 0, 2, 0.127625 },
608 { 16, 0, 0, 16, 0, 0, 0.132917 },
609 { 8, 0, 0, 8, 0, 0, 0.141667 },
610 { 25, 0, 0, 25, 0, 0, 0.152500 },
611 { 58, 1, 0, 29, 3, 0, 0.157946 },
612 { 29, 1, 0, 29, 1, 0, 0.158393 },
613 { 116, 0, 0, 29, 4, 0, 0.162946 },
614 { 3, 1, 0, 3, 1, 0, 0.163125 },
615 { 29, 0, 0, 29, 0, 0, 0.169286 },
616 { 5, 0, 0, 5, 0, 0, 0.170000 },
617 { 6, 0, 0, 6, 0, 0, 0.170000 },
618 { 31, 1, 0, 31, 1, 0, 0.191333 },
619 { 62, 1, 0, 31, 3, 0, 0.192083 },
620 { 124, 0, 0, 31, 4, 0, 0.196389 },
621 { 31, 0, 0, 31, 0, 0, 0.205000 },
622 { 3, 0, 0, 3, 0, 0, 0.255000 },
623 { 4, 0, 0, 4, 0, 0, 0.262500 },
624 { 2, 1, 1, 2, 3, 1, 0.325000 },
625 { 1, 1, 0, 1, 1, 0, 0.385000 },
626 { 2, 0, 0, 2, 0, 0, 0.420000 },
627 };
628
629 #define TWIST_DOUBLE_RATIO (9.0/16.0)
630
631 static long
torsion_constraint(struct torctab_rec * torctab,long ltorc,double tormod[],long n,long m)632 torsion_constraint(struct torctab_rec *torctab, long ltorc, double tormod[], long n, long m)
633 {
634 long i, b = -1;
635 double rb = -1.;
636 for (i = 0 ; i < ltorc ; i++)
637 {
638 struct torctab_rec *ti = torctab + i;
639 if ( ! (n%ti->m) && ( !ti->fix2 || (n%(2*ti->m)) ) && ( ! ti->fix3 || (n%(3*ti->m)) ) )
640 if ( n == m || ( ! (m%ti->m) && ( !ti->fix2 || (m%(2*ti->m)) ) && ( ! ti->fix3 || (m%(3*ti->m)) ) ) )
641 {
642 double ri = ti->rating*tormod[ti->N];
643 if ( b < 0 || ri < rb ) { b = i; rb = ri; }
644 }
645 }
646 if (b < 0) pari_err_BUG("find_rating");
647 return b;
648 }
649
650 static void
best_torsion_constraint(ulong p,long t,int * ptwist,ulong * ptor,int * ps2,int * pt3)651 best_torsion_constraint(ulong p, long t, int *ptwist, ulong *ptor, int *ps2, int *pt3)
652 {
653 struct torctab_rec *torctab;
654 double tormod[32];
655 long ltorc;
656 long n1, n2;
657 long b, b1, b2, b12;
658 long i;
659
660 if ( (p%3)==2 ) {
661 if ( (p&3)==3 ) {
662 torctab = torctab1;
663 ltorc = sizeof(torctab1)/sizeof(*torctab1);
664 } else {
665 torctab = torctab2;
666 ltorc = sizeof(torctab2)/sizeof(*torctab2);
667 }
668 } else {
669 if ( (p&3)==3 ) {
670 torctab = torctab3;
671 ltorc = sizeof(torctab3)/sizeof(*torctab3);
672 } else {
673 torctab = torctab4;
674 ltorc = sizeof(torctab4)/sizeof(*torctab4);
675 }
676 }
677 for ( i = 0 ; i < 32 ; i++ ) tormod[i] = 1.0;
678 if ( (p%5)==1 ) tormod[5] = tormod[10] = tormod[15] = 6.0/5.0;
679 if ( (p%7)==1 ) tormod[7] = tormod[14] = 8.0/7.0;
680 if ( (p%11)== 1 ) tormod[11] = 12.0/11.0;
681 if ( (p%13)==1 ) tormod[13] = 14.0/13.0;
682 if ( (p%17)==1 ) tormod[17] = 18.0/17.0;
683 if ( (p%19)==1 ) tormod[19] = 20.0/19.0;
684 if ( (p%23)==1 ) tormod[23] = 24.0/23.0;
685 if ( (p%29)==1 ) tormod[29] = 30.0/29.0;
686 if ( (p%31)==1 ) tormod[31] = 32.0/31.0;
687
688 n1 = p+1-t;
689 n2 = p+1+t;
690 b12 = -1;
691 b1 = torsion_constraint(torctab, ltorc, tormod, n1, n1);
692 b2 = torsion_constraint(torctab, ltorc, tormod, n2, n2);
693 b12 = torsion_constraint(torctab, ltorc, tormod, n1, n2);
694 if ( b1 > b2 ) {
695 if ( torctab[b2].rating / TWIST_DOUBLE_RATIO > torctab[b12].rating )
696 *ptwist = 3;
697 else
698 *ptwist = 2;
699 } else
700 if ( torctab[b1].rating / TWIST_DOUBLE_RATIO > torctab[b12].rating )
701 *ptwist = 3;
702 else
703 *ptwist = 1;
704 b = *ptwist ==1 ? b1: *ptwist ==2 ? b2: b12;
705 *ptor = torctab[b].N; *ps2 = torctab[b].s2_flag; *pt3 = torctab[b].t3_flag;
706 }
707
708 /* This is Sutherland 2009 Algorithm 1.1 */
709 static long
find_j_inv_with_given_trace(ulong * j_t,norm_eqn_t ne,long rho_inv,long max_curves)710 find_j_inv_with_given_trace(
711 ulong *j_t, norm_eqn_t ne, long rho_inv, long max_curves)
712 {
713 pari_sp ltop = avma, av;
714 long curves_tested = 0, batch_size;
715 long N0, N1, hasse[2];
716 GEN n0, n1;
717 long i, found = 0;
718 ulong p = ne->p, pi = ne->pi;
719 long t = ne->t;
720 ulong p1 = p + 1, a4, a6, m, N;
721 GEN A4, A6, tx, ty;
722 int s2_flag, t3_flag, twist;
723
724 if (p == 2 || p == 3) {
725 if (t == 0) pari_err_BUG("find_j_inv_with_given_trace");
726 *j_t = t; return 1;
727 }
728
729 N0 = (long)p1 - t; n0 = factoru(N0);
730 N1 = (long)p1 + t; n1 = factoru(N1);
731
732 best_torsion_constraint(p, t, &twist, &m, &s2_flag, &t3_flag);
733 N = p1 - (twist<3 ? (twist==1 ? t: -t): 0);
734
735 /* Select batch size so that we have roughly a 50% chance of finding
736 * a good curve in a batch. */
737 batch_size = 1.0 + rho_inv / (2.0 * m);
738 A4 = cgetg(batch_size + 1, t_VECSMALL);
739 A6 = cgetg(batch_size + 1, t_VECSMALL);
740 tx = cgetg(batch_size + 1, t_VECSMALL);
741 ty = cgetg(batch_size + 1, t_VECSMALL);
742
743 dbg_printf(2)(" Selected torsion constraint m = %lu and batch "
744 "size = %ld\n", m, batch_size);
745 hasse_bounds(&hasse[0], &hasse[1], p);
746 av = avma;
747 while (!found && (max_curves <= 0 || curves_tested < max_curves))
748 {
749 GEN Pp1, Pt;
750 random_curves_with_m_torsion((ulong *)(A4 + 1), (ulong *)(A6 + 1),
751 (ulong *)(tx + 1), (ulong *)(ty + 1),
752 batch_size, m, p);
753 Pp1 = random_FleV(A4, A6, p, pi);
754 Pt = gcopy(Pp1);
755 FleV_mulu_pre_inplace(Pp1, N, A4, p, pi);
756 if (twist >= 3) FleV_mulu_pre_inplace(Pt, t, A4, p, pi);
757 for (i = 1; i <= batch_size; ++i) {
758 ++curves_tested;
759 a4 = A4[i];
760 a6 = A6[i]; if (a4 == 0 || a6 == 0) continue;
761
762 if (( (twist >= 3 && mael(Pp1,i,1) == mael(Pt,i,1))
763 || (twist < 3 && umael(Pp1,i,1) == p))
764 && test_curve_order(ne, a4, a6, N0, N1, n0, n1, hasse)) {
765 *j_t = Fl_ellj_pre(a4, a6, p, pi);
766 found = 1; break;
767 }
768 }
769 set_avma(av);
770 }
771 return gc_long(ltop, curves_tested);
772 }
773
774 /*
775 * SECTION: Functions for dealing with polycyclic presentations.
776 */
777
778 static GEN
next_generator(GEN DD,long D,ulong u,long filter,GEN * genred,long * P)779 next_generator(GEN DD, long D, ulong u, long filter, GEN *genred, long *P)
780 {
781 pari_sp av = avma;
782 ulong p = (ulong)*P;
783 while (1)
784 {
785 p = unextprime(p + 1);
786 if (p > LONG_MAX) pari_err_BUG("next_generator");
787 if (kross(D, (long)p) != -1 && u % p != 0 && filter % p != 0)
788 {
789 GEN gen = primeform_u(DD, p);
790 /* If gen is in the principal class, skip it */
791 *genred = redimag(gen);
792 if (!equali1(gel(*genred,1))) { *P = (long)p; return gen; }
793 set_avma(av);
794 }
795 }
796 }
797
798 INLINE long *
evec_ri_mutate(long r[],long i)799 evec_ri_mutate(long r[], long i)
800 { return r + (i * (i - 1) >> 1); }
801
802 INLINE const long *
evec_ri(const long r[],long i)803 evec_ri(const long r[], long i)
804 { return r + (i * (i - 1) >> 1); }
805
806 /* Reduces evec e so that e[i] < n[i] (assume e[i] >= 0) using pcp(n,r,k).
807 * No check for overflow, this could be an issue for large groups */
808 INLINE void
evec_reduce(long e[],const long n[],const long r[],long k)809 evec_reduce(long e[], const long n[], const long r[], long k)
810 {
811 long i, j, q;
812 const long *ri;
813 if (!k) return;
814 for (i = k - 1; i > 0; i--) {
815 if (e[i] >= n[i]) {
816 q = e[i] / n[i];
817 ri = evec_ri(r, i);
818 for (j = 0; j < i; j++) e[j] += q * ri[j];
819 e[i] -= q * n[i];
820 }
821 }
822 e[0] %= n[0];
823 }
824
825 /* Computes e3 = log(a^e1*a^e2) in terms of the given polycyclic
826 * presentation (here a denotes the implicit vector of generators) */
827 INLINE void
evec_compose(long e3[],const long e1[],const long e2[],const long n[],const long r[],long k)828 evec_compose(long e3[],
829 const long e1[], const long e2[], const long n[],const long r[], long k)
830 {
831 long i;
832 for (i = 0; i < k; i++) e3[i] = e1[i] + e2[i];
833 evec_reduce(e3, n, r, k);
834 }
835
836 /* Converts an evec to an integer index corresponding to the
837 * multi-radix representation of the evec with moduli corresponding to
838 * the subgroup orders m[i] */
839 INLINE long
evec_to_index(const long e[],const long m[],long k)840 evec_to_index(const long e[], const long m[], long k)
841 {
842 long i, index = e[0];
843 for (i = 1; i < k; i++) index += e[i] * m[i - 1];
844 return index;
845 }
846
847 INLINE void
evec_copy(long f[],const long e[],long k)848 evec_copy(long f[], const long e[], long k)
849 {
850 long i;
851 for (i = 0; i < k; ++i) f[i] = e[i];
852 }
853
854 INLINE void
evec_clear(long e[],long k)855 evec_clear(long e[], long k)
856 {
857 long i;
858 for (i = 0; i < k; ++i) e[i] = 0;
859 }
860
861 /* e1 and e2 may overlap */
862 /* Note that this function is not very efficient because it does not know the
863 * orders of the elements in the presentation, only the relative orders */
864 INLINE void
evec_inverse(long e2[],const long e1[],const long n[],const long r[],long k)865 evec_inverse(long e2[], const long e1[], const long n[], const long r[], long k)
866 {
867 pari_sp av = avma;
868 long i, *e3, *e4;
869
870 e3 = new_chunk(k);
871 e4 = new_chunk(k);
872 evec_clear(e4, k);
873 evec_copy(e3, e1, k);
874 /* We have e1 + e4 = e3 which we maintain throughout while making e1
875 * the zero vector */
876 for (i = k - 1; i >= 0; i--) if (e3[i])
877 {
878 e4[i] += n[i] - e3[i];
879 evec_reduce(e4, n, r, k);
880 e3[i] = n[i];
881 evec_reduce(e3, n, r, k);
882 }
883 evec_copy(e2, e4, k);
884 set_avma(av);
885 }
886
887 /* e1 and e2 may overlap */
888 /* This is a faster way to compute inverses, if the presentation
889 * element orders are known (these are specified in the array o, the
890 * array n holds the relative orders) */
891 INLINE void
evec_inverse_o(long e2[],const long e1[],const long n[],const long o[],const long r[],long k)892 evec_inverse_o(
893 long e2[],
894 const long e1[], const long n[], const long o[], const long r[], long k)
895 {
896 long j;
897 for (j = 0; j < k; j++) e2[j] = (e1[j] ? o[j] - e1[j] : 0);
898 evec_reduce(e2, n, r, k);
899 }
900
901 /* Computes the order of the group element a^e using the pcp (n,r,k) */
902 INLINE long
evec_order(const long e[],const long n[],const long r[],long k)903 evec_order(const long e[], const long n[], const long r[], long k)
904 {
905 pari_sp av = avma;
906 long *f = new_chunk(k);
907 long i, j, o, m;
908
909 evec_copy(f, e, k);
910 for (o = 1, i = k - 1; i >= 0; i--) if (f[i])
911 {
912 m = n[i] / ugcd(f[i], n[i]);
913 for (j = 0; j < k; j++) f[j] *= m;
914 evec_reduce(f, n, r, k);
915 o *= m;
916 }
917 return gc_long(av,o);
918 }
919
920 /* Computes orders o[] for each generator using relative orders n[]
921 * and power relations r[] */
922 INLINE void
evec_orders(long o[],const long n[],const long r[],long k)923 evec_orders(long o[], const long n[], const long r[], long k)
924 {
925 pari_sp av = avma;
926 long i, *e = new_chunk(k);
927
928 evec_clear(e, k);
929 for (i = 0; i < k; i++) {
930 e[i] = 1;
931 if (i) e[i - 1] = 0;
932 o[i] = evec_order(e, n, r, k);
933 }
934 set_avma(av);
935 }
936
937 INLINE int
evec_equal(const long e1[],const long e2[],long k)938 evec_equal(const long e1[], const long e2[], long k)
939 {
940 long j;
941 for (j = 0; j < k; ++j)
942 if (e1[j] != e2[j]) break;
943 return j == k;
944 }
945
946 INLINE void
index_to_evec(long e[],long index,const long m[],long k)947 index_to_evec(long e[], long index, const long m[], long k)
948 {
949 long i;
950 for (i = k - 1; i > 0; --i) {
951 e[i] = index / m[i - 1];
952 index -= e[i] * m[i - 1];
953 }
954 e[0] = index;
955 }
956
957 INLINE void
evec_n_to_m(long m[],const long n[],long k)958 evec_n_to_m(long m[], const long n[], long k)
959 {
960 long i;
961 m[0] = n[0];
962 for (i = 1; i < k; ++i) m[i] = m[i - 1] * n[i];
963 }
964
965 /* Based on logfac() in Sutherland's classpoly package.
966 * Ramanujan approximation to log(n!), accurate to O(1/n^3) */
967 INLINE double
logfac(long n)968 logfac(long n)
969 {
970 const double HALFLOGPI = 0.57236494292470008707171367567653;
971 return n * log((double) n) - (double) n +
972 log((double) n * (1.0 + 4.0 * n * (1.0 + 2.0 * n))) / 6.0 +
973 HALFLOGPI;
974 }
975
976 /* This is based on Sutherland 2009, Lemma 8 (p31). */
977 static double
upper_bound_on_classpoly_coeffs(long D,long h,GEN qfinorms)978 upper_bound_on_classpoly_coeffs(long D, long h, GEN qfinorms)
979 {
980 const double LOG2E = 1.44269504088896340735992468100189;
981 pari_sp ltop = avma;
982 GEN C = dbltor(2114.567);
983 double Mk, m, logbinom;
984 GEN tmp = mulrr(mppi(LOWDEFAULTPREC), sqrtr(utor(-D, LOWDEFAULTPREC)));
985 /* We treat this case separately since the table is not initialised when
986 * h = 1. This is the same as in the for loop below but with ak = 1. */
987 double log2Mk = dbllog2r(mpadd(mpexp(tmp), C));
988 double res = log2Mk;
989 ulong maxak = 1;
990 double log2Mh = log2Mk;
991
992 pari_sp btop = avma;
993 long k;
994 for (k = 2; k <= h; ++k) {
995 ulong ak = uel(qfinorms, k);
996 /* exp(tmp/a[k]) can overflow for even moderate discriminants, so we need
997 * to use t_REALs instead of doubles. Sutherland has a (more complicated)
998 * implementation in the classpoly package which should be consulted if
999 * this ever turns out to be a bottleneck.
1000 *
1001 * One idea to avoid t_REALs is the following: we have
1002 * log(e^x + C) - x <= log(2) ~ 0.69 for x >= log(C) ~ 0.44 and
1003 * the difference is basically zero for x slightly bigger than log(C).
1004 * Hence for large discriminants, we have x = \pi\sqrt{-D}/ak >> log(C)
1005 * and so we could approximate log(e^x + C) by x. */
1006 log2Mk = dbllog2r(mpadd(mpexp(divru(tmp, ak)), C));
1007 res += log2Mk;
1008 if (ak > maxak) { maxak = ak; log2Mh = log2Mk; }
1009 set_avma(btop);
1010 }
1011
1012 Mk = pow(2.0, log2Mh);
1013 m = floor((h + 1)/(Mk + 1.0));
1014 /* This line computes "log2(itos(binomialuu(h, m)))". The smallest
1015 * fundamental discriminant for which logbinom is not zero is
1016 * -1579751. */
1017 logbinom = (m > 0 && m < h)
1018 ? LOG2E * (logfac(h) - logfac(m) - logfac(h - m))
1019 : 0;
1020 set_avma(ltop);
1021 return res + logbinom - m * log2Mh + 2.0;
1022 }
1023
1024 INLINE long
distinct_inverses(const long f[],const long ef[],const long ei[],const long n[],const long o[],const long r[],long k,long L0,long i)1025 distinct_inverses(const long f[], const long ef[], const long ei[],
1026 const long n[], const long o[], const long r[], long k, long L0, long i)
1027 {
1028 pari_sp av = avma;
1029 long j, *e2, *e3;
1030
1031 if ( ! ef[i] || (L0 && ef[0])) return 0;
1032 for (j = i + 1; j < k; ++j)
1033 if (ef[j]) break;
1034 if (j < k) return 0;
1035
1036 e2 = new_chunk(k);
1037 evec_copy(e2, ef, i);
1038 e2[i] = o[i] - ef[i];
1039 for (j = i + 1; j < k; ++j) e2[j] = 0;
1040 evec_reduce(e2, n, r, k);
1041
1042 if (evec_equal(ef, e2, k)) return gc_long(av,0);
1043
1044 e3 = new_chunk(k);
1045 evec_inverse_o(e3, ef, n, o, r, k);
1046 if (evec_equal(e2, e3, k)) return gc_long(av,0);
1047
1048 if (f) {
1049 evec_compose(e3, f, ei, n, r, k);
1050 if (evec_equal(e2, e3, k)) return gc_long(av,0);
1051
1052 evec_inverse_o(e3, e3, n, o, r, k);
1053 if (evec_equal(e2, e3, k)) return gc_long(av,0);
1054 }
1055 return gc_long(av,1);
1056 }
1057
1058 INLINE long
next_prime_evec(long * qq,long f[],const long m[],long k,hashtable * tbl,long D,GEN DD,long u,long lvl,long ubound)1059 next_prime_evec(long *qq, long f[], const long m[], long k,
1060 hashtable *tbl, long D, GEN DD, long u, long lvl, long ubound)
1061 {
1062 pari_sp av = avma;
1063 hashentry *he;
1064 GEN P;
1065 long idx, q = *qq;
1066
1067 do q = unextprime(q + 1);
1068 while (!(u % q) || kross(D, q) == -1 || !(lvl % q) || !(D % (q * q)));
1069 if (q > ubound) return 0;
1070 *qq = q;
1071
1072 /* Get evec f corresponding to q */
1073 P = redimag(primeform_u(DD, q));
1074 he = hash_search(tbl, P);
1075 if (!he) pari_err_BUG("next_prime_evec");
1076 idx = itos((GEN) he->val);
1077 index_to_evec(f, idx, m, k);
1078 return gc_long(av,1);
1079 }
1080
1081 /* Return 1 on success, 0 on failure. */
1082 static int
orient_pcp(classgp_pcp_t G,long * ni,long D,long u,hashtable * tbl)1083 orient_pcp(classgp_pcp_t G, long *ni, long D, long u, hashtable *tbl)
1084 {
1085 pari_sp av = avma;
1086 /* 199 seems to suffice, but can be increased if necessary */
1087 enum { MAX_ORIENT_P = 199 };
1088 const long *L = G->L, *n = G->n, *r = G->r, *m = G->m, *o = G->o;
1089 long i, *ps = G->orient_p, *qs = G->orient_q, *reps = G->orient_reps;
1090 long *ef, *e, *ei, *f, k = G->k, lvl = modinv_level(G->inv);
1091 GEN DD = stoi(D);
1092
1093 memset(ps, 0, k * sizeof(long));
1094 memset(qs, 0, k * sizeof(long));
1095 memset(reps, 0, k * k * sizeof(long));
1096
1097 for (i = 0; i < k; ++i) { ps[i] = -1; if (o[i] > 2) break; }
1098 for (++i; i < k; ++i) ps[i] = (o[i] > 2) ? 0 : -1; /* ps[i] = -!(o[i] > 2); */
1099
1100 e = new_chunk(k);
1101 ei = new_chunk(k);
1102 f = new_chunk(k);
1103
1104 for (i = 0; i < k; ++i) {
1105 long p;
1106 if (ps[i]) continue;
1107 p = L[i];
1108 ef = &reps[i * k];
1109 while (!ps[i]) {
1110 if (!next_prime_evec(&p, ef, m, k, tbl, D, DD, u, lvl, MAX_ORIENT_P))
1111 break;
1112 evec_inverse_o(ei, ef, n, o, r, k);
1113 if (!distinct_inverses(NULL, ef, ei, n, o, r, k, G->L0, i)) continue;
1114 ps[i] = p;
1115 qs[i] = 1;
1116 }
1117 if (ps[i]) continue;
1118
1119 p = unextprime(L[i] + 1);
1120 while (!ps[i]) {
1121 long q;
1122
1123 if (!next_prime_evec(&p, e, m, k, tbl, D, DD, u, lvl, MAX_ORIENT_P))
1124 break;
1125 evec_inverse_o(ei, e, n, o, r, k);
1126
1127 q = L[i];
1128 while (!qs[i]) {
1129 if (!next_prime_evec(&q, f, m, k, tbl, D, DD, u, lvl, p - 1)) break;
1130 evec_compose(ef, e, f, n, r, k);
1131 if (!distinct_inverses(f, ef, ei, n, o, r, k, G->L0, i)) continue;
1132 ps[i] = p;
1133 qs[i] = q;
1134 }
1135 }
1136 if (!ps[i]) return 0;
1137 }
1138 if (ni) {
1139 GEN N = qfb_nform(D, *ni);
1140 hashentry *he = hash_search(tbl, N);
1141 if (!he) pari_err_BUG("orient_pcp");
1142 *ni = itos((GEN) he->val);
1143 }
1144 return gc_bool(av,1);
1145 }
1146
1147 /* We must avoid situations where L_i^{+/-2} = L_j^2 (or = L_0*L_j^2
1148 * if ell0 flag is set), with |L_i| = |L_j| = 4 (or have 4th powers in
1149 * <L0> but not 2nd powers in <L0>) and j < i */
1150 /* These cases cause problems when enumerating roots via gcds */
1151 /* returns the index of the first bad generator, or -1 if no bad
1152 * generators are found */
1153 static long
classgp_pcp_check_generators(const long * n,long * r,long k,long L0)1154 classgp_pcp_check_generators(const long *n, long *r, long k, long L0)
1155 {
1156 pari_sp av = avma;
1157 long *e1, i, i0, j, s;
1158 const long *ei;
1159
1160 s = !!L0;
1161 e1 = new_chunk(k);
1162
1163 for (i = s + 1; i < k; i++) {
1164 if (n[i] != 2) continue;
1165 ei = evec_ri(r, i);
1166 for (j = s; j < i; j++)
1167 if (ei[j]) break;
1168 if (j == i) continue;
1169 for (i0 = s; i0 < i; i0++) {
1170 if ((4 % n[i0])) continue;
1171 evec_clear(e1, k);
1172 e1[i0] = 4;
1173 evec_reduce(e1, n, r, k);
1174 for (j = s; j < i; j++)
1175 if (e1[j]) break;
1176 if (j < i) continue; /* L_i0^4 is not trivial or in <L_0> */
1177 evec_clear(e1, k);
1178 e1[i0] = 2;
1179 evec_reduce(e1, n, r, k); /* compute L_i0^2 */
1180 for (j = s; j < i; j++)
1181 if (e1[j] != ei[j]) break;
1182 if (j == i) return i;
1183 evec_inverse(e1, e1, n, r, k); /* compute L_i0^{-2} */
1184 for (j = s; j < i; j++)
1185 if (e1[j] != ei[j]) break;
1186 if (j == i) return i;
1187 }
1188 }
1189 return gc_long(av,-1);
1190 }
1191
1192 static void
pcp_alloc_and_set(classgp_pcp_t G,const long * L,const long * n,const long * r,long k)1193 pcp_alloc_and_set(
1194 classgp_pcp_t G, const long *L, const long *n, const long *r, long k)
1195 {
1196 /* classgp_pcp contains 6 arrays of length k (L, m, n, o, orient_p, orient_q),
1197 * one of length binom(k, 2) (r) and one of length k^2 (orient_reps) */
1198 long rlen = k * (k - 1) / 2, datalen = 6 * k + rlen + k * k;
1199 G->_data = newblock(datalen);
1200 G->L = G->_data;
1201 G->m = G->L + k;
1202 G->n = G->m + k;
1203 G->o = G->n + k;
1204 G->r = G->o + k;
1205 G->orient_p = G->r + rlen;
1206 G->orient_q = G->orient_p + k;
1207 G->orient_reps = G->orient_q + k;
1208 G->k = k;
1209
1210 evec_copy(G->L, L, k);
1211 evec_copy(G->n, n, k);
1212 evec_copy(G->r, r, rlen);
1213 evec_orders(G->o, n, r, k);
1214 evec_n_to_m(G->m, n, k);
1215 }
1216
1217 static void
classgp_pcp_clear(classgp_pcp_t G)1218 classgp_pcp_clear(classgp_pcp_t G)
1219 { if (G->_data) killblock(G->_data); }
1220
1221 /* This is Sutherland 2009, Algorithm 2.2 (p16). */
1222 static void
classgp_make_pcp(classgp_pcp_t G,double * height,long * ni,long h,long D,ulong u,long inv,long Lfilter,long orient)1223 classgp_make_pcp(
1224 classgp_pcp_t G, double *height, long *ni,
1225 long h, long D, ulong u, long inv, long Lfilter, long orient)
1226 {
1227 enum { MAX_GENS = 16, MAX_RLEN = MAX_GENS * (MAX_GENS - 1) / 2 };
1228 pari_sp av = avma, bv;
1229 long curr_p, h2, nelts, lvl = modinv_level(inv);
1230 GEN DD, ident, T, v;
1231 hashtable *tbl;
1232 long i, L1, L2;
1233 long k, L[MAX_GENS], n[MAX_GENS], r[MAX_RLEN];
1234
1235 memset(G, 0, sizeof *G);
1236
1237 G->D = D;
1238 G->h = h;
1239 G->inv = inv;
1240 G->L0 = (modinv_is_double_eta(inv) && modinv_ramified(D, inv))
1241 ? modinv_degree(NULL, NULL, inv) : 0;
1242 G->enum_cnt = h / (1 + !!G->L0);
1243 G->Lfilter = ulcm(Lfilter, lvl);
1244
1245 if (h == 1) {
1246 if (G->L0) pari_err_BUG("classgp_pcp");
1247 G->k = 0;
1248 G->_data = NULL;
1249 v = const_vecsmall(1, 1);
1250 *height = upper_bound_on_classpoly_coeffs(D, h, v);
1251 /* NB: No need to set *ni when h = 1 */
1252 set_avma(av); return;
1253 }
1254
1255 DD = stoi(D);
1256 bv = avma;
1257 while (1) {
1258 k = 0;
1259 /* Hash table has a QFI as a key and the (boxed) index of that QFI
1260 * in T as its value */
1261 tbl = hash_create(h, (ulong(*)(void*)) hash_GEN,
1262 (int(*)(void*,void*))&gequal, 1);
1263 ident = redimag(primeform_u(DD, 1));
1264 hash_insert(tbl, ident, gen_0);
1265
1266 T = vectrunc_init(h + 1);
1267 vectrunc_append(T, ident);
1268 nelts = 1;
1269 curr_p = 1;
1270
1271 while (nelts < h) {
1272 GEN gamma_i, beta;
1273 hashentry *e;
1274 long N = glength(T), Tlen = N, ri = 1;
1275
1276 if (k == MAX_GENS) pari_err_IMPL("classgp_pcp");
1277
1278 if (nelts == 1 && G->L0) {
1279 curr_p = G->L0;
1280 gamma_i = qfb_nform(D, curr_p);
1281 beta = redimag(gamma_i);
1282 if (equali1(gel(beta, 1)))
1283 {
1284 curr_p = 1;
1285 gamma_i = next_generator(DD, D, u, G->Lfilter, &beta, &curr_p);
1286 }
1287 } else
1288 gamma_i = next_generator(DD, D, u, G->Lfilter, &beta, &curr_p);
1289 while ((e = hash_search(tbl, beta)) == NULL) {
1290 long j;
1291 for (j = 1; j <= N; ++j) {
1292 GEN t = qficomp(beta, gel(T, j));
1293 vectrunc_append(T, t);
1294 hash_insert(tbl, t, stoi(Tlen++));
1295 }
1296 beta = qficomp(beta, gamma_i);
1297 ++ri;
1298 }
1299 if (ri > 1) {
1300 long j, si;
1301 L[k] = curr_p;
1302 n[k] = ri;
1303 nelts *= ri;
1304
1305 /* This is to reset the curr_p counter when we have G->L0 != 0
1306 * in the first position of L. */
1307 if (curr_p == G->L0) curr_p = 1;
1308
1309 N = 1;
1310 si = itos((GEN) e->val);
1311 for (j = 0; j < k; ++j) {
1312 evec_ri_mutate(r, k)[j] = (si / N) % n[j];
1313 N *= n[j];
1314 }
1315 ++k;
1316 }
1317 }
1318
1319 if ((i = classgp_pcp_check_generators(n, r, k, G->L0)) < 0) {
1320 pcp_alloc_and_set(G, L, n, r, k);
1321 if (!orient || orient_pcp(G, ni, D, u, tbl)) break;
1322 G->Lfilter *= G->L[0];
1323 classgp_pcp_clear(G);
1324 } else if (log2(G->Lfilter) + log2(L[i]) >= BITS_IN_LONG)
1325 pari_err_IMPL("classgp_pcp");
1326 else
1327 G->Lfilter *= L[i];
1328 set_avma(bv);
1329 }
1330
1331 v = cgetg(h + 1, t_VECSMALL);
1332 v[1] = 1;
1333 for (i = 2; i <= h; ++i) uel(v,i) = itou(gmael(T,i,1));
1334
1335 h2 = G->L0 ? h / 2 : h;
1336 *height = upper_bound_on_classpoly_coeffs(D, h2, v);
1337
1338 /* The norms of the last one or two generators. */
1339 L1 = L[k - 1];
1340 L2 = k > 1 ? L[k - 2] : 1;
1341 /* 4 * L1^2 * L2^2 must fit in a ulong */
1342 if (2 * (1 + log2(L1) + log2(L2)) >= BITS_IN_LONG)
1343 pari_err_IMPL("classgp_pcp");
1344
1345 if (G->L0 && (G->L[0] != G->L0 || G->o[0] != 2))
1346 pari_err_BUG("classgp_pcp");
1347
1348 set_avma(av); return;
1349 }
1350
1351 INLINE ulong
classno_wrapper(long D)1352 classno_wrapper(long D)
1353 {
1354 pari_sp av = avma;
1355 GEN G = quadclassunit0(stoi(D), 0, NULL, DEFAULTPREC);
1356 return gc_ulong(av, itou(abgrp_get_no(G)));
1357 }
1358
1359 /*
1360 * SECTION: Functions for calculating class polynomials.
1361 */
1362
1363 #define NSMALL_PRIMES 11
1364 static const long SMALL_PRIMES[11] = {
1365 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31
1366 };
1367
1368 static long
is_smooth_enough(ulong * factors,long v)1369 is_smooth_enough(ulong *factors, long v)
1370 {
1371 long i;
1372 *factors = 0;
1373 for (i = 0; i < NSMALL_PRIMES; ++i) {
1374 long p = SMALL_PRIMES[i];
1375 if (v % p == 0) *factors |= 1UL << i;
1376 while (v % p == 0) v /= p;
1377 if (v == 1) break;
1378 }
1379 return v == 1;
1380 }
1381
1382 /* Hurwitz class number of |D| assuming hclassno() and attached
1383 * conversion to double costs much more than unegisfundamental(). */
1384 INLINE double
hclassno_wrapper(long D,long h)1385 hclassno_wrapper(long D, long h)
1386 {
1387 /* TODO: Can probably calculate hurwitz faster using -D, factor(u)
1388 * and classno(D). */
1389 pari_sp av = avma;
1390 ulong abs_D = D < 0 ? -D : D;
1391 double hurwitz;
1392
1393 if (h && unegisfundamental(abs_D))
1394 hurwitz = (double) h;
1395 else
1396 hurwitz = rtodbl(gtofp(hclassno(utoi(abs_D)), DEFAULTPREC));
1397 set_avma(av); return hurwitz;
1398 }
1399
1400 /* This is Sutherland 2009, Algorithm 2.1 (p8); delta > 0 */
1401 static GEN
select_classpoly_prime_pool(double min_bits,double delta,classgp_pcp_t G)1402 select_classpoly_prime_pool(double min_bits, double delta, classgp_pcp_t G)
1403 { /* Sutherland defines V_MAX to be 1200 without saying why */
1404 const long V_MAX = 1200;
1405 pari_sp av = avma;
1406 double bits = 0.0, hurwitz, z;
1407 ulong t_size_lim;
1408 long ires, D = G->D, inv = G->inv;
1409 GEN res, t_min; /* t_min[v] = lower bound for the t we look at for that v */
1410
1411 hurwitz = hclassno_wrapper(D, G->h);
1412
1413 res = cgetg(128+1, t_VEC);
1414 ires = 1;
1415 /* Initialise t_min to be all 2's. This avoids trace 0 and trace 1 curves */
1416 t_min = const_vecsmall(V_MAX-1, 2);
1417
1418 /* maximum possible trace = sqrt(2^BIL - D) */
1419 t_size_lim = 2.0 * sqrt((double)((1UL << (BITS_IN_LONG - 2)) - (((ulong)-D) >> 2)));
1420
1421 for (z = -D / (2.0 * hurwitz); ; z *= delta + 1.0)
1422 { /* v_bound_aux = -4 z H(-D) */
1423 double v_bound_aux = -4.0 * z * hurwitz;
1424 ulong v;
1425 dbg_printf(1)("z = %.2f\n", z);
1426 for (v = 1; v < V_MAX; v++)
1427 {
1428 ulong p, t, t_max, vfactors, m_vsqr_D = v * v * (ulong)(-D);
1429 /* hurwitz_ratio_bound = 11 * log(log(v + 4))^2 */
1430 double hurwitz_ratio_bound = log(log(v + 4.0)), max_p, H;
1431 long ires0;
1432 hurwitz_ratio_bound *= 11.0 * hurwitz_ratio_bound;
1433
1434 if (v >= v_bound_aux * hurwitz_ratio_bound / D) break;
1435 if (!is_smooth_enough(&vfactors, v)) continue;
1436 H = hclassno_wrapper(m_vsqr_D, 0);
1437
1438 /* t <= 2 sqrt(p) and p <= z H(-v^2 D) and
1439 * H(-v^2 D) < vH(-D) (11 log(log(v + 4))^2)
1440 * This last term is v * hurwitz * hurwitz_ratio_bound. */
1441 max_p = z * v * hurwitz * hurwitz_ratio_bound;
1442 t_max = 2.0 * sqrt(mindd((1UL<<(BITS_IN_LONG-2)) - (m_vsqr_D>>2), max_p));
1443 t = t_min[v]; if ((t & 1) != (m_vsqr_D & 1)) t++;
1444 p = (t * t + m_vsqr_D) >> 2;
1445 ires0 = ires;
1446 for (; t <= t_max; p += t+1, t += 2) /* 4p = t^2 - v^2*D */
1447 if (modinv_good_prime(inv,p) && uisprime(p))
1448 {
1449 if (ires == lg(res)) res = vec_lengthen(res, lg(res) << 1);
1450 gel(res, ires++) = mkvecsmall5(p, t, v, (long)(p / H), vfactors);
1451 bits += log2(p);
1452 }
1453 t_min[v] = t;
1454
1455 if (ires - ires0) {
1456 dbg_printf(2)(" Found %lu primes for v = %lu.\n", ires - ires0, v);
1457 }
1458 if (bits > min_bits) {
1459 dbg_printf(1)("Found %ld primes; total size %.2f bits.\n", ires-1,bits);
1460 setlg(res, ires); return gerepilecopy(av, res);
1461 }
1462 }
1463 if (uel(t_min,1) >= t_size_lim) {
1464 /* exhausted all solutions that fit in ulong */
1465 char *err = stack_sprintf("class polynomial of discriminant %ld", D);
1466 pari_err(e_ARCH, err);
1467 }
1468 }
1469 }
1470
1471 static int
primecmp(void * data,GEN v1,GEN v2)1472 primecmp(void *data, GEN v1, GEN v2)
1473 { (void)data; return cmpss(v1[4], v2[4]); }
1474
1475 static long
height_margin(long inv,long D)1476 height_margin(long inv, long D)
1477 {
1478 (void)D;
1479 /* NB: avs just uses a height margin of 256 for everyone and everything. */
1480 if (inv == INV_F) return 64; /* checked for discriminants up to -350000 */
1481 if (inv == INV_G2) return 5;
1482 if (inv != INV_J) return 256; /* TODO: This should be made more accurate */
1483 return 0;
1484 }
1485
1486 static GEN
select_classpoly_primes(ulong * vfactors,ulong * biggest_v,long k,double delta,classgp_pcp_t G,double height)1487 select_classpoly_primes(
1488 ulong *vfactors, ulong *biggest_v,
1489 long k, double delta, classgp_pcp_t G, double height)
1490 {
1491 pari_sp av = avma;
1492 long i, s, D = G->D, inv = G->inv;
1493 ulong biggest_p;
1494 double prime_bits, min_prime_bits, b;
1495 GEN prime_pool;
1496
1497 if (k < 2) pari_err_BUG("select_suitable_primes");
1498
1499 s = modinv_height_factor(inv);
1500 b = height / s + height_margin(inv, D);
1501 dbg_printf(1)("adjusted height = %.2f\n", b);
1502 min_prime_bits = k * b;
1503
1504 prime_pool = select_classpoly_prime_pool(min_prime_bits, delta, G);
1505
1506 /* FIXME: Apply torsion constraints */
1507 /* FIXME: Rank elts of res according to cost/benefit ratio */
1508 gen_sort_inplace(prime_pool, NULL, primecmp, NULL);
1509
1510 prime_bits = 0.0;
1511 biggest_p = gel(prime_pool, 1)[1];
1512 *biggest_v = gel(prime_pool, 1)[3];
1513 *vfactors = 0;
1514 for (i = 1; i < lg(prime_pool); ++i) {
1515 ulong p = gel(prime_pool, i)[1];
1516 ulong v = gel(prime_pool, i)[3];
1517 prime_bits += log2(p);
1518 *vfactors |= gel(prime_pool, i)[5];
1519 if (p > biggest_p) biggest_p = p;
1520 if (v > *biggest_v) *biggest_v = v;
1521 if (prime_bits > b) break;
1522 }
1523 dbg_printf(1)("Selected %ld primes; largest is %lu ~ 2^%.2f\n",
1524 i, biggest_p, log2(biggest_p));
1525 return gerepilecopy(av, vecslice0(prime_pool, 1, i));
1526 }
1527
1528 /* This is Sutherland 2009 Algorithm 1.2. */
1529 static long
oneroot_of_classpoly(ulong * j_endo,int * endo_cert,ulong j,norm_eqn_t ne,GEN jdb)1530 oneroot_of_classpoly(
1531 ulong *j_endo, int *endo_cert, ulong j, norm_eqn_t ne, GEN jdb)
1532 {
1533 pari_sp av = avma;
1534 long nfactors, L_bound, i;
1535 ulong p = ne->p, pi = ne->pi;
1536 GEN factw, factors, u_levels, vdepths;
1537
1538 if (j == 0 || j == 1728 % p) pari_err_BUG("oneroot_of_classpoly");
1539
1540 *endo_cert = 1;
1541 if (ne->u * ne->v == 1) { *j_endo = j; return 1; }
1542
1543 /* TODO: Precalculate all this data further up */
1544 factw = factoru(ne->u * ne->v);
1545 factors = gel(factw, 1);
1546 nfactors = lg(factors) - 1;
1547 u_levels = cgetg(nfactors + 1, t_VECSMALL);
1548 for (i = 1; i <= nfactors; ++i)
1549 u_levels[i] = z_lval(ne->u, gel(factw, 1)[i]);
1550 vdepths = gel(factw, 2);
1551
1552 /* FIXME: This should be bigger */
1553 L_bound = maxdd(log((double) -ne->D), (double)ne->v);
1554
1555 /* Iterate over the primes L dividing w */
1556 for (i = 1; i <= nfactors; ++i) {
1557 pari_sp bv = avma;
1558 GEN phi;
1559 long jlvl, lvl_diff, depth = vdepths[i];
1560 long L = factors[i];
1561 if (L > L_bound) { *endo_cert = 0; break; }
1562
1563 phi = polmodular_db_getp(jdb, L, p);
1564
1565 /* TODO: See if I can reuse paths created in j_level_in_volcano()
1566 * later in {ascend,descend}_volcano(), perhaps by combining the
1567 * functions into one "adjust_level" function. */
1568 jlvl = j_level_in_volcano(phi, j, p, pi, L, depth);
1569 lvl_diff = u_levels[i] - jlvl;
1570
1571 if (lvl_diff < 0)
1572 /* j's level is less than v(u) so we must ascend */
1573 j = ascend_volcano(phi, j, p, pi, jlvl, L, depth, -lvl_diff);
1574 else if (lvl_diff > 0)
1575 /* otherwise j's level is greater than v(u) so we descend */
1576 j = descend_volcano(phi, j, p, pi, jlvl, L, depth, lvl_diff);
1577 set_avma(bv);
1578 }
1579 set_avma(av);
1580 /* At this point the probability that j has the wrong endomorphism
1581 * ring is about \sum_{p|u_compl} 1/p (and u_compl must be bigger
1582 * than L_bound, so pretty big), so just return it and rely on
1583 * detection code in enum_j_with_endo_ring(). Detection is that we
1584 * hit a previously found j-invariant earlier than expected. OR, we
1585 * evaluate class polynomials of the suborders at j and if any are
1586 * zero then j must be chosen again. */
1587 *j_endo = j;
1588 return j != 0 && j != 1728 % p;
1589 }
1590
1591 INLINE long
vecsmall_isin_skip(GEN v,long x,long k)1592 vecsmall_isin_skip(GEN v, long x, long k)
1593 {
1594 long i, l = lg(v);
1595 for (i = k; i < l; ++i)
1596 if (v[i] == x) return i;
1597 return 0;
1598 }
1599
1600 INLINE ulong
select_twisting_param(ulong p)1601 select_twisting_param(ulong p)
1602 {
1603 ulong T;
1604 do T = random_Fl(p); while (krouu(T, p) != -1);
1605 return T;
1606 }
1607
1608 INLINE void
setup_norm_eqn(norm_eqn_t ne,long D,long u,GEN norm_eqn)1609 setup_norm_eqn(norm_eqn_t ne, long D, long u, GEN norm_eqn)
1610 {
1611 ne->D = D;
1612 ne->u = u;
1613 ne->t = norm_eqn[2];
1614 ne->v = norm_eqn[3];
1615 ne->p = (ulong) norm_eqn[1];
1616 ne->pi = get_Fl_red(ne->p);
1617 ne->s2 = Fl_2gener_pre(ne->p, ne->pi);
1618 ne->T = select_twisting_param(ne->p);
1619 }
1620
1621 INLINE ulong
Flv_powsum_pre(GEN v,ulong n,ulong p,ulong pi)1622 Flv_powsum_pre(GEN v, ulong n, ulong p, ulong pi)
1623 {
1624 long i, l = lg(v);
1625 ulong psum = 0;
1626 for (i = 1; i < l; ++i)
1627 psum = Fl_add(psum, Fl_powu_pre(uel(v,i), n, p, pi), p);
1628 return psum;
1629 }
1630
1631 INLINE int
modinv_has_sign_ambiguity(long inv)1632 modinv_has_sign_ambiguity(long inv)
1633 {
1634 switch (inv) {
1635 case INV_F:
1636 case INV_F3:
1637 case INV_W2W3E2:
1638 case INV_W2W7E2:
1639 case INV_W2W3:
1640 case INV_W2W5:
1641 case INV_W2W7:
1642 case INV_W3W3:
1643 case INV_W2W13:
1644 case INV_W3W7: return 1;
1645 }
1646 return 0;
1647 }
1648
1649 INLINE int
modinv_units(int inv)1650 modinv_units(int inv)
1651 { return modinv_is_double_eta(inv) || modinv_is_Weber(inv); }
1652
1653 INLINE int
adjust_signs(GEN js,ulong p,ulong pi,long inv,GEN T,long e)1654 adjust_signs(GEN js, ulong p, ulong pi, long inv, GEN T, long e)
1655 {
1656 long negate = 0;
1657 long h = lg(js) - 1;
1658 if ((h & 1) && modinv_units(inv)) {
1659 ulong prod = Flv_prod_pre(js, p, pi);
1660 if (prod != p - 1) {
1661 if (prod != 1) pari_err_BUG("adjust_signs: constant term is not +/-1");
1662 negate = 1;
1663 }
1664 } else {
1665 ulong tp, t;
1666 tp = umodiu(T, p);
1667 t = Flv_powsum_pre(js, e, p, pi);
1668 if (t == 0) return 0;
1669 if (t != tp) {
1670 if (Fl_neg(t, p) != tp) pari_err_BUG("adjust_signs: incorrect trace");
1671 negate = 1;
1672 }
1673 }
1674 if (negate) Flv_neg_inplace(js, p);
1675 return 1;
1676 }
1677
1678 static ulong
find_jinv(long * trace_tries,long * endo_tries,int * cert,norm_eqn_t ne,long inv,long rho_inv,GEN jdb)1679 find_jinv(
1680 long *trace_tries, long *endo_tries, int *cert,
1681 norm_eqn_t ne, long inv, long rho_inv, GEN jdb)
1682 {
1683 long found, ok = 1;
1684 ulong j, r;
1685 do {
1686 do {
1687 long tries;
1688 ulong j_t = 0;
1689 /* TODO: Set batch size according to expected number of tries and
1690 * experimental cost/benefit analysis. */
1691 tries = find_j_inv_with_given_trace(&j_t, ne, rho_inv, 0);
1692 if (j_t == 0)
1693 pari_err_BUG("polclass0: Couldn't find j-invariant with given trace.");
1694 dbg_printf(2)(" j-invariant %ld has trace +/-%ld (%ld tries, 1/rho = %ld)\n",
1695 j_t, ne->t, tries, rho_inv);
1696 *trace_tries += tries;
1697
1698 found = oneroot_of_classpoly(&j, cert, j_t, ne, jdb);
1699 ++*endo_tries;
1700 } while (!found);
1701
1702 if (modinv_is_double_eta(inv))
1703 ok = modfn_unambiguous_root(&r, inv, j, ne, jdb);
1704 else
1705 r = modfn_root(j, ne, inv);
1706 } while (!ok);
1707 return r;
1708 }
1709
1710 static GEN
polclass_roots_modp(long * n_trace_curves,norm_eqn_t ne,long rho_inv,classgp_pcp_t G,GEN db)1711 polclass_roots_modp(
1712 long *n_trace_curves,
1713 norm_eqn_t ne, long rho_inv, classgp_pcp_t G, GEN db)
1714 {
1715 pari_sp av = avma;
1716 ulong j = 0;
1717 long inv = G->inv, endo_tries = 0;
1718 int endo_cert;
1719 GEN res, jdb, fdb;
1720
1721 jdb = polmodular_db_for_inv(db, INV_J);
1722 fdb = polmodular_db_for_inv(db, inv);
1723
1724 dbg_printf(2)("p = %ld, t = %ld, v = %ld\n", ne->p, ne->t, ne->v);
1725
1726 do {
1727 j = find_jinv(n_trace_curves, &endo_tries, &endo_cert, ne, inv, rho_inv, jdb);
1728
1729 res = enum_roots(j, ne, fdb, G);
1730 if ( ! res && endo_cert) pari_err_BUG("polclass_roots_modp");
1731 if (res && ! endo_cert && vecsmall_isin_skip(res, res[1], 2))
1732 {
1733 set_avma(av);
1734 res = NULL;
1735 }
1736 } while (!res);
1737
1738 dbg_printf(2)(" j-invariant %ld has correct endomorphism ring "
1739 "(%ld tries)\n", j, endo_tries);
1740 dbg_printf(4)(" all such j-invariants: %Ps\n", res);
1741 return gerepileupto(av, res);
1742 }
1743
1744 INLINE int
modinv_inverted_involution(long inv)1745 modinv_inverted_involution(long inv)
1746 { return modinv_is_double_eta(inv); }
1747
1748 INLINE int
modinv_negated_involution(long inv)1749 modinv_negated_involution(long inv)
1750 { /* determined by trial and error */
1751 return inv == INV_F || inv == INV_W3W5 || inv == INV_W3W7
1752 || inv == INV_W3W3 || inv == INV_W5W7;
1753 }
1754
1755 /* Return true iff Phi_L(j0, j1) = 0. */
1756 INLINE long
verify_edge(ulong j0,ulong j1,ulong p,ulong pi,long L,GEN fdb)1757 verify_edge(ulong j0, ulong j1, ulong p, ulong pi, long L, GEN fdb)
1758 {
1759 pari_sp av = avma;
1760 GEN phi = polmodular_db_getp(fdb, L, p);
1761 GEN f = Flm_Fl_polmodular_evalx(phi, L, j1, p, pi);
1762 return gc_long(av, Flx_eval_pre(f, j0, p, pi) == 0);
1763 }
1764
1765 INLINE long
verify_2path(ulong j1,ulong j2,ulong p,ulong pi,long L1,long L2,GEN fdb)1766 verify_2path(
1767 ulong j1, ulong j2, ulong p, ulong pi, long L1, long L2, GEN fdb)
1768 {
1769 pari_sp av = avma;
1770 GEN phi1 = polmodular_db_getp(fdb, L1, p);
1771 GEN phi2 = polmodular_db_getp(fdb, L2, p);
1772 GEN f = Flm_Fl_polmodular_evalx(phi1, L1, j1, p, pi);
1773 GEN g = Flm_Fl_polmodular_evalx(phi2, L2, j2, p, pi);
1774 GEN d = Flx_gcd(f, g, p);
1775 long n = degpol(d);
1776 if (n >= 2) n = Flx_nbroots(d, p);
1777 return gc_long(av, n);
1778 }
1779
1780 static long
oriented_n_action(const long * ni,classgp_pcp_t G,GEN v,ulong p,ulong pi,GEN fdb)1781 oriented_n_action(
1782 const long *ni, classgp_pcp_t G, GEN v, ulong p, ulong pi, GEN fdb)
1783 {
1784 pari_sp av = avma;
1785 long i, j, k = G->k;
1786 long nr = k * (k - 1) / 2;
1787 const long *n = G->n, *m = G->m, *o = G->o, *r = G->r,
1788 *ps = G->orient_p, *qs = G->orient_q, *reps = G->orient_reps;
1789 long *signs = new_chunk(k);
1790 long *e = new_chunk(k);
1791 long *rels = new_chunk(nr);
1792
1793 evec_copy(rels, r, nr);
1794
1795 for (i = 0; i < k; ++i) {
1796 /* If generator doesn't require orientation, continue; power rels already
1797 * copied to *rels in initialisation */
1798 if (ps[i] <= 0) { signs[i] = 1; continue; }
1799 /* Get rep of orientation element and express it in terms of the
1800 * (partially) oriented presentation */
1801 for (j = 0; j < i; ++j) {
1802 long t = reps[i * k + j];
1803 e[j] = (signs[j] < 0 ? o[j] - t : t);
1804 }
1805 e[j] = reps[i * k + j];
1806 for (++j; j < k; ++j) e[j] = 0;
1807 evec_reduce(e, n, rels, k);
1808 j = evec_to_index(e, m, k);
1809
1810 /* FIXME: These calls to verify_edge recalculate powers of v[0]
1811 * and v[j] over and over again, they also reduce Phi_{ps[i]} modulo p over
1812 * and over again. Need to cache these things! */
1813 if (qs[i] > 1)
1814 signs[i] =
1815 (verify_2path(uel(v,1), uel(v,j+1), p, pi, ps[i], qs[i], fdb) ? 1 : -1);
1816 else
1817 /* Verify ps[i]-edge to orient ith generator */
1818 signs[i] =
1819 (verify_edge(uel(v,1), uel(v,j+1), p, pi, ps[i], fdb) ? 1 : -1);
1820 /* Update power relation */
1821 for (j = 0; j < i; ++j) {
1822 long t = evec_ri(r, i)[j];
1823 e[j] = (signs[i] * signs[j] < 0 ? o[j] - t : t);
1824 }
1825 while (j < k) e[j++] = 0;
1826 evec_reduce(e, n, rels, k);
1827 for (j = 0; j < i; ++j) evec_ri_mutate(rels, i)[j] = e[j];
1828 /* TODO: This is a sanity check, can be removed if everything is working */
1829 for (j = 0; j <= i; ++j) {
1830 long t = reps[i * k + j];
1831 e[j] = (signs[j] < 0 ? o[j] - t : t);
1832 }
1833 while (j < k) e[j++] = 0;
1834 evec_reduce(e, n, rels, k);
1835 j = evec_to_index(e, m, k);
1836 if (qs[i] > 1) {
1837 if (!verify_2path(uel(v,1), uel(v, j+1), p, pi, ps[i], qs[i], fdb))
1838 pari_err_BUG("oriented_n_action");
1839 } else {
1840 if (!verify_edge(uel(v,1), uel(v, j+1), p, pi, ps[i], fdb))
1841 pari_err_BUG("oriented_n_action");
1842 }
1843 }
1844
1845 /* Orient representation of [N] relative to the torsor <signs, rels> */
1846 for (i = 0; i < k; ++i) e[i] = (signs[i] < 0 ? o[i] - ni[i] : ni[i]);
1847 evec_reduce(e, n, rels, k);
1848 return gc_long(av, evec_to_index(e,m,k));
1849 }
1850
1851 /* F = double_eta_raw(inv) */
1852 INLINE void
adjust_orientation(GEN F,long inv,GEN v,long e,ulong p,ulong pi)1853 adjust_orientation(GEN F, long inv, GEN v, long e, ulong p, ulong pi)
1854 {
1855 ulong j0 = uel(v, 1), je = uel(v, e);
1856
1857 if (!modinv_j_from_2double_eta(F, inv, j0, je, p, pi)) {
1858 if (modinv_inverted_involution(inv)) Flv_inv_pre_inplace(v, p, pi);
1859 if (modinv_negated_involution(inv)) Flv_neg_inplace(v, p);
1860 }
1861 }
1862
1863 static void
polclass_psum(GEN * psum,long * d,GEN roots,GEN primes,GEN pilist,ulong h,long inv)1864 polclass_psum(
1865 GEN *psum, long *d, GEN roots, GEN primes, GEN pilist, ulong h, long inv)
1866 {
1867 /* Number of consecutive CRT stabilisations before we assume we have
1868 * the correct answer. */
1869 enum { MIN_STAB_CNT = 3 };
1870 pari_sp av = avma, btop;
1871 GEN ps, psum_sqr, P;
1872 long i, e, stabcnt, nprimes = lg(primes) - 1;
1873
1874 if ((h & 1) && modinv_units(inv)) { *psum = gen_1; *d = 0; return; }
1875 e = -1;
1876 ps = cgetg(nprimes+1, t_VECSMALL);
1877 do {
1878 e += 2;
1879 for (i = 1; i <= nprimes; ++i)
1880 {
1881 GEN roots_modp = gel(roots, i);
1882 ulong p = uel(primes, i), pi = uel(pilist, i);
1883 uel(ps, i) = Flv_powsum_pre(roots_modp, e, p, pi);
1884 }
1885 btop = avma;
1886 psum_sqr = Z_init_CRT(0, 1);
1887 P = gen_1;
1888 for (i = 1, stabcnt = 0; stabcnt < MIN_STAB_CNT && i <= nprimes; ++i)
1889 {
1890 ulong p = uel(primes, i), pi = uel(pilist, i);
1891 ulong ps2 = Fl_sqr_pre(uel(ps, i), p, pi);
1892 ulong stab = Z_incremental_CRT(&psum_sqr, ps2, &P, p);
1893 /* stabcnt = stab * (stabcnt + 1) */
1894 if (stab) ++stabcnt; else stabcnt = 0;
1895 if (gc_needed(av, 2)) gerepileall(btop, 2, &psum_sqr, &P);
1896 }
1897 if (stabcnt == 0 && nprimes >= MIN_STAB_CNT)
1898 pari_err_BUG("polclass_psum");
1899 } while (!signe(psum_sqr));
1900
1901 if ( ! Z_issquareall(psum_sqr, psum)) pari_err_BUG("polclass_psum");
1902
1903 dbg_printf(1)("Classpoly power sum (e = %ld) is %Ps; found with %.2f%% of the primes\n",
1904 e, *psum, 100 * (i - 1) / (double) nprimes);
1905 *psum = gerepileupto(av, *psum);
1906 *d = e;
1907 }
1908
1909 static GEN
polclass_small_disc(long D,long inv,long xvar)1910 polclass_small_disc(long D, long inv, long xvar)
1911 {
1912 if (D == -3) return pol_x(xvar);
1913 if (D == -4) {
1914 switch (inv) {
1915 case INV_J: return deg1pol(gen_1, stoi(-1728), xvar);
1916 case INV_G2:return deg1pol(gen_1, stoi(-12), xvar);
1917 default: /* no other invariants for which we can calculate H_{-4}(X) */
1918 pari_err_BUG("polclass_small_disc");
1919 }
1920 }
1921 return NULL;
1922 }
1923
1924 GEN
polclass0(long D,long inv,long xvar,GEN * db)1925 polclass0(long D, long inv, long xvar, GEN *db)
1926 {
1927 pari_sp av = avma;
1928 GEN primes;
1929 long n_curves_tested = 0;
1930 long nprimes, s, i, j, del, ni, orient;
1931 GEN P, H, plist, pilist;
1932 ulong u, L, maxL, vfactors, biggest_v;
1933 long h, p1, p2, filter = 1;
1934 classgp_pcp_t G;
1935 double height;
1936 static const long k = 2;
1937 static const double delta = 0.5;
1938
1939 if (D >= -4) return polclass_small_disc(D, inv, xvar);
1940
1941 (void) corediscs(D, &u);
1942 h = classno_wrapper(D);
1943
1944 dbg_printf(1)("D = %ld, conductor = %ld, inv = %ld\n", D, u, inv);
1945
1946 ni = modinv_degree(&p1, &p2, inv);
1947 orient = modinv_is_double_eta(inv) && kross(D, p1) && kross(D, p2);
1948
1949 classgp_make_pcp(G, &height, &ni, h, D, u, inv, filter, orient);
1950 primes = select_classpoly_primes(&vfactors, &biggest_v, k, delta, G, height);
1951
1952 /* Prepopulate *db with all the modpolys we might need */
1953 /* TODO: Clean this up; in particular, note that u is factored later on. */
1954 /* This comes from L_bound in oneroot_of_classpoly() above */
1955 maxL = maxdd(log((double) -D), (double)biggest_v);
1956 if (u > 1) {
1957 for (L = 2; L <= maxL; L = unextprime(L + 1))
1958 if (!(u % L)) polmodular_db_add_level(db, L, INV_J);
1959 }
1960 for (i = 0; vfactors; ++i) {
1961 if (vfactors & 1UL)
1962 polmodular_db_add_level(db, SMALL_PRIMES[i], INV_J);
1963 vfactors >>= 1;
1964 }
1965 if (p1 > 1) polmodular_db_add_level(db, p1, INV_J);
1966 if (p2 > 1) polmodular_db_add_level(db, p2, INV_J);
1967 s = !!G->L0;
1968 polmodular_db_add_levels(db, G->L + s, G->k - s, inv);
1969 if (orient) {
1970 for (i = 0; i < G->k; ++i)
1971 {
1972 if (G->orient_p[i] > 1) polmodular_db_add_level(db, G->orient_p[i], inv);
1973 if (G->orient_q[i] > 1) polmodular_db_add_level(db, G->orient_q[i], inv);
1974 }
1975 }
1976 nprimes = lg(primes) - 1;
1977 H = cgetg(nprimes + 1, t_VEC);
1978 plist = cgetg(nprimes + 1, t_VECSMALL);
1979 pilist = cgetg(nprimes + 1, t_VECSMALL);
1980 for (i = 1; i <= nprimes; ++i) {
1981 long rho_inv = gel(primes, i)[4];
1982 norm_eqn_t ne;
1983 setup_norm_eqn(ne, D, u, gel(primes, i));
1984
1985 gel(H, i) = polclass_roots_modp(&n_curves_tested, ne, rho_inv, G, *db);
1986 uel(plist, i) = ne->p;
1987 uel(pilist, i) = ne->pi;
1988 if (DEBUGLEVEL>2 && (i & 3L)==0) err_printf(" %ld%%", i*100/nprimes);
1989 }
1990 dbg_printf(0)("\n");
1991
1992 if (orient) {
1993 GEN nvec = new_chunk(G->k);
1994 GEN fdb = polmodular_db_for_inv(*db, inv);
1995 GEN F = double_eta_raw(inv);
1996 index_to_evec((long *)nvec, ni, G->m, G->k);
1997 for (i = 1; i <= nprimes; ++i) {
1998 GEN v = gel(H, i);
1999 ulong p = uel(plist, i), pi = uel(pilist, i);
2000 long oni = oriented_n_action(nvec, G, v, p, pi, fdb);
2001 adjust_orientation(F, inv, v, oni + 1, p, pi);
2002 }
2003 }
2004
2005 if (modinv_has_sign_ambiguity(inv)) {
2006 GEN psum;
2007 long e;
2008 polclass_psum(&psum, &e, H, plist, pilist, h, inv);
2009 for (i = 1; i <= nprimes; ++i) {
2010 GEN v = gel(H, i);
2011 ulong p = uel(plist, i), pi = uel(pilist, i);
2012 if (!adjust_signs(v, p, pi, inv, psum, e))
2013 uel(plist, i) = 0;
2014 }
2015 }
2016
2017 for (i = 1, j = 1, del = 0; i <= nprimes; ++i) {
2018 GEN v = gel(H, i), pol;
2019 ulong p = uel(plist, i);
2020 if (!p) { del++; continue; }
2021 pol = Flv_roots_to_pol(v, p, xvar);
2022 uel(plist, j) = p;
2023 gel(H, j++) = Flx_to_Flv(pol, lg(pol) - 2);
2024 }
2025 setlg(H,nprimes+1-del);
2026 setlg(plist,nprimes+1-del);
2027 classgp_pcp_clear(G);
2028
2029 dbg_printf(1)("Total number of curves tested: %ld\n", n_curves_tested);
2030 H = ncV_chinese_center(H, plist, &P);
2031 dbg_printf(1)("Result height: %.2f\n",
2032 dbllog2r(itor(gsupnorm(H, DEFAULTPREC), DEFAULTPREC)));
2033 return gerepilecopy(av, RgV_to_RgX(H, xvar));
2034 }
2035
2036 void
check_modinv(long inv)2037 check_modinv(long inv)
2038 {
2039 switch (inv) {
2040 case INV_J:
2041 case INV_F:
2042 case INV_F2:
2043 case INV_F3:
2044 case INV_F4:
2045 case INV_G2:
2046 case INV_W2W3:
2047 case INV_F8:
2048 case INV_W3W3:
2049 case INV_W2W5:
2050 case INV_W2W7:
2051 case INV_W3W5:
2052 case INV_W3W7:
2053 case INV_W2W3E2:
2054 case INV_W2W5E2:
2055 case INV_W2W13:
2056 case INV_W2W7E2:
2057 case INV_W3W3E2:
2058 case INV_W5W7:
2059 case INV_W3W13:
2060 break;
2061 default:
2062 pari_err_DOMAIN("polmodular", "inv", "invalid invariant", stoi(inv), gen_0);
2063 }
2064 }
2065
2066 GEN
polclass(GEN DD,long inv,long xvar)2067 polclass(GEN DD, long inv, long xvar)
2068 {
2069 GEN db, H;
2070 long dummy, D;
2071
2072 if (xvar < 0) xvar = 0;
2073 check_quaddisc_imag(DD, &dummy, "polclass");
2074 check_modinv(inv);
2075
2076 D = itos(DD);
2077 if (!modinv_good_disc(inv, D))
2078 pari_err_DOMAIN("polclass", "D", "incompatible with given invariant", stoi(inv), DD);
2079
2080 db = polmodular_db_init(inv);
2081 H = polclass0(D, inv, xvar, &db);
2082 gunclone_deep(db); return H;
2083 }
2084