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 9351 : pcp_get_L(GEN G) { return gmael(G,1,1)+1; }
24 : static GEN
25 6495 : pcp_get_n(GEN G) { return gmael(G,1,2)+1; }
26 : static GEN
27 6768 : pcp_get_m(GEN G) { return gmael(G,1,4)+1; }
28 : static GEN
29 6495 : pcp_get_o(GEN G) { return gmael(G,1,3)+1; }
30 : static GEN
31 6495 : pcp_get_r(GEN G) { return gmael(G,1,5)+1; }
32 : static GEN
33 6768 : pcp_get_orient_p(GEN G) { return gmael(G,3,1)+1; }
34 : static GEN
35 6768 : pcp_get_orient_q(GEN G) { return gmael(G,3,2)+1; }
36 : static GEN
37 6495 : pcp_get_orient_reps(GEN G) { return gmael(G,3,3)+1; }
38 : static long
39 9351 : pcp_get_L0(GEN G) { return mael(G,2,1); }
40 : static long
41 11384 : pcp_get_k(GEN G) { return mael(G,2,2); }
42 : static long
43 194617 : pcp_get_h(GEN G) { return mael(G,4,1); }
44 : static long
45 149552 : pcp_get_inv(GEN G) { return mael(G,4,2); }
46 : static long
47 149391 : pcp_get_D(GEN G) { return mael(G,4,3); }
48 : static long
49 194617 : pcp_get_D0(GEN G) { return mael(G,4,4); }
50 : static long
51 139624 : pcp_get_u(GEN G) { return mael(G,4,5); }
52 : static GEN
53 4889 : 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 140884 : hasse_bounds(long *low, long *high, long p)
60 140884 : { 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 3036 : famatsmall_divexact(GEN a, GEN b)
65 : {
66 3036 : GEN a1 = gel(a,1), a2 = gel(a,2), c1, c2;
67 3036 : GEN b1 = gel(b,1), b2 = gel(b,2);
68 3036 : long i, j, k, la = lg(a1);
69 3036 : c1 = cgetg(la, t_VECSMALL);
70 3036 : c2 = cgetg(la, t_VECSMALL);
71 10969 : for (i = j = k = 1; j < la; j++)
72 : {
73 7933 : c1[k] = a1[j];
74 7933 : c2[k] = a2[j];
75 7933 : if (a1[j] == b1[i]) { c2[k] -= b2[i++]; if (!c2[k]) continue; }
76 4950 : k++;
77 : }
78 3036 : setlg(c1, k);
79 3036 : 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 168989 : 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 168989 : pari_sp ltop = avma, av;
91 168989 : ulong a4t, a6t, p = ne->p, pi = ne->pi, T = ne->T;
92 : long m0, m1;
93 168989 : int swapped = 0, first = 1;
94 :
95 168989 : if (p <= 11) {
96 433 : long card = (long)p + 1 - Fl_elltrace(a4, a6, p);
97 433 : return card == N0 || card == N1;
98 : }
99 168556 : m0 = m1 = 1;
100 168556 : for (av = avma;;)
101 132116 : {
102 : long a1, x, n_s;
103 300672 : GEN Q = random_Flj_pre(a4, a6, p, pi);
104 300670 : if (m0 == 1)
105 297535 : n_s = Flj_order_ufact(Q, N0, n0, a4, p, pi);
106 3135 : else if (N0 % m0) n_s = 0;
107 : else
108 : { /* m0 | N0 */
109 3036 : GEN fa0 = famatsmall_divexact(n0, factoru(m0));
110 3036 : Q = Flj_mulu_pre(Q, m0, a4, p, pi);
111 3036 : n_s = Flj_order_ufact(Q, N0 / m0, fa0, a4, p, pi);
112 : }
113 300669 : if (n_s == 0)
114 : { /* If m0 divides N1 and m1 divides N0 and N0 < N1, then swap */
115 131418 : if (swapped || N1 % m0 || N0 % m1) return gc_long(ltop,0);
116 103571 : swapspec(n0, n1, N0, N1); swapped = 1; continue;
117 : }
118 :
119 169251 : m0 *= n_s; a1 = (2 * p + 2) % m1;
120 332994 : for (x = ceildivuu(hasse_low, m0) * m0; x <= hasse_high; x += m0)
121 192291 : if ((x % m1) == a1 && x != N0 && x != N1) break;
122 : /* every x was either N0 or N1, so we return true */
123 169251 : 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 28547 : if (first) { Fl_elltwist_disc(a4, a6, T, p, &a4t, &a6t); first = 0; }
127 28547 : lswap(a4, a4t);
128 28547 : lswap(a6, a6t);
129 28547 : lswap(m0, m1); set_avma(av);
130 : }
131 : }
132 :
133 : static GEN
134 868171 : random_FleV(GEN x, GEN a6, ulong p, ulong pi)
135 3997483 : { 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 420819 : torsion_constraint(struct torctab_rec *torctab, long ltorc, double tormod[], long n, long m)
648 : {
649 420819 : long i, b = -1;
650 420819 : double rb = -1.;
651 50267544 : for (i = 0 ; i < ltorc ; i++)
652 : {
653 49846725 : struct torctab_rec *ti = torctab + i;
654 49846725 : if ( ! (n%ti->m) && ( !ti->fix2 || (n%(2*ti->m)) ) && ( ! ti->fix3 || (n%(3*ti->m)) ) )
655 4523166 : if ( n == m || ( ! (m%ti->m) && ( !ti->fix2 || (m%(2*ti->m)) ) && ( ! ti->fix3 || (m%(3*ti->m)) ) ) )
656 : {
657 3632559 : double ri = ti->rating*tormod[ti->N];
658 3632559 : if ( b < 0 || ri < rb ) { b = i; rb = ri; }
659 : }
660 : }
661 420819 : if (b < 0) pari_err_BUG("find_rating");
662 419948 : return b;
663 : }
664 :
665 : /* p > 3 prime */
666 : static void
667 140777 : 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 140777 : switch(p % 12)
674 : {
675 33111 : case 11:torctab = torctab1; ltorc = numberof(torctab1); break;
676 32491 : case 5: torctab = torctab2; ltorc = numberof(torctab2); break;
677 41042 : case 7: torctab = torctab3; ltorc = numberof(torctab3); break;
678 34133 : default/*1*/: torctab = torctab4; ltorc = numberof(torctab4); break;
679 : }
680 4645237 : for ( i = 0 ; i < 32 ; i++ ) tormod[i] = 1.0;
681 140777 : if (p%5 == 1) tormod[5] = tormod[10] = tormod[15] = 6.0/5.0;
682 140777 : if (p%7 == 1) tormod[7] = tormod[14] = 8.0/7.0;
683 140777 : if (p%11 == 1) tormod[11] = 12.0/11.0;
684 140777 : if (p%13 == 1) tormod[13] = 14.0/13.0;
685 140777 : if (p%17 == 1) tormod[17] = 18.0/17.0;
686 140777 : if (p%19 == 1) tormod[19] = 20.0/19.0;
687 140777 : if (p%23 == 1) tormod[23] = 24.0/23.0;
688 140777 : if (p%29 == 1) tormod[29] = 30.0/29.0;
689 140777 : if (p%31 == 1) tormod[31] = 32.0/31.0;
690 :
691 140777 : n1 = p+1-t;
692 140777 : n2 = p+1+t;
693 140777 : b1 = torsion_constraint(torctab, ltorc, tormod, n1, n1);
694 140477 : b2 = torsion_constraint(torctab, ltorc, tormod, n2, n2);
695 140569 : b12 = torsion_constraint(torctab, ltorc, tormod, n1, n2);
696 140983 : if (b1 > b2) {
697 65654 : if (torctab[b2].rating / TWIST_DOUBLE_RATIO > torctab[b12].rating)
698 17619 : *ptwist = 3;
699 : else
700 48035 : *ptwist = 2;
701 : } else
702 75329 : if (torctab[b1].rating / TWIST_DOUBLE_RATIO > torctab[b12].rating)
703 24008 : *ptwist = 3;
704 : else
705 51321 : *ptwist = 1;
706 140983 : b = *ptwist ==1 ? b1: *ptwist ==2 ? b2: b12;
707 140983 : *ptor = torctab[b].N;
708 140983 : }
709 :
710 : /* This is Sutherland 2009 Algorithm 1.1 */
711 : static long
712 140938 : find_j_inv_with_given_trace(
713 : ulong *j_t, norm_eqn_t ne, long rho_inv, long max_curves)
714 : {
715 140938 : pari_sp ltop = avma, av;
716 140938 : long tested = 0, t = ne->t, batch_size, N0, N1, hasse_low, hasse_high, i;
717 : GEN n0, n1, A4, A6, tx, ty;
718 140938 : ulong p = ne->p, pi = ne->pi, p1 = p + 1, a4, a6, m, N;
719 : int twist;
720 :
721 140938 : 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 140887 : N0 = (long)p1 - t; n0 = factoru(N0);
727 140774 : N1 = (long)p1 + t; n1 = factoru(N1);
728 140765 : best_torsion_constraint(p, t, &twist, &m);
729 140952 : switch(twist)
730 : {
731 51311 : case 1: N = N0; break;
732 48035 : case 2: N = N1; break;
733 41606 : 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 140952 : batch_size = 1.0 + rho_inv / (2.0 * m);
739 140952 : A4 = cgetg(batch_size + 1, t_VECSMALL);
740 140723 : A6 = cgetg(batch_size + 1, t_VECSMALL);
741 140851 : tx = cgetg(batch_size + 1, t_VECSMALL);
742 140862 : ty = cgetg(batch_size + 1, t_VECSMALL);
743 :
744 140875 : dbg_printf(2)(" Selected torsion constraint m = %lu and batch "
745 : "size = %ld\n", m, batch_size);
746 140875 : hasse_bounds(&hasse_low, &hasse_high, p);
747 140886 : av = avma;
748 867618 : while (max_curves <= 0 || tested < max_curves)
749 : {
750 : GEN Pp1, Pt;
751 867618 : random_curves_with_m_torsion((ulong *)(A4 + 1), (ulong *)(A6 + 1),
752 867618 : (ulong *)(tx + 1), (ulong *)(ty + 1),
753 : batch_size, m, p, pi);
754 868183 : Pp1 = random_FleV(A4, A6, p, pi);
755 868408 : Pt = gcopy(Pp1);
756 868568 : FleV_mulu_pre_inplace(Pp1, N, A4, p, pi);
757 867807 : if (twist >= 3) FleV_mulu_pre_inplace(Pt, t, A4, p, pi);
758 3627014 : for (i = 1; i <= batch_size; ++i) {
759 2900223 : ++tested;
760 2900223 : a4 = A4[i];
761 2900223 : a6 = A6[i]; if (a4 == 0 || a6 == 0) continue;
762 :
763 2898016 : if (( (twist >= 3 && mael(Pp1,i,1) == mael(Pt,i,1))
764 2851950 : || (twist < 3 && umael(Pp1,i,1) == p))
765 168883 : && test_curve_order(ne, a4, a6, N0,N1, n0,n1, hasse_low,hasse_high)) {
766 141132 : *j_t = Fl_ellj_pre(a4, a6, p, pi);
767 141141 : return gc_long(ltop, tested);
768 : }
769 : }
770 726791 : 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 5571 : next_generator(GEN DD, long D, ulong u, long filter, GEN *genred, long *P)
779 : {
780 5571 : pari_sp av = avma;
781 5571 : ulong p = (ulong)*P;
782 : while (1)
783 : {
784 8807 : p = unextprime(p + 1);
785 8807 : if (p > LONG_MAX) pari_err_BUG("next_generator");
786 8807 : if (kross(D, (long)p) != -1 && u % p != 0 && filter % p != 0)
787 : {
788 5571 : GEN gen = primeform_u(DD, p);
789 : /* If gen is in the principal class, skip it */
790 5571 : *genred = qfi_red(gen);
791 5571 : if (!equali1(gel(*genred,1))) { *P = (long)p; return gen; }
792 0 : set_avma(av);
793 : }
794 : }
795 : }
796 :
797 : INLINE long *
798 5641 : evec_ri_mutate(long r[], long i)
799 5641 : { return r + (i * (i - 1) >> 1); }
800 :
801 : INLINE const long *
802 13191 : evec_ri(const long r[], long i)
803 13191 : { 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 21047 : evec_reduce(long e[], const long n[], const long r[], long k)
809 : {
810 : long i, j, q;
811 : const long *ri;
812 21047 : if (!k) return;
813 50168 : for (i = k - 1; i > 0; i--) {
814 29121 : if (e[i] >= n[i]) {
815 8510 : q = e[i] / n[i];
816 8510 : ri = evec_ri(r, i);
817 21495 : for (j = 0; j < i; j++) e[j] += q * ri[j];
818 8510 : e[i] -= q * n[i];
819 : }
820 : }
821 21047 : 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 10543 : evec_to_index(const long e[], const long m[], long k)
840 : {
841 10543 : long i, index = e[0];
842 25237 : for (i = 1; i < k; i++) index += e[i] * m[i - 1];
843 10543 : return index;
844 : }
845 :
846 : INLINE void
847 12220 : evec_copy(long f[], const long e[], long k)
848 : {
849 : long i;
850 28759 : for (i = 0; i < k; ++i) f[i] = e[i];
851 12220 : }
852 :
853 : INLINE void
854 4979 : evec_clear(long e[], long k)
855 : {
856 : long i;
857 12748 : for (i = 0; i < k; ++i) e[i] = 0;
858 4979 : }
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 5550 : evec_order(const long e[], const long n[], const long r[], long k)
903 : {
904 5550 : pari_sp av = avma;
905 5550 : long *f = new_chunk(k);
906 : long i, j, o, m;
907 :
908 5550 : evec_copy(f, e, k);
909 14272 : for (o = 1, i = k - 1; i >= 0; i--) if (f[i])
910 : {
911 6499 : m = n[i] / ugcd(f[i], n[i]);
912 17805 : for (j = 0; j < k; j++) f[j] *= m;
913 6499 : evec_reduce(f, n, r, k);
914 6499 : o *= m;
915 : }
916 5550 : 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 4349 : evec_orders(long o[], const long n[], const long r[], long k)
923 : {
924 4349 : pari_sp av = avma;
925 4349 : long i, *e = new_chunk(k);
926 :
927 4349 : evec_clear(e, k);
928 9899 : for (i = 0; i < k; i++) {
929 5550 : e[i] = 1;
930 5550 : if (i) e[i - 1] = 0;
931 5550 : o[i] = evec_order(e, n, r, k);
932 : }
933 4349 : set_avma(av);
934 4349 : }
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 1372 : index_to_evec(long e[], long index, const long m[], long k)
947 : {
948 : long i;
949 4130 : for (i = k - 1; i > 0; --i) {
950 2758 : e[i] = index / m[i - 1];
951 2758 : index -= e[i] * m[i - 1];
952 : }
953 1372 : e[0] = index;
954 1372 : }
955 :
956 : INLINE void
957 4349 : evec_n_to_m(long m[], const long n[], long k)
958 : {
959 : long i;
960 4349 : m[0] = n[0];
961 5550 : for (i = 1; i < k; ++i) m[i] = m[i - 1] * n[i];
962 4349 : }
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 4889 : upper_bound_on_classpoly_coeffs(long D, long h, GEN qfinorms)
977 : {
978 4889 : double B, logbinom, lnMk, lnMh = 0, C = 2114.567, t = M_PI * sqrt((double)-D);
979 4889 : ulong maxak = 0;
980 : long k, m;
981 :
982 80360 : for (k = 1, B = 0.0; k <= h; ++k)
983 : {
984 75471 : ulong ak = uel(qfinorms, k);
985 75471 : double tk = t / ak;
986 75471 : lnMk = tk + log(1.0 + C * exp(-tk));
987 75471 : B += lnMk;
988 75471 : if (ak > maxak) { maxak = ak; lnMh = lnMk; }
989 : }
990 4889 : m = floor((h + 1)/(exp(lnMh) + 1.0));
991 : /* log(binom(h, m)); 0 unless D <= -1579751 */
992 4889 : logbinom = (m > 0 && m < h)? logfac(h) - logfac(m) - logfac(h - m): 0;
993 4889 : 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 = qfi_red(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 280 : orient_pcp(GEN G, long *ni, long D, long u, hashtable *tbl)
1056 : {
1057 280 : pari_sp av = avma;
1058 : /* 199 seems to suffice, but can be increased if necessary */
1059 : enum { MAX_ORIENT_P = 199 };
1060 280 : 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 280 : long i, *ps = pcp_get_orient_p(G), *qs = pcp_get_orient_q(G), *reps = pcp_get_orient_reps(G);
1062 280 : long *ef, *e, *ei, *f, k = pcp_get_k(G), lvl = modinv_level(pcp_get_inv(G));
1063 280 : GEN DD = stoi(D);
1064 280 : long L0 = pcp_get_L0(G);
1065 :
1066 280 : memset(ps, 0, k * sizeof(long));
1067 280 : memset(qs, 0, k * sizeof(long));
1068 280 : memset(reps, 0, k * k * sizeof(long));
1069 :
1070 364 : for (i = 0; i < k; ++i) { ps[i] = -1; if (o[i] > 2) break; }
1071 378 : for (++i; i < k; ++i) ps[i] = (o[i] > 2) ? 0 : -1; /* ps[i] = -!(o[i] > 2); */
1072 :
1073 280 : e = new_chunk(k);
1074 280 : ei = new_chunk(k);
1075 280 : f = new_chunk(k);
1076 :
1077 728 : for (i = 0; i < k; ++i) {
1078 : long p;
1079 525 : 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 273 : if (ni) {
1112 273 : GEN N = qfb_nform(D, *ni);
1113 273 : hashentry *he = hash_search(tbl, N);
1114 273 : if (!he) pari_err_BUG("orient_pcp");
1115 273 : *ni = (long)he->val;
1116 : }
1117 273 : 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 4370 : classgp_pcp_check_generators(const long *n, long *r, long k, long L0)
1128 : {
1129 4370 : pari_sp av = avma;
1130 : long *e1, i, i0, j, s;
1131 : const long *ei;
1132 :
1133 4370 : s = !!L0;
1134 4370 : e1 = new_chunk(k);
1135 :
1136 5431 : for (i = s + 1; i < k; i++) {
1137 1082 : if (n[i] != 2) continue;
1138 647 : ei = evec_ri(r, i);
1139 906 : for (j = s; j < i; j++)
1140 703 : if (ei[j]) break;
1141 647 : if (j == i) continue;
1142 1042 : for (i0 = s; i0 < i; i0++) {
1143 619 : 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 4349 : return gc_long(av,-1);
1163 : }
1164 :
1165 : /* This is Sutherland 2009, Algorithm 2.2 (p16). */
1166 : static GEN
1167 4889 : 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 4889 : 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 4889 : L = zero_zv(MAX_GENS);
1177 4889 : m = zero_zv(MAX_GENS);
1178 4889 : n = zero_zv(MAX_GENS);
1179 4889 : o = zero_zv(MAX_GENS);
1180 4889 : r = zero_zv(MAX_GENS * (MAX_GENS-1) / 2);
1181 4889 : orient_p = zero_zv(MAX_GENS);
1182 4889 : orient_q = zero_zv(MAX_GENS);
1183 4889 : orient_reps = zero_zv(MAX_GENS*MAX_GENS);
1184 4889 : 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 4889 : av = avma;
1188 4889 : if (h == 1)
1189 : {
1190 547 : *height = upper_bound_on_classpoly_coeffs(D, h, mkvecsmall(1));
1191 547 : return gc_const(av, G); /* no need to set *ni when h = 1 */
1192 : }
1193 4342 : if (!modinv_is_double_eta(inv) || !modinv_ramified(D, inv, &L0)) L0 = 0;
1194 4342 : enum_cnt = h / (1 + !!L0);
1195 4342 : GLfilter = ulcm(Lfilter, lvl);
1196 4342 : DD = stoi(D); av2 = avma;
1197 : while (1) {
1198 4370 : 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 4370 : tbl = hash_create_GEN(h, 1);
1202 4370 : ident = primeform(DD, gen_1);
1203 4370 : hash_insert(tbl, ident, (void*)0);
1204 :
1205 4370 : T = vectrunc_init(h + 1);
1206 4370 : vectrunc_append(T, ident);
1207 4370 : nelts = 1;
1208 4370 : curr_p = 1;
1209 :
1210 10137 : while (nelts < h) {
1211 5767 : GEN gamma_i = NULL, beta;
1212 : hashentry *e;
1213 5767 : long N = lg(T)-1, ri = 1;
1214 :
1215 5767 : if (k == MAX_GENS) pari_err_IMPL("classgp_pcp");
1216 :
1217 5767 : if (nelts == 1 && L0) {
1218 196 : curr_p = L0;
1219 196 : gamma_i = qfb_nform(D, curr_p);
1220 196 : beta = qfi_red(gamma_i);
1221 196 : if (equali1(gel(beta, 1))) { curr_p = 1; gamma_i = NULL; }
1222 : }
1223 5767 : if (!gamma_i)
1224 5571 : gamma_i = next_generator(DD, D, u, GLfilter, &beta, &curr_p);
1225 62683 : while ((e = hash_search(tbl, beta)) == NULL) {
1226 : long j;
1227 131075 : for (j = 1; j <= N; ++j) {
1228 74159 : GEN t = qfbcomp_i(beta, gel(T, j));
1229 74159 : vectrunc_append(T, t);
1230 74159 : hash_insert(tbl, t, (void*)tbl->nb);
1231 : }
1232 56916 : beta = qfbcomp_i(beta, gamma_i);
1233 56916 : ++ri;
1234 : }
1235 5767 : if (ri > 1) {
1236 : long j, si;
1237 5592 : L[k+1] = curr_p;
1238 5592 : n[k+1] = ri;
1239 5592 : nelts *= ri;
1240 :
1241 : /* This is to reset the curr_p counter when we have L0 != 0
1242 : * in the first position of L. */
1243 5592 : if (curr_p == L0) curr_p = 1;
1244 :
1245 5592 : N = 1;
1246 5592 : si = (long)e->val;
1247 7199 : for (j = 1; j <= k; ++j) {
1248 1607 : evec_ri_mutate(r, k)[j] = (si / N) % n[j];
1249 1607 : N *= n[j];
1250 : }
1251 5592 : ++k;
1252 : }
1253 : }
1254 :
1255 4370 : if ((i = classgp_pcp_check_generators(n+1, r+1, k, L0)) < 0) {
1256 4349 : mael(G,2,1) = L0;
1257 4349 : mael(G,2,2) = k;
1258 4349 : mael(G,2,3) = enum_cnt;
1259 4349 : evec_orders(o+1, n+1, r+1, k);
1260 4349 : evec_n_to_m(m+1, n+1, k);
1261 4349 : if (!orient || orient_pcp(G, ni, D, u, tbl)) break;
1262 7 : GLfilter *= L[1];
1263 : } else {
1264 21 : GLfilter = umuluu_or_0(GLfilter, L[i+1]);
1265 21 : if (!GLfilter) pari_err_IMPL("classgp_pcp");
1266 : }
1267 28 : set_avma(av2);
1268 : }
1269 4342 : v = cgetg(h + 1, t_VECSMALL);
1270 4342 : v[1] = 1;
1271 76681 : for (i = 2; i <= h; ++i) uel(v,i) = itou(gmael(T,i,1));
1272 4342 : *height = upper_bound_on_classpoly_coeffs(D, enum_cnt, v);
1273 :
1274 : /* The norms of the last one or two generators. */
1275 4342 : L1 = L[k];
1276 4342 : L2 = k > 1 ? L[k - 1] : 1;
1277 : /* 4 * L1^2 * L2^2 must fit in a ulong */
1278 4342 : if (2 * (1 + log2(L1) + log2(L2)) >= BITS_IN_LONG)
1279 0 : pari_err_IMPL("classgp_pcp");
1280 4342 : if (L0 && (L[1] != L0 || o[1] != 2)) pari_err_BUG("classgp_pcp");
1281 4342 : return gc_const(av, G);
1282 : }
1283 :
1284 : /* SECTION: Functions for calculating class polynomials. */
1285 : static const long SMALL_PRIMES[11] = { 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31 };
1286 :
1287 : /* Encodes prime divisors of smooth integers <= 1200
1288 : P = primes(11); V = vector(31, i, vecsearch(P, i));
1289 : vfactor(v) =
1290 : { if (v == 1, return (0));
1291 : my(q = factor(v)[,1]);
1292 : if (vecmax(q) > 31, return (-1)); \\ not smooth
1293 : sum(i = 1, #P, 1 << (V[q[i]]-1)); }
1294 : vector(1200, v, vfactor(v)) */
1295 : static const long
1296 : SMOOTH_INTS[] = { -1, /* 0 */
1297 : 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,
1298 : 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,
1299 : -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,
1300 : 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,
1301 : 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,
1302 : -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,
1303 : 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,
1304 : 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,
1305 : 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,
1306 : 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,
1307 : -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,
1308 : -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,
1309 : -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,
1310 : 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,
1311 : 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,
1312 : -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,
1313 : 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,
1314 : -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,
1315 : 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,
1316 : 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,
1317 : -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,
1318 : 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,
1319 : -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,
1320 : 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,
1321 : 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,
1322 : -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,
1323 : -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,
1324 : 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,
1325 : 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,
1326 : -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,
1327 : -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,
1328 : 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,
1329 : -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,
1330 : -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,
1331 : -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,
1332 : 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,
1333 : 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,
1334 : -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,
1335 : -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,
1336 : -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,
1337 : -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,
1338 : 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,
1339 : -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,
1340 : -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,
1341 : -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,
1342 : -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,
1343 : -1,1153,-1,-1,-1,-1,40,-1,-1,-1,-1,19,-1,77,-1,-1,-1,-1,-1,289,138,-1,-1,7
1344 : };
1345 :
1346 : /* Upper bound for H(v^2 d) / H(d) = \prod_{p | v} (p + 1) / (p - 1)
1347 : * We actually store ceil(128 * bound)
1348 : P = primes(11); V = vector(31, i, vecsearch(P, i));
1349 : hbound(v) =
1350 : { my(q = factor(v)[,1]); if (q && vecmax(q) > 31, return (0));
1351 : ceil(prod(i = 1, #q, (q[i] + 1)/(q[i] - 1), 1.) * 128); }
1352 : vector(1200, v, hbound(v))
1353 : */
1354 : static const long HURWITZ_RATIO[] = {
1355 : 128,384,256,384,192,768,171,384,256,576,154,768,150,512,384,
1356 : 384,144,768,143,576,342,461,140,768,192,448,256,512,138,1152,
1357 : 137,384,308,432,256,768,0,427,299,576,0,1024,0,461,384,419,0,
1358 : 768,171,576,288,448,0,768,231,512,285,412,0,1152,0,410,342,
1359 : 384,224,922,0,432,280,768,0,768,0,0,384,427,205,896,0,576,
1360 : 256,0,0,1024,216,0,275,461,0,1152,200,419,274,0,214,768,0,
1361 : 512,308,576,0,864,0,448,512,0,0,768,0,692,0,512,0,854,210,
1362 : 412,299,0,192,1152,154,0,0,410,192,1024,0,384,0,672,0,922,
1363 : 190,0,384,432,0,838,0,768,0,0,180,768,206,0,342,0,0,1152,0,
1364 : 427,288,615,205,896,0,0,0,576,187,768,0,0,461,0,0,1024,150,
1365 : 648,285,0,0,823,256,461,0,0,0,1152,0,598,0,419,0,820,173,0,
1366 : 342,640,0,768,0,0,448,512,0,922,0,576,0,0,183,864,0,0,280,
1367 : 448,171,1536,0,0,0,0,0,768,183,0,0,692,168,0,0,512,384,0,
1368 : 0,854,0,629,410,412,0,896,0,0,0,576,0,1152,0,461,256,0,256,
1369 : 0,166,410,0,576,0,1024,168,0,432,384,0,0,0,672,275,0,0,922,
1370 : 0,569,0,0,0,1152,0,432,399,0,231,838,0,0,274,768,0,0,0,0,
1371 : 427,538,0,768,144,618,0,0,0,1024,0,0,308,0,163,1152,0,0,0,
1372 : 427,0,864,0,615,0,615,0,896,0,0,512,0,0,0,165,576,0,559,
1373 : 160,768,224,0,0,0,0,1383,0,0,0,0,0,1024,0,448,0,648,164,
1374 : 854,171,0,419,0,0,823,0,768,299,461,0,0,0,0,384,0,0,1152,
1375 : 143,0,308,598,0,0,0,419,0,0,0,820,0,519,384,0,160,1024,0,
1376 : 640,0,0,0,768,308,0,0,0,0,1344,158,512,0,0,0,922,0,0,380,
1377 : 576,0,0,160,0,384,549,0,864,0,0,0,0,0,838,0,448,0,512,0,
1378 : 1536,0,0,0,0,216,0,0,0,359,0,0,768,0,547,412,0,156,0,0,
1379 : 692,342,504,0,0,0,0,0,512,0,1152,0,0,0,0,299,854,0,0,288,
1380 : 629,0,1229,0,412,410,0,0,896,0,0,0,0,0,0,214,576,0,0,0,
1381 : 1152,0,0,373,461,0,768,0,0,0,768,0,0,155,498,461,410,0,0,
1382 : 0,576,0,0,0,1024,0,503,299,0,0,1296,0,384,285,0,0,0,0,0,
1383 : 0,672,0,823,0,0,512,0,154,922,140,0,0,569,0,0,0,0,0,0,
1384 : 205,1152,0,0,0,432,0,1195,0,0,0,692,153,838,0,0,0,0,0,820,
1385 : 0,768,346,0,0,0,0,0,342,0,0,1280,0,538,0,0,210,768,0,432,
1386 : 0,618,0,0,0,0,448,0,0,1024,152,0,0,0,0,922,288,0,0,489,0,
1387 : 1152,0,0,0,0,231,0,0,427,366,0,0,864,0,0,0,615,0,0,0,615,
1388 : 280,0,0,896,192,0,342,0,0,1536,0,0,0,0,0,0,200,494,0,576,
1389 : 0,0,0,559,0,480,0,768,0,672,365,0,0,0,0,0,0,0,0,1383,0,
1390 : 0,336,0,285,0,150,0,0,0,0,1024,0,0,384,448,0,0,0,648,0,
1391 : 492,0,854,0,512,0,0,0,1257,0,0,410,0,0,823,0,0,0,768,0,
1392 : 896,0,461,0,0,0,0,0,0,0,0,149,1152,269,0,0,0,0,1152,0,
1393 : 427,0,0,206,922,0,598,256,0,0,0,0,0,512,419,0,0,0,0,332,
1394 : 0,0,820,0,0,0,519,0,1152,0,0,0,480,0,1024,0,0,336,640,0,
1395 : 0,0,0,432,0,0,768,0,922,0,0,0,0,205,0,0,0,0,1344,0,472,
1396 : 275,512,0,0,0,0,0,0,0,922,0,0,0,0,0,1138,0,576,0,0,0,0,
1397 : 280,478,0,0,0,1152,0,549,0,0,0,864,0,0,399,0,0,0,0,0,461,
1398 : 0,0,838,0,0,0,448,192,0,0,512,274,0,0,1536,138,0,0,0,224,
1399 : 0,205,0,0,648,0,0,0,0,427,0,0,1076,0,0,0,0,0,768,0,0,
1400 : 288,547,0,1235,0,0,0,466,256,0,0,0,0,692,0,1024,0,504,0,0,
1401 : 0,0,0,0,308,0,0,0,0,512,326,0,147,1152,0,0,0,0,0,0,0,0,
1402 : 0,896,0,854,0,0,0,0,0,864,0,629,0,0,0,1229,0,0,0,412,0,
1403 : 1229,190,0,0,0,260,896,0,0,0,0,0,0,0,0,512,0,0,0,0,640,
1404 : 0,576,0,0,0,0,330,0,0,1152,137,0,0,0,0,1118,0,461,320,0,
1405 : 0,768,0,0,448,0,0,0,0,768,0,0,0,0,0,463,0,498,0,1383,0,
1406 : 410,0,0,0,0,0,0,0,576,239,0,0,0,0,0,0,1024,0,0,0,503,0,
1407 : 896,275,0,0,0,0,1296,0,0,328,384,0,854,0,0,342,0,0,0,0,0,
1408 : 419,0,0,0,0,672,0,0,0,823,256,0,0,0,0,1536,0,0,299,461,0,
1409 : 922,0,419,0,0,0,0,0,569,0,0,0,0,0,0,384,0,0,0,0,0,0,
1410 : 615,0,1152,0,0,285,0,274,0,0,432,308,0,0,1195,0,0,0,0,0,
1411 : 0,0,692,0,458,0,838,252,0,0,0,0,0,0,0,0,0,0,820,0,0,0,
1412 : 768,0,1037,0,0,384,0,187,0,0,0,320,0,0,1024,0,0,0,0,0,
1413 : 1280,0,0,0,538,0,0,0,0,0,629,0,768,0,0,615,432,0,0,0,618,
1414 : 0,0,0,0,0,0,0,0,0,1344,0,0,315,0,0,1024,0,456,0,0,0,0,
1415 : 200,0,0,0,0,922,0,864,0,0,0,0,0,489,380,0,0,1152
1416 : };
1417 :
1418 : /* Hurwitz class number of Df^2, D < 0; h = classno(D), Faf = factoru(f) */
1419 : static double
1420 194617 : hclassno_wrapper(long h, long D, GEN Faf)
1421 : {
1422 : pari_sp av;
1423 194617 : if (lg(gel(Faf,1)) == 1) switch(D)
1424 : {
1425 0 : case -3: return 1. / 3;
1426 0 : case -4: return 1. / 2;
1427 28017 : default: return (double)h;
1428 : }
1429 166600 : av = avma;
1430 166600 : return (double)gc_long(av, h * uhclassnoF_fact(Faf, D));
1431 : }
1432 :
1433 : /* return factor(u*v) */
1434 : static GEN
1435 189728 : factor_uv(GEN fau, ulong v, ulong vfactors)
1436 : {
1437 : GEN P, E;
1438 : long i;
1439 189728 : if (!vfactors) return fau;
1440 165283 : P = gel(fau,1);
1441 165283 : E = gel(fau,2);
1442 368138 : for (i = 0; vfactors; i++, vfactors >>= 1)
1443 368138 : if (vfactors & 1UL)
1444 : {
1445 193086 : long p = SMALL_PRIMES[i];
1446 193086 : P = vecsmall_append(P, p);
1447 193086 : E = vecsmall_append(E, u_lvalrem(v, p, &v));
1448 193086 : if (v == 1) break;
1449 : }
1450 165283 : return famatsmall_reduce(mkmat2(P, E));
1451 : }
1452 :
1453 : /* This is Sutherland 2009, Algorithm 2.1 (p8); delta > 0 */
1454 : static GEN
1455 4889 : select_classpoly_prime_pool(double min_bits, double delta, GEN G)
1456 : { /* Sutherland defines V_MAX to be 1200 without saying why */
1457 4889 : const long V_MAX = 1200;
1458 4889 : double bits = 0.0, hurwitz, z;
1459 4889 : ulong t_size_lim, d = (ulong)-pcp_get_D(G);
1460 4889 : long ires, inv = pcp_get_inv(G);
1461 4889 : GEN fau = pcp_get_fau(G);
1462 : GEN res, t_min; /* t_min[v] = lower bound for the t we look at for that v */
1463 : #ifdef LONG_IS_64BIT
1464 4182 : long L = pcp_get_L(G)[!!pcp_get_L0(G)];
1465 : #endif
1466 :
1467 4889 : hurwitz = hclassno_wrapper(pcp_get_h(G), pcp_get_D0(G), fau);
1468 :
1469 4889 : res = cgetg(128+1, t_VEC);
1470 4889 : ires = 1;
1471 : /* Initialise t_min to be all 2's. This avoids trace 0 and trace 1 curves */
1472 4889 : t_min = const_vecsmall(V_MAX-1, 2);
1473 :
1474 : /* maximum possible trace = sqrt(2^BIL + D) */
1475 4889 : t_size_lim = 2.0 * sqrt((double)((1UL << (BITS_IN_LONG - 2)) - (d >> 2)));
1476 :
1477 4889 : for (z = d / (2.0 * hurwitz); ; z *= delta + 1.0)
1478 19556 : {
1479 24445 : double v_bound_aux = 4.0 * z * hurwitz; /* = 4 z H(d) */
1480 : ulong v;
1481 24445 : dbg_printf(1)("z = %.2f\n", z);
1482 214060 : for (v = 1; v < V_MAX; v++)
1483 : #ifdef LONG_IS_64BIT
1484 183738 : if (L<=2 || v%L)
1485 : #endif
1486 : {
1487 : ulong p, t, t_max, vfactors, v2d, vd;
1488 : double hurwitz_ratio_bound, max_p, H;
1489 : long ires0;
1490 : GEN faw;
1491 :
1492 209284 : if ((long)(vfactors = SMOOTH_INTS[v]) < 0) continue;
1493 209284 : hurwitz_ratio_bound = HURWITZ_RATIO[v] / 128.0;
1494 209284 : vd = v * d;
1495 209284 : if (vd >= v_bound_aux * hurwitz_ratio_bound) break;
1496 189728 : v2d = v * vd;
1497 189728 : faw = factor_uv(fau, v, vfactors);
1498 189728 : H = hclassno_wrapper(pcp_get_h(G), pcp_get_D0(G), faw);
1499 : /* t <= 2 sqrt(p) and p <= z H(v^2 d) and
1500 : * H(v^2 d) < v H(d) \prod_{p | v} (p+1)/(p-1)
1501 : * This last term is v * hurwitz * hurwitz_ratio_bound. */
1502 189728 : max_p = z * v * hurwitz * hurwitz_ratio_bound;
1503 189728 : t_max = 2.0 * sqrt(mindd((1UL<<(BITS_IN_LONG-2)) - (v2d>>2), max_p));
1504 189728 : t = t_min[v]; if ((t & 1) != (v2d & 1)) t++;
1505 189728 : p = (t * t + v2d) >> 2;
1506 189728 : ires0 = ires;
1507 10386567 : for (; t <= t_max; p += t+1, t += 2) /* 4p = t^2 + v^2*d */
1508 10196839 : if (modinv_good_prime(inv,p) && uisprime(p))
1509 : {
1510 289745 : if (ires == lg(res)) res = vec_lengthen(res, lg(res) << 1);
1511 289745 : gel(res, ires++) = mkvec2(mkvecsmall5(p,t,v,(long)(p/H),vfactors),
1512 : faw);
1513 289745 : bits += log2(p);
1514 : }
1515 189728 : t_min[v] = t;
1516 :
1517 189728 : if (ires - ires0) {
1518 70910 : dbg_printf(2)(" Found %lu primes for v = %lu.\n", ires - ires0, v);
1519 : }
1520 189728 : if (bits > min_bits) {
1521 4889 : dbg_printf(1)("Found %ld primes; total size %.2f bits.\n", ires-1,bits);
1522 4889 : setlg(res, ires); return res;
1523 : }
1524 : }
1525 19556 : if (uel(t_min,1) >= t_size_lim) {
1526 : /* exhausted all solutions that fit in ulong */
1527 0 : char *err = stack_sprintf("class polynomial of discriminant %ld", pcp_get_D(G));
1528 0 : pari_err(e_ARCH, err);
1529 : }
1530 : }
1531 : }
1532 :
1533 : static int
1534 1237943 : primecmp(void *data, GEN v1, GEN v2)
1535 1237943 : { (void)data; return cmpss(gel(v1,1)[4], gel(v2,1)[4]); }
1536 :
1537 : static long
1538 4889 : height_margin(long inv, long D)
1539 : {
1540 : (void)D;
1541 : /* NB: avs just uses a height margin of 256 for everyone and everything. */
1542 4889 : if (inv == INV_F) return 64; /* checked for discriminants up to -350000 */
1543 4735 : if (inv == INV_G2) return 5;
1544 4497 : if (inv != INV_J) return 256; /* TODO: This should be made more accurate */
1545 3895 : return 0;
1546 : }
1547 :
1548 : static GEN
1549 4889 : select_classpoly_primes(ulong *vfactors, ulong *biggest_v, double delta,
1550 : GEN G, double height)
1551 : {
1552 4889 : const long k = 2;
1553 4889 : pari_sp av = avma;
1554 4889 : long i, s, D = pcp_get_D(G), inv = pcp_get_inv(G);
1555 : ulong biggest_p;
1556 : double prime_bits, min_prime_bits, b;
1557 : GEN prime_pool;
1558 :
1559 4889 : s = modinv_height_factor(inv);
1560 4889 : b = height / s + height_margin(inv, D);
1561 4889 : dbg_printf(1)("adjusted height = %.2f\n", b);
1562 4889 : min_prime_bits = k * b;
1563 :
1564 4889 : prime_pool = select_classpoly_prime_pool(min_prime_bits, delta, G);
1565 :
1566 : /* FIXME: Apply torsion constraints */
1567 : /* FIXME: Rank elts of res according to cost/benefit ratio */
1568 4889 : gen_sort_inplace(prime_pool, NULL, primecmp, NULL);
1569 4889 : prime_bits = 0.0; biggest_p = *biggest_v = *vfactors = 0;
1570 139726 : for (i = 1; i < lg(prime_pool); i++)
1571 : {
1572 139726 : GEN q = gmael(prime_pool, i, 1);
1573 139726 : ulong p = q[1], v = q[3];
1574 139726 : *vfactors |= q[5];
1575 139726 : prime_bits += log2(p);
1576 139726 : if (p > biggest_p) biggest_p = p;
1577 139726 : if (v > *biggest_v) *biggest_v = v;
1578 139726 : if (prime_bits > b) break;
1579 : }
1580 4889 : dbg_printf(1)("Selected %ld primes; largest is %lu ~ 2^%.2f\n",
1581 : i, biggest_p, log2(biggest_p));
1582 4889 : return gc_GEN(av, vecslice(prime_pool, 1, i));
1583 : }
1584 :
1585 : /* This is Sutherland 2009 Algorithm 1.2. */
1586 : static long
1587 141239 : oneroot_of_classpoly(
1588 : ulong *j_endo, int *endo_cert, ulong j, norm_eqn_t ne, GEN jdb)
1589 : {
1590 141239 : pari_sp av = avma;
1591 : long nfactors, L_bound, i;
1592 141239 : ulong p = ne->p, pi = ne->pi;
1593 : GEN factors, vdepths;
1594 :
1595 141239 : if (j == 0 || j == 1728 % p) pari_err_BUG("oneroot_of_classpoly");
1596 :
1597 141239 : *endo_cert = 1;
1598 141239 : factors = gel(ne->faw, 1); nfactors = lg(factors) - 1;
1599 141239 : if (!nfactors) { *j_endo = j; return 1; }
1600 137429 : vdepths = gel(ne->faw, 2);
1601 :
1602 : /* FIXME: This should be bigger */
1603 137429 : L_bound = maxdd(log((double) -ne->D), (double)ne->v);
1604 :
1605 : /* Iterate over the primes L dividing w */
1606 343064 : for (i = 1; i <= nfactors; ++i) {
1607 214303 : pari_sp av2 = avma;
1608 : GEN phi;
1609 214303 : long jlvl, lvl_diff, depth = vdepths[i], L = factors[i];
1610 214303 : if (L > L_bound) { *endo_cert = 0; break; }
1611 :
1612 205635 : phi = polmodular_db_getp(jdb, L, p);
1613 :
1614 : /* TODO: See if I can reuse paths created in j_level_in_volcano()
1615 : * later in {ascend,descend}_volcano(), perhaps by combining the
1616 : * functions into one "adjust_level" function. */
1617 205631 : jlvl = j_level_in_volcano(phi, j, p, pi, L, depth);
1618 205621 : lvl_diff = z_lval(ne->u, L) - jlvl;
1619 205639 : if (lvl_diff < 0)
1620 : /* j's level is less than v(u) so we must ascend */
1621 127181 : j = ascend_volcano(phi, j, p, pi, jlvl, L, depth, -lvl_diff);
1622 78458 : else if (lvl_diff > 0)
1623 : /* otherwise j's level is greater than v(u) so we descend */
1624 1304 : j = descend_volcano(phi, j, p, pi, jlvl, L, depth, lvl_diff);
1625 205638 : set_avma(av2);
1626 : }
1627 : /* Prob(j has the wrong endomorphism ring) ~ \sum_{p|u_compl} 1/p
1628 : * (and u_compl > L_bound), so just return it and rely on detection code in
1629 : * enum_j_with_endo_ring(). Detection is that we hit a previously found
1630 : * j-invariant earlier than expected. OR, we evaluate class polynomials of
1631 : * the suborders at j and if any are zero then j must be chosen again. */
1632 137429 : set_avma(av); *j_endo = j; return j != 0 && j != 1728 % p;
1633 : }
1634 :
1635 : INLINE long
1636 7627 : vecsmall_isin_skip(GEN v, long x, long k)
1637 : {
1638 7627 : long i, l = lg(v);
1639 458535 : for (i = k; i < l; ++i)
1640 451384 : if (v[i] == x) return i;
1641 7151 : return 0;
1642 : }
1643 :
1644 : void
1645 178447 : norm_eqn_set(norm_eqn_t ne, long D, long t, long u, long v, GEN faw, ulong p)
1646 : {
1647 178447 : ne->D = D;
1648 178447 : ne->u = u;
1649 178447 : ne->t = t;
1650 178447 : ne->v = v;
1651 178447 : ne->faw = faw;
1652 178447 : ne->p = p;
1653 178447 : ne->pi = get_Fl_red(ne->p);
1654 178463 : ne->s2 = Fl_2gener_pre(ne->p, ne->pi);
1655 : /* select twisting parameter */
1656 354815 : do ne->T = random_Fl(p); while (krouu(ne->T, p) != -1);
1657 178310 : }
1658 :
1659 : INLINE ulong
1660 15230 : Flv_powsum_pre(GEN v, ulong n, ulong p, ulong pi)
1661 : {
1662 15230 : long i, l = lg(v);
1663 15230 : ulong psum = 0;
1664 374028 : for (i = 1; i < l; ++i)
1665 358798 : psum = Fl_add(psum, Fl_powu_pre(uel(v,i), n, p, pi), p);
1666 15230 : return psum;
1667 : }
1668 :
1669 : INLINE int
1670 4889 : modinv_has_sign_ambiguity(long inv)
1671 : {
1672 4889 : switch (inv) {
1673 574 : case INV_F:
1674 : case INV_F3:
1675 : case INV_W2W3E2:
1676 : case INV_W2W7E2:
1677 : case INV_W2W3:
1678 : case INV_W2W5:
1679 : case INV_W2W7:
1680 : case INV_W3W3:
1681 : case INV_W2W13:
1682 574 : case INV_W3W7: return 1;
1683 : }
1684 4315 : return 0;
1685 : }
1686 :
1687 : INLINE int
1688 3879 : modinv_units(int inv)
1689 3879 : { return modinv_is_double_eta(inv) || modinv_is_Weber(inv); }
1690 :
1691 : INLINE int
1692 10841 : adjust_signs(GEN js, ulong p, ulong pi, long inv, GEN T, long e)
1693 : {
1694 10841 : long negate = 0;
1695 10841 : long h = lg(js) - 1;
1696 14468 : if ((h & 1) && modinv_units(inv)) {
1697 3627 : ulong prod = Flv_prod_pre(js, p, pi);
1698 3627 : if (prod != p - 1) {
1699 2102 : if (prod != 1) pari_err_BUG("adjust_signs: constant term is not +/-1");
1700 2102 : negate = 1;
1701 : }
1702 : } else {
1703 : ulong tp, t;
1704 7214 : tp = umodiu(T, p);
1705 7214 : t = Flv_powsum_pre(js, e, p, pi);
1706 7214 : if (t == 0) return 0;
1707 7200 : if (t != tp) {
1708 3691 : if (Fl_neg(t, p) != tp) pari_err_BUG("adjust_signs: incorrect trace");
1709 3691 : negate = 1;
1710 : }
1711 : }
1712 10827 : if (negate) Flv_neg_inplace(js, p);
1713 10827 : return 1;
1714 : }
1715 :
1716 : static ulong
1717 140731 : find_jinv(
1718 : long *trace_tries, long *endo_tries, int *cert,
1719 : norm_eqn_t ne, long inv, long rho_inv, GEN jdb)
1720 : {
1721 140731 : long found, ok = 1;
1722 : ulong j, r;
1723 : do {
1724 : do {
1725 : long tries;
1726 140885 : ulong j_t = 0;
1727 : /* TODO: Set batch size according to expected number of tries and
1728 : * experimental cost/benefit analysis. */
1729 140885 : tries = find_j_inv_with_given_trace(&j_t, ne, rho_inv, 0);
1730 141239 : if (j_t == 0)
1731 0 : pari_err_BUG("polclass0: Couldn't find j-invariant with given trace.");
1732 141239 : dbg_printf(2)(" j-invariant %ld has trace +/-%ld (%ld tries, 1/rho = %ld)\n",
1733 : j_t, ne->t, tries, rho_inv);
1734 141239 : *trace_tries += tries;
1735 :
1736 141239 : found = oneroot_of_classpoly(&j, cert, j_t, ne, jdb);
1737 141240 : ++*endo_tries;
1738 141240 : } while (!found);
1739 :
1740 141189 : if (modinv_is_double_eta(inv))
1741 14788 : ok = modfn_unambiguous_root(&r, inv, j, ne, jdb);
1742 : else
1743 126401 : r = modfn_root(j, ne, inv);
1744 141189 : } while (!ok);
1745 141086 : return r;
1746 : }
1747 :
1748 : static GEN
1749 139490 : polclass_roots_modp(long *n_trace_curves,
1750 : norm_eqn_t ne, long rho_inv, GEN G, GEN db)
1751 : {
1752 : pari_sp av;
1753 139490 : ulong j = 0;
1754 139490 : long inv = pcp_get_inv(G), endo_tries = 0;
1755 : int endo_cert;
1756 139489 : GEN res, jdb, fdb, vshape = factoru(ne->v);
1757 :
1758 139328 : jdb = polmodular_db_for_inv(db, INV_J);
1759 139333 : fdb = polmodular_db_for_inv(db, inv);
1760 139343 : dbg_printf(2)("p = %ld, t = %ld, v = %ld\n", ne->p, ne->t, ne->v);
1761 139346 : av = avma;
1762 : do {
1763 140709 : j = find_jinv(n_trace_curves,&endo_tries,&endo_cert, ne, inv, rho_inv, jdb);
1764 141086 : res = enum_roots(j, ne, fdb, G, vshape);
1765 141085 : if (!res && endo_cert) pari_err_BUG("polclass_roots_modp");
1766 141085 : if (res && !endo_cert && vecsmall_isin_skip(res, res[1], 2))
1767 : {
1768 476 : set_avma(av);
1769 476 : res = NULL;
1770 : }
1771 141085 : } while (!res);
1772 :
1773 139722 : dbg_printf(2)(" j-invariant %ld has correct endomorphism ring "
1774 : "(%ld tries)\n", j, endo_tries);
1775 139722 : dbg_printf(4)(" all such j-invariants: %Ps\n", res);
1776 139722 : return res;
1777 : }
1778 :
1779 : INLINE int
1780 2083 : modinv_inverted_involution(long inv)
1781 2083 : { return modinv_is_double_eta(inv); }
1782 :
1783 : INLINE int
1784 2083 : modinv_negated_involution(long inv)
1785 : { /* determined by trial and error */
1786 2083 : return inv == INV_F || inv == INV_W3W5 || inv == INV_W3W7
1787 4166 : || inv == INV_W3W3 || inv == INV_W5W7;
1788 : }
1789 :
1790 : /* Return true iff Phi_L(j0, j1) = 0. */
1791 : INLINE long
1792 3656 : verify_edge(ulong j0, ulong j1, ulong p, ulong pi, long L, GEN fdb)
1793 : {
1794 3656 : pari_sp av = avma;
1795 3656 : GEN phi = polmodular_db_getp(fdb, L, p);
1796 3656 : GEN f = Flm_Fl_polmodular_evalx(phi, L, j1, p, pi);
1797 3656 : return gc_long(av, Flx_eval_pre(f, j0, p, pi) == 0);
1798 : }
1799 :
1800 : INLINE long
1801 672 : verify_2path(
1802 : ulong j1, ulong j2, ulong p, ulong pi, long L1, long L2, GEN fdb)
1803 : {
1804 672 : pari_sp av = avma;
1805 672 : GEN phi1 = polmodular_db_getp(fdb, L1, p);
1806 672 : GEN phi2 = polmodular_db_getp(fdb, L2, p);
1807 672 : GEN f = Flm_Fl_polmodular_evalx(phi1, L1, j1, p, pi);
1808 672 : GEN g = Flm_Fl_polmodular_evalx(phi2, L2, j2, p, pi);
1809 672 : GEN d = Flx_gcd(f, g, p);
1810 672 : long n = degpol(d);
1811 672 : if (n >= 2) n = Flx_nbroots(d, p);
1812 672 : return gc_long(av, n);
1813 : }
1814 :
1815 : static long
1816 6215 : oriented_n_action(
1817 : const long *ni, GEN G, GEN v, ulong p, ulong pi, GEN fdb)
1818 : {
1819 6215 : pari_sp av = avma;
1820 6215 : long i, j, k = pcp_get_k(G);
1821 6215 : long nr = k * (k - 1) / 2;
1822 6215 : const long *n = pcp_get_n(G), *m = pcp_get_m(G), *o = pcp_get_o(G), *r = pcp_get_r(G),
1823 6215 : *ps = pcp_get_orient_p(G), *qs = pcp_get_orient_q(G), *reps = pcp_get_orient_reps(G);
1824 6215 : long *signs = new_chunk(k);
1825 6215 : long *e = new_chunk(k);
1826 6215 : long *rels = new_chunk(nr);
1827 :
1828 6215 : evec_copy(rels, r, nr);
1829 :
1830 16746 : for (i = 0; i < k; ++i) {
1831 : /* If generator doesn't require orientation, continue; power rels already
1832 : * copied to *rels in initialisation */
1833 10531 : if (ps[i] <= 0) { signs[i] = 1; continue; }
1834 : /* Get rep of orientation element and express it in terms of the
1835 : * (partially) oriented presentation */
1836 6198 : for (j = 0; j < i; ++j) {
1837 4034 : long t = reps[i * k + j];
1838 4034 : e[j] = (signs[j] < 0 ? o[j] - t : t);
1839 : }
1840 2164 : e[j] = reps[i * k + j];
1841 3319 : for (++j; j < k; ++j) e[j] = 0;
1842 2164 : evec_reduce(e, n, rels, k);
1843 2164 : j = evec_to_index(e, m, k);
1844 :
1845 : /* FIXME: These calls to verify_edge recalculate powers of v[0]
1846 : * and v[j] over and over again, they also reduce Phi_{ps[i]} modulo p over
1847 : * and over again. Need to cache these things! */
1848 2164 : if (qs[i] > 1)
1849 336 : signs[i] =
1850 336 : (verify_2path(uel(v,1), uel(v,j+1), p, pi, ps[i], qs[i], fdb) ? 1 : -1);
1851 : else
1852 : /* Verify ps[i]-edge to orient ith generator */
1853 1828 : signs[i] =
1854 1828 : (verify_edge(uel(v,1), uel(v,j+1), p, pi, ps[i], fdb) ? 1 : -1);
1855 : /* Update power relation */
1856 6198 : for (j = 0; j < i; ++j) {
1857 4034 : long t = evec_ri(r, i)[j];
1858 4034 : e[j] = (signs[i] * signs[j] < 0 ? o[j] - t : t);
1859 : }
1860 5483 : while (j < k) e[j++] = 0;
1861 2164 : evec_reduce(e, n, rels, k);
1862 6198 : for (j = 0; j < i; ++j) evec_ri_mutate(rels, i)[j] = e[j];
1863 : /* TODO: This is a sanity check, can be removed if everything is working */
1864 8362 : for (j = 0; j <= i; ++j) {
1865 6198 : long t = reps[i * k + j];
1866 6198 : e[j] = (signs[j] < 0 ? o[j] - t : t);
1867 : }
1868 3319 : while (j < k) e[j++] = 0;
1869 2164 : evec_reduce(e, n, rels, k);
1870 2164 : j = evec_to_index(e, m, k);
1871 2164 : if (qs[i] > 1) {
1872 336 : if (!verify_2path(uel(v,1), uel(v, j+1), p, pi, ps[i], qs[i], fdb))
1873 0 : pari_err_BUG("oriented_n_action");
1874 : } else {
1875 1828 : if (!verify_edge(uel(v,1), uel(v, j+1), p, pi, ps[i], fdb))
1876 0 : pari_err_BUG("oriented_n_action");
1877 : }
1878 : }
1879 :
1880 : /* Orient representation of [N] relative to the torsor <signs, rels> */
1881 16746 : for (i = 0; i < k; ++i) e[i] = (signs[i] < 0 ? o[i] - ni[i] : ni[i]);
1882 6215 : evec_reduce(e, n, rels, k);
1883 6215 : return gc_long(av, evec_to_index(e,m,k));
1884 : }
1885 :
1886 : /* F = double_eta_raw(inv) */
1887 : INLINE void
1888 6215 : adjust_orientation(GEN F, long inv, GEN v, long e, ulong p, ulong pi)
1889 : {
1890 6215 : ulong j0 = uel(v, 1), je = uel(v, e);
1891 :
1892 6215 : if (!modinv_j_from_2double_eta(F, inv, j0, je, p, pi)) {
1893 2083 : if (modinv_inverted_involution(inv)) Flv_inv_pre_inplace(v, p, pi);
1894 2083 : if (modinv_negated_involution(inv)) Flv_neg_inplace(v, p);
1895 : }
1896 6215 : }
1897 :
1898 : static void
1899 574 : polclass_psum(
1900 : GEN *psum, long *d, GEN roots, GEN primes, GEN pilist, ulong h, long inv)
1901 : {
1902 : /* Number of consecutive CRT stabilisations before we assume we have
1903 : * the correct answer. */
1904 : enum { MIN_STAB_CNT = 3 };
1905 574 : pari_sp av = avma, btop;
1906 : GEN ps, psum_sqr, P;
1907 574 : long i, e, stabcnt, nprimes = lg(primes) - 1;
1908 :
1909 574 : if ((h & 1) && modinv_units(inv)) { *psum = gen_1; *d = 0; return; }
1910 322 : e = -1;
1911 322 : ps = cgetg(nprimes+1, t_VECSMALL);
1912 : do {
1913 364 : e += 2;
1914 8380 : for (i = 1; i <= nprimes; ++i)
1915 : {
1916 8016 : GEN roots_modp = gel(roots, i);
1917 8016 : ulong p = uel(primes, i), pi = uel(pilist, i);
1918 8016 : uel(ps, i) = Flv_powsum_pre(roots_modp, e, p, pi);
1919 : }
1920 364 : btop = avma;
1921 364 : psum_sqr = Z_init_CRT(0, 1);
1922 364 : P = gen_1;
1923 2079 : for (i = 1, stabcnt = 0; stabcnt < MIN_STAB_CNT && i <= nprimes; ++i)
1924 : {
1925 1715 : ulong p = uel(primes, i), pi = uel(pilist, i);
1926 1715 : ulong ps2 = Fl_sqr_pre(uel(ps, i), p, pi);
1927 1715 : ulong stab = Z_incremental_CRT(&psum_sqr, ps2, &P, p);
1928 : /* stabcnt = stab * (stabcnt + 1) */
1929 1715 : if (stab) ++stabcnt; else stabcnt = 0;
1930 1715 : if (gc_needed(av, 2)) (void)gc_all(btop, 2, &psum_sqr, &P);
1931 : }
1932 364 : if (stabcnt == 0 && nprimes >= MIN_STAB_CNT)
1933 0 : pari_err_BUG("polclass_psum");
1934 364 : } while (!signe(psum_sqr));
1935 :
1936 322 : if ( ! Z_issquareall(psum_sqr, psum)) pari_err_BUG("polclass_psum");
1937 :
1938 322 : dbg_printf(1)("Classpoly power sum (e = %ld) is %Ps; found with %.2f%% of the primes\n",
1939 0 : e, *psum, 100 * (i - 1) / (double) nprimes);
1940 322 : *psum = gc_upto(av, *psum);
1941 322 : *d = e;
1942 : }
1943 :
1944 : static GEN
1945 1505 : polclass_small_disc(long D, long inv, long vx)
1946 : {
1947 1505 : if (D == -3) return pol_x(vx);
1948 525 : if (D == -4) {
1949 525 : switch (inv) {
1950 364 : case INV_J: return deg1pol(gen_1, stoi(-1728), vx);
1951 161 : case INV_G2:return deg1pol(gen_1, stoi(-12), vx);
1952 0 : default: /* no other invariants for which we can calculate H_{-4}(X) */
1953 0 : pari_err_BUG("polclass_small_disc");
1954 : }
1955 : }
1956 0 : return NULL;
1957 : }
1958 :
1959 : static ulong
1960 4889 : quadnegclassnou(long D, long *pD0, GEN *pP, GEN *pE)
1961 : {
1962 4889 : ulong d = (ulong)-D, d0 = coredisc2u_fact(factoru(d), -1, pP, pE);
1963 4889 : ulong h = uquadclassnoF_fact(d0, -1, *pP, *pE) * quadclassnos(-(long)d0);
1964 4889 : *pD0 = -(long)d0; return h;
1965 : }
1966 :
1967 : GEN
1968 139631 : polclass_worker(GEN q, GEN G, GEN db)
1969 : {
1970 139631 : GEN T = cgetg(3, t_VEC), z, P = gel(q,1), faw = gel(q,2);
1971 139617 : long n_curves_tested = 0, t = P[2], v = P[3], rho_inv = P[4];
1972 139617 : ulong p = (ulong)P[1];
1973 139617 : pari_sp av = avma;
1974 139617 : norm_eqn_t ne; norm_eqn_set(ne, pcp_get_D(G), t, pcp_get_u(G), v, faw, p);
1975 139484 : z = polclass_roots_modp(&n_curves_tested, ne, rho_inv, G, db);
1976 139722 : gel(T,1) = gc_leaf(av, z);
1977 139725 : gel(T,2) = mkvecsmall3(ne->p, ne->pi, n_curves_tested); return T;
1978 : }
1979 :
1980 : GEN
1981 6394 : polclass0(long D, long inv, long vx, GEN *db)
1982 : {
1983 6394 : pari_sp av = avma;
1984 : GEN primes, H, plist, pilist, Pu, Eu;
1985 6394 : long n_curves_tested = 0, filter = 1, pending = 0, cnt = 0;
1986 : long D0, nprimes, s, i, j, del, ni, orient, h, p1, p2, k;
1987 : ulong u, vfactors, biggest_v;
1988 : GEN G, worker, vec;
1989 : double height;
1990 : static const double delta = 0.5;
1991 : struct pari_mt pt;
1992 :
1993 6394 : if (D >= -4) return polclass_small_disc(D, inv, vx);
1994 :
1995 4889 : h = quadnegclassnou(D, &D0, &Pu, &Eu);
1996 4889 : u = D == D0? 1: (ulong)sqrt(D / D0); /* u = \prod Pu[i]^Eu[i] */
1997 4889 : dbg_printf(1)("D = %ld, conductor = %ld, inv = %ld\n", D, u, inv);
1998 :
1999 4889 : ni = modinv_degree(&p1, &p2, inv);
2000 4889 : orient = modinv_is_double_eta(inv) && kross(D, p1) && p2>1 && kross(D, p2);
2001 :
2002 4889 : G = classgp_make_pcp(&height, &ni, h, D, D0, u, Pu, Eu, inv, filter, orient);
2003 4889 : primes = select_classpoly_primes(&vfactors, &biggest_v, delta, G, height);
2004 :
2005 : /* Prepopulate *db with all the modpolys we might need */
2006 4889 : if (u > 1)
2007 : { /* L_bound in oneroot_of_classpoly() */
2008 322 : long l = lg(Pu), maxL = maxdd(log((double) -D), (double)biggest_v);
2009 550 : for (i = 1; i < l; i++)
2010 : {
2011 364 : long L = Pu[i];
2012 364 : if (L > maxL) break;
2013 228 : polmodular_db_add_level(db, L, INV_J);
2014 : }
2015 : }
2016 19686 : for (i = 0; vfactors; ++i) {
2017 14797 : if (vfactors & 1UL)
2018 14051 : polmodular_db_add_level(db, SMALL_PRIMES[i], INV_J);
2019 14797 : vfactors >>= 1;
2020 : }
2021 4889 : if (p1 > 1) polmodular_db_add_level(db, p1, INV_J);
2022 4889 : if (p2 > 1) polmodular_db_add_level(db, p2, INV_J);
2023 4889 : s = !!pcp_get_L0(G); k = pcp_get_k(G);
2024 4889 : polmodular_db_add_levels(db, pcp_get_L(G) + s, k - s, inv);
2025 4889 : if (orient)
2026 : {
2027 273 : GEN orient_p = pcp_get_orient_p(G);
2028 273 : GEN orient_q = pcp_get_orient_q(G);
2029 714 : for (i = 0; i < k; ++i)
2030 : {
2031 441 : if (orient_p[i] > 1) polmodular_db_add_level(db, orient_p[i], inv);
2032 441 : if (orient_q[i] > 1) polmodular_db_add_level(db, orient_q[i], inv);
2033 : }
2034 : }
2035 4889 : nprimes = lg(primes) - 1;
2036 4889 : worker = snm_closure(is_entry("_polclass_worker"), mkvec2(G, *db));
2037 4889 : H = cgetg(nprimes + 1, t_VEC);
2038 4889 : plist = cgetg(nprimes + 1, t_VECSMALL);
2039 4889 : pilist = cgetg(nprimes + 1, t_VECSMALL);
2040 4889 : vec = cgetg(2, t_VEC);
2041 4889 : dbg_printf(0)("Calculating class polynomial of disc %ld and inv %ld:", D, inv);
2042 4889 : mt_queue_start_lim(&pt, worker, nprimes);
2043 153367 : for (i = 1; i <= nprimes || pending; i++)
2044 : {
2045 : long workid;
2046 : GEN done;
2047 148478 : if (i <= nprimes) gel(vec,1) = gel(primes,i);
2048 148478 : mt_queue_submit(&pt, i, i <= nprimes? vec: NULL);
2049 148478 : done = mt_queue_get(&pt, &workid, &pending);
2050 148478 : if (done)
2051 : {
2052 139726 : gel(H, workid) = gel(done,1);
2053 139726 : uel(plist, workid) = mael(done,2,1);
2054 139726 : uel(pilist, workid) = mael(done,2,2);
2055 139726 : n_curves_tested += mael(done,2,3);
2056 139726 : cnt++;
2057 139726 : if (DEBUGLEVEL>2 && (cnt & 3L)==0)
2058 0 : err_printf(" %ld%%", cnt*100/nprimes);
2059 : }
2060 : }
2061 4889 : mt_queue_end(&pt);
2062 4889 : dbg_printf(0)(" done\n");
2063 :
2064 4889 : if (orient) {
2065 273 : GEN nvec = new_chunk(k);
2066 273 : GEN fdb = polmodular_db_for_inv(*db, inv);
2067 273 : GEN F = double_eta_raw(inv);
2068 273 : index_to_evec((long *)nvec, ni, pcp_get_m(G), k);
2069 6488 : for (i = 1; i <= nprimes; ++i) {
2070 6215 : GEN v = gel(H, i);
2071 6215 : ulong p = uel(plist, i), pi = uel(pilist, i);
2072 6215 : long oni = oriented_n_action(nvec, G, v, p, pi, fdb);
2073 6215 : adjust_orientation(F, inv, v, oni + 1, p, pi);
2074 : }
2075 : }
2076 :
2077 4889 : if (modinv_has_sign_ambiguity(inv)) {
2078 : GEN psum;
2079 : long e;
2080 574 : polclass_psum(&psum, &e, H, plist, pilist, h, inv);
2081 11415 : for (i = 1; i <= nprimes; ++i) {
2082 10841 : GEN v = gel(H, i);
2083 10841 : ulong p = uel(plist, i), pi = uel(pilist, i);
2084 10841 : if (!adjust_signs(v, p, pi, inv, psum, e))
2085 14 : uel(plist, i) = 0;
2086 : }
2087 : }
2088 :
2089 144615 : for (i = 1, j = 1, del = 0; i <= nprimes; ++i)
2090 : {
2091 139726 : pari_sp av2 = avma;
2092 139726 : ulong p = uel(plist, i);
2093 : long l;
2094 : GEN v;
2095 139726 : if (!p) { del++; continue; }
2096 139712 : v = Flv_roots_to_pol(gel(H,i), p, vx); l = lg(v);
2097 139712 : *++v = evaltyp(t_VECSMALL) | _evallg(l-1); /* Flx_to_Flv inplace */
2098 139712 : uel(plist, j) = p;
2099 139712 : gel(H, j++) = gc_leaf(av2, v);
2100 : }
2101 4889 : setlg(H,nprimes+1-del);
2102 4889 : setlg(plist,nprimes+1-del);
2103 :
2104 4889 : dbg_printf(1)("Total number of curves tested: %ld\n", n_curves_tested);
2105 4889 : H = ncV_chinese_center(H, plist, NULL);
2106 4889 : dbg_printf(1)("Result height: %.2f\n", dbllog2(gsupnorm(H, DEFAULTPREC)));
2107 4889 : return gc_GEN(av, RgV_to_RgX(H, vx));
2108 : }
2109 :
2110 : void
2111 2800 : check_modinv(long inv)
2112 : {
2113 2800 : switch (inv) {
2114 2786 : case INV_J:
2115 : case INV_F:
2116 : case INV_F2:
2117 : case INV_F3:
2118 : case INV_F4:
2119 : case INV_G2:
2120 : case INV_W2W3:
2121 : case INV_F8:
2122 : case INV_W3W3:
2123 : case INV_W2W5:
2124 : case INV_W2W7:
2125 : case INV_W3W5:
2126 : case INV_W3W7:
2127 : case INV_W2W3E2:
2128 : case INV_W2W5E2:
2129 : case INV_W2W13:
2130 : case INV_W2W7E2:
2131 : case INV_W3W3E2:
2132 : case INV_W5W7:
2133 : case INV_W3W13:
2134 : case INV_ATKIN3:
2135 : case INV_ATKIN5:
2136 : case INV_ATKIN7:
2137 2786 : break;
2138 14 : default:
2139 14 : pari_err_DOMAIN("polmodular", "inv", "invalid invariant", stoi(inv), gen_0);
2140 : }
2141 2786 : }
2142 :
2143 : GEN
2144 2177 : polclass(GEN DD, long inv, long vx)
2145 : {
2146 : GEN db, H;
2147 : long D;
2148 :
2149 2177 : if (vx < 0) vx = 0;
2150 2177 : check_quaddisc_imag(DD, NULL, "polclass");
2151 2170 : check_modinv(inv);
2152 :
2153 2163 : D = itos(DD);
2154 2163 : if (!modinv_good_disc(inv, D))
2155 0 : pari_err_DOMAIN("polclass", "D", "incompatible with given invariant", stoi(inv), DD);
2156 :
2157 2163 : db = polmodular_db_init(inv);
2158 2163 : H = polclass0(D, inv, vx, &db);
2159 2163 : gunclone_deep(db); return H;
2160 : }
|
