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