Line data Source code
1 : /* Copyright (C) 2008 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 : /* This file is a C version by Bill Allombert of the 'ellsea' GP package
16 : * whose copyright statement is as follows:
17 : Authors:
18 : Christophe Doche <cdoche@math.u-bordeaux.fr>
19 : Sylvain Duquesne <duquesne@math.u-bordeaux.fr>
20 :
21 : Universite Bordeaux I, Laboratoire A2X
22 : For the AREHCC project, see http://www.arehcc.com/
23 :
24 : Contributors:
25 : Karim Belabas (code cleanup and package release, faster polynomial arithmetic)
26 :
27 : 'ellsea' is free software; you can redistribute it and/or modify it under the
28 : terms of the GNU General Public License as published by the Free Software
29 : Foundation. It is distributed in the hope that it will be useful, but WITHOUT
30 : ANY WARRANTY WHATSOEVER. */
31 :
32 : /* Extension to non prime finite fields by Bill Allombert 2012 */
33 :
34 : #include "pari.h"
35 : #include "paripriv.h"
36 :
37 : #define DEBUGLEVEL DEBUGLEVEL_ellsea
38 :
39 : static THREAD GEN modular_eqn;
40 :
41 : void
42 373050 : pari_set_seadata(GEN mod) { modular_eqn = mod; }
43 : GEN
44 371173 : pari_get_seadata(void) { return modular_eqn; }
45 :
46 : static char *
47 91 : seadata_filename(ulong ell)
48 91 : { return stack_sprintf("%s/seadata/sea%ld", pari_datadir, ell); }
49 :
50 : static GEN
51 91 : get_seadata(ulong ell)
52 : {
53 91 : pari_sp av = avma;
54 : GEN eqn;
55 91 : char *s = seadata_filename(ell);
56 91 : pariFILE *F = pari_fopengz(s);
57 91 : if (!F) return NULL;
58 35 : if (ell) /* large single polynomial */
59 7 : eqn = gp_read_stream(F->file);
60 : else
61 : { /* table of polynomials of small level */
62 28 : eqn = gp_readvec_stream(F->file);
63 28 : modular_eqn = eqn = gclone(eqn);
64 28 : set_avma(av);
65 : }
66 35 : pari_fclose(F);
67 35 : return eqn;
68 : }
69 :
70 : /*Builds the modular equation corresponding to the vector list. Shallow */
71 : static GEN
72 9968 : list_to_pol(GEN list, long vx, long vy)
73 : {
74 9968 : long i, l = lg(list);
75 9968 : GEN P = cgetg(l, t_VEC);
76 205121 : for (i = 1; i < l; i++)
77 : {
78 195153 : GEN L = gel(list,i);
79 195153 : if (typ(L) == t_VEC) L = RgV_to_RgX_reverse(L, vy);
80 195153 : gel(P, i) = L;
81 : }
82 9968 : return RgV_to_RgX_reverse(P, vx);
83 : }
84 :
85 : struct meqn {
86 : char type;
87 : GEN eq, eval;
88 : long vx,vy;
89 : };
90 :
91 : static GEN
92 10024 : seadata_cache(ulong ell)
93 : {
94 10024 : long n = uprimepi(ell)-1;
95 : GEN C;
96 10024 : if (!modular_eqn && !get_seadata(0))
97 56 : C = NULL;
98 9968 : else if (n && n < lg(modular_eqn))
99 9961 : C = gel(modular_eqn, n);
100 : else
101 7 : C = get_seadata(ell);
102 10024 : return C;
103 : }
104 : /* C = [prime level, type "A" or "C", pol. coeffs] */
105 : static void
106 9968 : seadata_parse(struct meqn *M, GEN C, long vx, long vy)
107 : {
108 9968 : M->type = *GSTR(gel(C,2));
109 9968 : M->eq = list_to_pol(gel(C,3), vx, vy);
110 9968 : }
111 : static void
112 10003 : get_modular_eqn(struct meqn *M, ulong ell, long vx, long vy)
113 : {
114 10003 : GEN C = seadata_cache(ell);
115 10003 : M->vx = vx;
116 10003 : M->vy = vy;
117 10003 : M->eval = gen_0;
118 10003 : if (C) seadata_parse(M, C, vx, vy);
119 : else
120 : {
121 56 : M->type = 'J'; /* j^(1/3) for ell != 3, j for 3 */
122 56 : M->eq = polmodular_ZXX(ell, ell==3? 0: 5, vx, vy);
123 : }
124 10003 : }
125 :
126 : GEN
127 35 : ellmodulareqn(long ell, long vx, long vy)
128 : {
129 35 : pari_sp av = avma;
130 : struct meqn meqn;
131 : GEN C;
132 35 : if (vx < 0) vx = 0;
133 35 : if (vy < 0) vy = 1;
134 35 : if (varncmp(vx,vy) >= 0)
135 7 : pari_err_PRIORITY("ellmodulareqn", pol_x(vx), ">=", vy);
136 28 : if (ell < 2 || !uisprime(ell))
137 7 : pari_err_PRIME("ellmodulareqn (level)", stoi(ell));
138 21 : C = seadata_cache(ell);
139 21 : if (!C) pari_err_FILE("seadata file", seadata_filename(ell));
140 21 : seadata_parse(&meqn, C, vx, vy);
141 21 : return gc_GEN(av, mkvec2(meqn.eq, meqn.type=='A'? gen_1: gen_0));
142 : }
143 :
144 : /***********************************************************************/
145 : /** **/
146 : /** FqE_group **/
147 : /** **/
148 : /***********************************************************************/
149 :
150 : static GEN
151 122 : Fq_to_Flx(GEN a4, GEN T, ulong p)
152 122 : { return typ(a4)==t_INT ? Z_to_Flx(a4, p, get_Flx_var(T)): ZX_to_Flx(a4, p); }
153 :
154 : /*FIXME: the name of the function does not quite match what it does*/
155 : static const struct bb_group *
156 980 : get_FqE_group(void ** pt_E, GEN a4, GEN a6, GEN T, GEN p)
157 : {
158 980 : if (!T) return get_FpE_group(pt_E,a4,a6,p);
159 77 : else if (lgefint(p)==3)
160 : {
161 61 : ulong pp = uel(p,2);
162 61 : GEN Tp = ZXT_to_FlxT(T,pp);
163 61 : return get_FlxqE_group(pt_E, Fq_to_Flx(a4, Tp, pp), Fq_to_Flx(a6, Tp, pp),
164 : Tp, pp);
165 : }
166 16 : return get_FpXQE_group(pt_E,a4,a6,T,p);
167 : }
168 :
169 : /***********************************************************************/
170 : /** **/
171 : /** Handle curves with CM by small order **/
172 : /** **/
173 : /***********************************************************************/
174 :
175 : /* l odd prime. Return the list of discriminants D such that
176 : * polclass(D) | poldisc(polmodular(l)) */
177 : static GEN
178 14 : list_singular_discs(long l)
179 : {
180 14 : const long _4l2 = 4*l*l;
181 : long v;
182 14 : GEN V = zero_F2v(_4l2);
183 : /* special cased for efficiency + remove factor l^2 from conductor */
184 14 : F2v_set(V, 4); /* v = 0 */
185 14 : F2v_set(V, 3); /* v = l */
186 1232 : for (v = 1; v < 2*l; v++)
187 1218 : if (v != l)
188 : { /* l does not divide _4l2 - v*v */
189 1204 : GEN F = factoru(_4l2 - v*v), P, E, c;
190 1204 : ulong d = coredisc2u_fact(F, -1, &P, &E);
191 : long i, lc;
192 1204 : c = divisorsu_fact(mkvec2(P,E));
193 1204 : lc = lg(c);
194 3528 : for (i = 1; i < lc; i++)
195 2324 : F2v_set(V, d * uel(c,i)*uel(c,i));
196 : }
197 14 : return V;
198 : }
199 :
200 : /* l odd prime. Find D such that j has CM by D, assuming
201 : * subst(polmodular(l),x,j) has a double root */
202 : static long
203 14 : find_CM(long l, GEN j, GEN T, GEN p)
204 : {
205 14 : const long inv = 0;
206 14 : GEN v = list_singular_discs(l);
207 14 : long i, n = v[1];
208 14 : GEN db = polmodular_db_init(inv);
209 861 : for (i = 1; i < n; i++)
210 861 : if (F2v_coeff(v,i))
211 : {
212 161 : GEN C = polclass0(-i, inv, 0, &db);
213 161 : GEN F = FqX_eval(C, j, T, p);
214 161 : if (signe(F)==0) break;
215 : }
216 14 : gunclone_deep(db); return i < n ? -i: 0;
217 : }
218 :
219 : static GEN
220 14 : vecpoints_to_vecx(GEN x, GEN q1)
221 : {
222 42 : pari_APPLY_type(t_COL, gadd(q1, signe(gmael(x,i,2)) > 0 ? gmael(x,i,1)
223 : : negi(gmael(x,i,1))));
224 : }
225 :
226 : static GEN
227 14 : Fq_ellcard_CM(long disc, GEN a4, GEN a6, GEN T, GEN p)
228 : {
229 : const struct bb_group *grp;
230 : void *E;
231 14 : long d = T ? degpol(T): 1;
232 14 : GEN q = powiu(p, d), q1 = addiu(q, 1), Q, S;
233 14 : Q = qfbsolve(Qfb0(gen_1,gen_0,stoi(-disc)), mkmat22(gen_2, gen_2, p, utoi(d)), 3);
234 14 : if (lg(Q)==1) return q1;
235 14 : S = vecpoints_to_vecx(Q, q1);
236 14 : grp = get_FqE_group(&E, a4, a6, T, p);
237 14 : return gen_select_order(S, E, grp);
238 : }
239 :
240 : /***********************************************************************/
241 : /** **/
242 : /** n-division polynomial **/
243 : /** **/
244 : /***********************************************************************/
245 :
246 : static GEN divpol(GEN t, GEN r2, long n, void *E, const struct bb_algebra *ff);
247 :
248 : /* f_n^2, return ff->(zero|one) or a clone */
249 : static GEN
250 145208 : divpol_f2(GEN t, GEN r2, long n, void *E, const struct bb_algebra *ff)
251 : {
252 145208 : if (n==0) return ff->zero(E);
253 145208 : if (n<=2) return ff->one(E);
254 120428 : if (gmael(t,2,n)) return gmael(t,2,n);
255 44408 : gmael(t,2,n) = ff->sqr(E,divpol(t,r2,n,E,ff));
256 44408 : return gmael(t,2,n);
257 : }
258 :
259 : /* f_n f_{n-2}, return ff->zero or a clone */
260 : static GEN
261 88214 : divpol_ff(GEN t, GEN r2, long n, void *E, const struct bb_algebra *ff)
262 : {
263 88214 : if (n<=2) return ff->zero(E);
264 88214 : if (gmael(t,3,n)) return gmael(t,3,n);
265 56784 : if (n<=4) return divpol(t,r2,n,E,ff);
266 25011 : gmael(t,3,n) = ff->mul(E,divpol(t,r2,n,E,ff), divpol(t,r2,n-2,E,ff));
267 25011 : return gmael(t,3,n);
268 : }
269 :
270 : /* f_n, return ff->zero or a clone */
271 : static GEN
272 188048 : divpol(GEN t, GEN r2, long n, void *E, const struct bb_algebra *ff)
273 : {
274 188048 : long m = n/2;
275 : pari_sp av;
276 : GEN f;
277 188048 : if (n==0) return ff->zero(E);
278 184240 : if (gmael(t,1,n)) return gmael(t,1,n);
279 51331 : switch(n)
280 : {
281 7224 : case 1:
282 : case 2:
283 7224 : f = ff->one(E);
284 7224 : break;
285 44107 : default:
286 44107 : if (odd(n))
287 : {
288 25627 : GEN a = divpol_ff(t,r2,m+2,E,ff);
289 25627 : GEN b = divpol_f2(t,r2,m,E,ff);
290 25627 : GEN c = divpol_ff(t,r2,m+1,E,ff);
291 25627 : GEN d = divpol_f2(t,r2,m+1,E,ff);
292 25627 : av = avma;
293 25627 : if (odd(m))
294 11382 : f = ff->sub(E, ff->mul(E, a,b), ff->mul(E, r2, ff->mul(E, c,d)));
295 : else
296 14245 : f = ff->sub(E, ff->mul(E, r2, ff->mul(E, a,b)), ff->mul(E, c,d));
297 : }
298 : else
299 : {
300 18480 : GEN a = divpol_ff(t,r2,m+2,E,ff);
301 18480 : GEN b = divpol_f2(t,r2,m-1,E,ff);
302 18480 : GEN c = divpol_ff(t,r2,m,E,ff);
303 18480 : GEN d = divpol_f2(t,r2,m+1,E,ff);
304 18480 : av = avma;
305 18480 : f = ff->sub(E, ff->mul(E, a,b), ff->mul(E, c,d));
306 : }
307 44107 : f = gc_upto(av, f);
308 : }
309 51331 : gmael(t,1,n) = f;
310 51331 : return f;
311 : }
312 :
313 : static GEN
314 1522 : Flxq_elldivpol34(long n, GEN a4, GEN a6, GEN S, GEN T, ulong p)
315 : {
316 : GEN res;
317 1522 : long vs = T[1];
318 1522 : switch(n)
319 : {
320 761 : case 3:
321 761 : res = mkpoln(5, Fl_to_Flx(3%p,vs), pol0_Flx(vs), Flx_mulu(a4, 6, p),
322 : Flx_mulu(a6, 12, p), Flx_neg(Flxq_sqr(a4, T, p), p));
323 761 : break;
324 761 : case 4:
325 : {
326 761 : GEN a42 = Flxq_sqr(a4, T, p);
327 1522 : res = mkpoln(7, pol1_Flx(vs), pol0_Flx(vs), Flx_mulu(a4, 5, p),
328 : Flx_mulu(a6, 20, p), Flx_mulu(a42,p-5, p),
329 : Flx_mulu(Flxq_mul(a4, a6, T, p), p-4, p),
330 761 : Flx_sub(Flx_mulu(Flxq_sqr(a6, T, p), p-8%p, p),
331 : Flxq_mul(a4, a42, T, p), p));
332 761 : res = FlxX_double(res, p);
333 : }
334 761 : break;
335 0 : default:
336 0 : pari_err_BUG("Flxq_elldivpol34");
337 : return NULL;/*LCOV_EXCL_LINE*/
338 : }
339 1522 : if(S)
340 : {
341 1522 : setvarn(res, get_FlxqX_var(S));
342 1522 : res = FlxqX_rem(res, S, T, p);
343 : }
344 1522 : return res;
345 : }
346 :
347 : static GEN
348 21102 : Fq_elldivpol34(long n, GEN a4, GEN a6, GEN S, GEN T, GEN p)
349 : {
350 : GEN res;
351 21102 : switch(n)
352 : {
353 10551 : case 3:
354 10551 : res = mkpoln(5, utoi(3), gen_0, Fq_mulu(a4, 6, T, p),
355 : Fq_mulu(a6, 12, T, p), Fq_neg(Fq_sqr(a4, T, p), T, p));
356 10551 : break;
357 10551 : case 4:
358 : {
359 10551 : GEN a42 = Fq_sqr(a4, T, p);
360 10551 : res = mkpoln(7, gen_1, gen_0, Fq_mulu(a4, 5, T, p),
361 : Fq_mulu(a6, 20, T, p), Fq_Fp_mul(a42,stoi(-5), T, p),
362 : Fq_Fp_mul(Fq_mul(a4, a6, T, p), stoi(-4), T, p),
363 : Fq_sub(Fq_Fp_mul(Fq_sqr(a6, T, p), stoi(-8), T, p),
364 : Fq_mul(a4,a42, T, p), T, p));
365 10551 : res = FqX_mulu(res, 2, T, p);
366 : }
367 10551 : break;
368 0 : default:
369 0 : pari_err_BUG("Fq_elldivpol34");
370 : return NULL;/*LCOV_EXCL_LINE*/
371 : }
372 21102 : if (S)
373 : {
374 21102 : setvarn(res, get_FpXQX_var(S));
375 21102 : res = FqX_rem(res, S, T, p);
376 : }
377 21102 : return res;
378 : }
379 :
380 : static GEN
381 17670 : rhs(GEN a4, GEN a6, long v)
382 : {
383 17670 : GEN RHS = mkpoln(4, gen_1, gen_0, a4, a6);
384 17670 : setvarn(RHS, v); return RHS;
385 : }
386 :
387 : static GEN
388 1132 : Flxq_rhs(GEN a4, GEN a6, long v, long vs)
389 : {
390 1132 : GEN RHS = mkpoln(4, pol1_Flx(vs), pol0_Flx(vs), a4, a6);
391 1132 : setvarn(RHS, v); return RHS;
392 : }
393 :
394 : struct divpolmod_red
395 : {
396 : const struct bb_algebra *ff;
397 : void *E;
398 : GEN t, r2;
399 : };
400 :
401 : static void
402 11312 : divpolmod_init(struct divpolmod_red *d, GEN D3, GEN D4, GEN RHS, long n,
403 : void *E, const struct bb_algebra *ff)
404 : {
405 11312 : long k = n+2;
406 11312 : d->ff = ff; d->E = E;
407 11312 : d->t = mkvec3(const_vec(k, NULL),const_vec(k, NULL),const_vec(k, NULL));
408 11312 : if (k>=3) gmael(d->t,1,3) = D3;
409 11312 : if (k>=4) gmael(d->t,1,4) = D4;
410 11312 : d->r2 = ff->sqr(E, RHS);
411 11312 : }
412 :
413 : static void
414 10551 : Fq_elldivpolmod_init(struct divpolmod_red *d, GEN a4, GEN a6, long n, GEN h, GEN T, GEN p)
415 : {
416 : void *E;
417 : const struct bb_algebra *ff;
418 10551 : GEN RHS, D3 = NULL, D4 = NULL;
419 10551 : long v = h ? get_FpXQX_var(h): 0;
420 10551 : D3 = n>=0 ? Fq_elldivpol34(3, a4, a6, h, T, p): NULL;
421 10551 : D4 = n>=1 ? Fq_elldivpol34(4, a4, a6, h, T, p): NULL;
422 10551 : RHS = rhs(a4, a6, v);
423 10551 : RHS = h ? FqX_rem(RHS, h, T, p): RHS;
424 10551 : RHS = FqX_mulu(RHS, 4, T, p);
425 10551 : ff = h ? T ? get_FpXQXQ_algebra(&E, h, T, p): get_FpXQ_algebra(&E, h, p):
426 0 : T ? get_FpXQX_algebra(&E, T, p, v): get_FpX_algebra(&E, p, v);
427 10551 : divpolmod_init(d, D3, D4, RHS, n, E, ff);
428 10551 : }
429 :
430 : static void
431 761 : Flxq_elldivpolmod_init(struct divpolmod_red *d, GEN a4, GEN a6, long n, GEN h, GEN T, ulong p)
432 : {
433 : void *E;
434 : const struct bb_algebra *ff;
435 761 : GEN RHS, D3 = NULL, D4 = NULL;
436 761 : long v = h ? get_FlxqX_var(h) : -1, vT = get_Flx_var(T);
437 761 : D3 = n>=0 ? Flxq_elldivpol34(3, a4, a6, h, T, p): NULL;
438 761 : D4 = n>=1 ? Flxq_elldivpol34(4, a4, a6, h, T, p): NULL;
439 761 : RHS = Flxq_rhs(a4, a6, v, vT);
440 761 : if (h) RHS = FlxqX_rem(RHS, h, T, p);
441 761 : RHS = FlxX_Fl_mul(RHS, 4, p);
442 761 : ff = h ? get_FlxqXQ_algebra(&E, h, T, p) : get_FlxqX_algebra(&E, T, p, 0);
443 761 : divpolmod_init(d, D3, D4, RHS, n, E, ff);
444 761 : }
445 :
446 : /*Computes the n-division polynomial modulo the polynomial h \in Fq[x] */
447 : GEN
448 390 : Flxq_elldivpolmod(GEN a4, GEN a6, long n, GEN h, GEN T, ulong p)
449 : {
450 : struct divpolmod_red d;
451 390 : pari_sp ltop = avma;
452 : GEN res;
453 390 : Flxq_elldivpolmod_init(&d, a4, a6, n, h, T, p);
454 390 : res = divpol(d.t,d.r2,n,d.E,d.ff);
455 390 : return gc_GEN(ltop, res);
456 : }
457 :
458 : /*Computes the n-division polynomial modulo the polynomial h \in Fq[x] */
459 : GEN
460 4851 : Fq_elldivpolmod(GEN a4, GEN a6, long n, GEN h, GEN T, GEN p)
461 : {
462 : struct divpolmod_red d;
463 4851 : pari_sp ltop = avma;
464 : GEN res;
465 4851 : if (lgefint(p)==3 && T)
466 : {
467 390 : ulong pp = p[2];
468 390 : GEN a4p = ZX_to_Flx(a4,pp), a6p = ZX_to_Flx(a6,pp);
469 390 : GEN hp = h ? ZXX_to_FlxX(h, pp, get_FpX_var(T)) : NULL;
470 390 : GEN Tp = ZXT_to_FlxT(T, pp);
471 390 : res = Flxq_elldivpolmod(a4p, a6p, n, hp, Tp, pp);
472 390 : return gc_upto(ltop, FlxX_to_ZXX(res));
473 : }
474 4461 : Fq_elldivpolmod_init(&d, a4, a6, n, h, T, p);
475 4461 : res = divpol(d.t,d.r2,n,d.E,d.ff);
476 4461 : return gc_GEN(ltop, res);
477 : }
478 :
479 : GEN
480 0 : FpXQ_elldivpol(GEN a4, GEN a6, long n, GEN T, GEN p)
481 0 : { return Fq_elldivpolmod(a4,a6,n,NULL,T,p); }
482 :
483 : GEN
484 0 : Fp_elldivpol(GEN a4, GEN a6, long n, GEN p)
485 0 : { return Fq_elldivpolmod(a4,a6,n,NULL,NULL,p); }
486 :
487 : static GEN
488 24451 : Fq_ellyn(struct divpolmod_red *d, long k)
489 : {
490 24451 : void *E = d->E;
491 24451 : const struct bb_algebra *ff = d->ff;
492 24451 : if (k==1) return mkvec2(ff->one(E), ff->one(E));
493 : else
494 : {
495 18998 : GEN t = d->t, r2 = d->r2;
496 18998 : GEN pn2 = divpol(t,r2,k-2,E,ff);
497 18998 : GEN pp2 = divpol(t,r2,k+2,E,ff);
498 18998 : GEN pn12 = divpol_f2(t,r2,k-1,E,ff);
499 18998 : GEN pp12 = divpol_f2(t,r2,k+1,E,ff);
500 18998 : GEN on = ff->red(E,ff->sub(E, ff->mul(E,pp2,pn12), ff->mul(E,pn2,pp12)));
501 18998 : GEN f = divpol(t,r2,k,E,ff);
502 18998 : GEN f2 = divpol_f2(t,r2,k,E,ff);
503 18998 : GEN f3 = ff->mul(E,f,f2);
504 18998 : if (!odd(k)) f3 = ff->mul(E,f3,r2);
505 18998 : return mkvec2(on, f3);
506 : }
507 : }
508 :
509 : static GEN
510 1540 : Fq_elldivpol2(GEN a4, GEN a6, GEN T, GEN p)
511 1540 : { return mkpoln(4, utoi(4), gen_0, Fq_mulu(a4, 4, T, p), Fq_mulu(a6, 4, T, p)); }
512 :
513 : static GEN
514 1540 : Fq_elldivpol2d(GEN a4, GEN T, GEN p)
515 1540 : { return mkpoln(3, utoi(6), gen_0, Fq_mulu(a4, 2, T, p)); }
516 :
517 : static GEN
518 1540 : FqX_numer_isog_abscissa(GEN h, GEN a4, GEN a6, GEN T, GEN p, long vx)
519 : {
520 : GEN mp1, dh, ddh, t, u, t1, t2, t3, t4, f0;
521 1540 : long m = degpol(h);
522 1540 : mp1 = gel(h, m + 1); /* negative of first power sum */
523 1540 : dh = FqX_deriv(h, T, p);
524 1540 : ddh = FqX_deriv(dh, T, p);
525 1540 : t = Fq_elldivpol2(a4, a6, T, p);
526 1540 : u = Fq_elldivpol2d(a4, T, p);
527 1540 : t1 = FqX_sub(FqX_sqr(dh, T, p), FqX_mul(ddh, h, T, p), T, p);
528 1540 : t2 = FqX_mul(u, FqX_mul(h, dh, T, p), T, p);
529 1540 : t3 = FqX_mul(FqX_sqr(h, T, p),
530 : deg1pol_shallow(stoi(2*m), Fq_mulu(mp1, 2, T, p), vx), T, p);
531 1540 : f0 = FqX_add(FqX_sub(FqX_mul(t, t1, T, p), t2, T, p), t3, T, p);
532 1540 : t4 = FqX_mul(pol_x(vx), FqX_sqr(h, T, p), T, p);
533 1540 : return FqX_add(t4, f0, T, p);
534 : }
535 :
536 : static GEN
537 1036 : Zq_inv(GEN b, GEN T, GEN p, long e)
538 : {
539 2023 : return e==1 ? Fq_inv(b, T, p):
540 987 : typ(b)==t_INT ? Zp_inv(b, p, e): ZpXQ_inv(b, T, p, e);
541 : }
542 :
543 : static GEN
544 98441 : Zq_div(GEN a, GEN b, GEN T, GEN q, GEN p, long e)
545 : {
546 98441 : if (e==1) return Fq_div(a, b, T, p);
547 987 : return Fq_mul(a, Zq_inv(b, T, p, e), T, q);
548 : }
549 :
550 : static GEN
551 0 : Zq_sqrt(GEN b, GEN T, GEN p, long e)
552 : {
553 0 : return e==1 ? Fq_sqrt(b, T, p):
554 0 : typ(b)==t_INT ? Zp_sqrt(b, p, e): ZpXQ_sqrt(b, T, p, e);
555 : }
556 :
557 : static GEN
558 14 : Zq_divexact(GEN a, GEN b)
559 14 : { return typ(a)==t_INT ? diviiexact(a, b): ZX_Z_divexact(a, b); }
560 :
561 : static long
562 14 : Zq_pval(GEN a, GEN p)
563 14 : { return typ(a)==t_INT ? Z_pval(a, p): ZX_pval(a, p); }
564 :
565 : static GEN
566 120204 : Zq_divu_safe(GEN a, ulong b, GEN T, GEN q, GEN p, long e)
567 : {
568 : long v, w;
569 120204 : if (e==1) return Fq_div(a, utoi(b), T, q);
570 2611 : v = u_pvalrem(b, p, &b);
571 2611 : if (v > 0)
572 : {
573 14 : if (signe(a)==0) return gen_0;
574 14 : w = Zq_pval(a, p);
575 14 : if (v > w) return NULL;
576 14 : a = Zq_divexact(a, powiu(p,v));
577 : }
578 2611 : return Fq_Fp_mul(a, Fp_inv(utoi(b), q), T, q);
579 : }
580 :
581 : static GEN
582 164381 : FqX_shift(GEN P,long n)
583 164381 : { return RgX_shift_shallow(P, n); }
584 :
585 : static GEN
586 38822 : FqX_mulhigh_i(GEN f, GEN g, long n, GEN T, GEN p)
587 38822 : { return FqX_shift(FqX_mul(f,g,T, p),-n); }
588 :
589 : static GEN
590 38822 : FqX_mulhigh(GEN f, GEN g, long n2, long n, GEN T, GEN p)
591 : {
592 38822 : GEN F = RgX_blocks(f, n2, 2), fl = gel(F,1), fh = gel(F,2);
593 38822 : return FqX_add(FqX_mulhigh_i(fl, g, n2, T, p), FqXn_mul(fh, g, n - n2, T, p), T, p);
594 : }
595 :
596 : static GEN
597 19411 : FqX_invlift1(GEN Q, GEN P, long t1, long t2, GEN T, GEN p)
598 : {
599 19411 : GEN H = FqXn_mul(FqX_mulhigh(Q, P, t1, t2, T, p), Q, t2-t1, T, p);
600 19411 : return FqX_sub(Q, FqX_shift(H, t1), T, p);
601 : }
602 :
603 : static GEN
604 19411 : FqX_invsqrtlift1(GEN Q, GEN P, long t1, long t2, GEN T, GEN p)
605 : {
606 19411 : GEN D = FqX_mulhigh(P, FqX_sqr(Q, T, p), t1, t2, T, p);
607 19411 : GEN H = FqXn_mul(Q, FqX_halve(D, T, p), t2-t1, T, p);
608 19411 : return FqX_sub(Q, FqX_shift(H, t1), T, p);
609 : }
610 :
611 : /* Q(x^2) = intformal(subst(x^N*P,x,x^2)) */
612 : static GEN
613 26537 : ZqX_integ2Xn(GEN P, long N, GEN T, GEN p, GEN pp, long e)
614 : {
615 26537 : long d = degpol(P), v = varn(P);
616 : long k;
617 : GEN Q;
618 26537 : if(d==-1) return pol_0(v);
619 19411 : Q = cgetg(d+3,t_POL);
620 19411 : Q[1] = evalsigne(1) | evalvarn(v);
621 83076 : for (k = 0; k <= d; k++)
622 : {
623 63665 : GEN q = Zq_divu_safe(gel(P,2+k), 2*(k+N)+1, T, p, pp, e);
624 63665 : if (!q) return NULL;
625 63665 : gel(Q, 2+k) = q;
626 : }
627 19411 : return ZXX_renormalize(Q,d+3);
628 : }
629 :
630 : /* solution of G*(S'^2)=(S/x)*(HoS) mod x^m */
631 : static GEN
632 7126 : Zq_Weierstrass(GEN a4, GEN a6, GEN b4, GEN b6, long m, GEN T, GEN p, GEN pp, long n)
633 : {
634 7126 : pari_sp av = avma;
635 7126 : long v = 0;
636 7126 : ulong mask = quadratic_prec_mask(m);
637 7126 : GEN iGdS2 = pol_1(v);
638 7126 : GEN G = mkpoln(4, a6, a4, gen_0, gen_1);
639 7126 : GEN GdS2 = G, S = pol_x(v), sG = pol_1(v), isG = sG, dS = sG;
640 7126 : long N = 1;
641 26537 : for (;mask>1;)
642 : {
643 : GEN S2, HS, K, dK, E;
644 26537 : long N2 = N, d;
645 26537 : N<<=1; if (mask & 1) N--;
646 26537 : mask >>= 1;
647 26537 : d = N-N2;
648 26537 : S2 = FqX_sqr(S, T, p);
649 26537 : HS = FqX_Fq_add(FqX_Fq_mul(S, b6, T, p), b4, T, p);
650 26537 : HS = FqX_Fq_add(FqXn_mul(S2, HS, N, T, p), gen_1, T, p);
651 26537 : HS = FqXn_mul(HS, FqX_shift(S,-1), N, T, p);
652 26537 : sG = FqXn_mul(G, isG, N2, T, p);
653 : /* (HS-Gds2)/(Gds2*sG) */
654 26537 : dK = FqXn_mul(FqX_shift(FqX_sub(HS, GdS2, T, p), -N2),
655 : FqXn_mul(iGdS2, isG, d, T, p), d, T, p);
656 26537 : K = ZqX_integ2Xn(dK, N2, T, p, pp, n);
657 26537 : if (!K) return gc_NULL(av);
658 26537 : E = FqXn_mul(FqXn_mul(K, sG, d, T, p), dS, d, T, p);
659 26537 : S = FqX_add(S, FqX_shift(E, N2+1), T, p);
660 26537 : if (mask <= 1) break;
661 19411 : isG = FqX_invsqrtlift1(isG, G, N2, N, T, p);
662 19411 : dS = FqX_deriv(S, T, p);
663 19411 : GdS2 = FqX_mul(G, FqX_sqr(dS, T, p), T, p);
664 19411 : iGdS2 = FqX_invlift1(iGdS2, GdS2, N2, N, T, p);
665 : }
666 7126 : return gc_upto(av, S);
667 : }
668 :
669 : static GEN
670 7126 : ZqXn_WNewton(GEN S, long l, GEN a4, GEN a6, GEN pp1, GEN T, GEN p, GEN pp, long e)
671 : {
672 7126 : long d = degpol(S);
673 : long k;
674 7126 : GEN Ge = cgetg(2+d,t_POL);
675 7126 : Ge[1] = evalsigne(1);
676 7126 : gel(Ge,2) = pp1;
677 7126 : if (d >= 2)
678 : {
679 7126 : GEN g = Zq_divu_safe(Fq_sub(gel(S,4), Fq_mulu(a4,(l-1),T,p),T,p), 6,T,p,pp,e);
680 7126 : if (!g) return NULL;
681 7126 : gel(Ge, 3) = g;
682 : }
683 7126 : if (d >= 3)
684 : {
685 7126 : GEN g = Zq_divu_safe(Fq_sub(Fq_sub(gel(S,5),
686 : Fq_mul(a4,Fq_mulu(pp1,6,T,p),T,p),T,p),
687 7126 : Fq_mulu(a6,(l-1)*2,T,p),T,p),10,T,p,pp,e);
688 7126 : if (!g) return NULL;
689 7126 : gel(Ge, 4) = g;
690 : }
691 49413 : for (k = 4; k <= d; k++)
692 : {
693 84574 : GEN g = Zq_divu_safe(Fq_sub(Fq_sub(gel(S,4+k-2),
694 42287 : Fq_mul(a4,Fq_mulu(gel(Ge,k-1),4*k-6,T,p),T,p),T,p),
695 42287 : Fq_mul(a6,Fq_mulu(gel(Ge,k-2),4*k-8,T,p),T,p),T,p),
696 42287 : 4*k-2, T, p, pp, e);
697 42287 : if (!g) return NULL;
698 42287 : gel(Ge, k+1) = g;
699 : }
700 7126 : return ZXX_renormalize(Ge, 2+d);
701 : }
702 :
703 : /****************************************************************************/
704 : /* SIMPLE ELLIPTIC CURVE OVER Fq */
705 : /****************************************************************************/
706 :
707 : static GEN
708 2604 : Fq_ellj(GEN a4, GEN a6, GEN T, GEN p)
709 : {
710 2604 : pari_sp ltop=avma;
711 2604 : GEN a43 = Fq_mulu(Fq_powu(a4, 3, T, p), 4, T, p);
712 2604 : GEN j = Fq_div(Fq_mulu(a43, 1728, T, p),
713 : Fq_add(a43, Fq_mulu(Fq_sqr(a6, T, p), 27, T, p), T, p), T, p);
714 2604 : return gc_upto(ltop, j);
715 : }
716 :
717 : static GEN
718 2688 : Zq_ellj(GEN a4, GEN a6, GEN T, GEN p, GEN pp, long e)
719 : {
720 2688 : pari_sp ltop=avma;
721 2688 : GEN a43 = Fq_mulu(Fq_powu(a4, 3, T, p), 4, T, p);
722 2688 : GEN j = Zq_div(Fq_mulu(a43, 1728, T, p),
723 : Fq_add(a43, Fq_mulu(Fq_sqr(a6, T, p), 27, T, p), T, p), T, p, pp, e);
724 2688 : return gc_upto(ltop, j);
725 : }
726 : /****************************************************************************/
727 : /* EIGENVALUE */
728 : /****************************************************************************/
729 :
730 : static GEN
731 371 : Flxq_find_eigen_Frobenius(GEN a4, GEN a6, GEN h, GEN T, ulong p)
732 : {
733 371 : long v = get_FlxqX_var(h), vT = get_Flx_var(T);
734 371 : GEN RHS = FlxqX_rem(Flxq_rhs(a4, a6, v, vT), h, T, p);
735 371 : return FlxqXQ_halfFrobenius(RHS, h, T, p);
736 : }
737 :
738 : static GEN
739 6090 : Fq_find_eigen_Frobenius(GEN a4, GEN a6, GEN h, GEN T, GEN p)
740 : {
741 6090 : long v = T ? get_FpXQX_var(h): get_FpX_var(h);
742 6090 : GEN RHS = FqX_rem(rhs(a4, a6, v), h, T, p);
743 11942 : return T ? FpXQXQ_halfFrobenius(RHS, h, T, p):
744 5852 : FpXQ_pow(RHS, shifti(p, -1), h, p);
745 : }
746 : /*Finds the eigenvalue of the Frobenius given E, ell odd prime, h factor of the
747 : *ell-division polynomial, p and tr the possible values for the trace
748 : *(useful for primes with one root)*/
749 : static ulong
750 504 : find_eigen_value_oneroot(GEN a4, GEN a6, ulong ell, GEN tr, GEN h, GEN T, GEN p)
751 : {
752 504 : pari_sp ltop = avma;
753 : ulong t;
754 : struct divpolmod_red d;
755 : GEN f, Dy, Gy;
756 504 : h = FqX_get_red(h, T, p);
757 504 : Gy = Fq_find_eigen_Frobenius(a4, a6, h, T, p);
758 504 : t = Fl_div(tr[1], 2, ell);
759 504 : if (t < (ell>>1)) t = ell - t;
760 504 : Fq_elldivpolmod_init(&d, a4, a6, t, h, T, p);
761 504 : f = Fq_ellyn(&d, t);
762 504 : Dy = FqXQ_mul(Gy, gel(f,2), h, T, p);
763 504 : if (!gequal(gel(f,1), Dy)) t = ell-t;
764 504 : return gc_ulong(ltop, t);
765 : }
766 :
767 : static ulong
768 371 : Flxq_find_eigen_value_power(GEN a4, GEN a6, ulong ell, long k, ulong lambda,
769 : GEN h, GEN T, ulong p)
770 : {
771 371 : pari_sp ltop = avma;
772 371 : ulong t, ellk1 = upowuu(ell, k-1), ellk = ell*ellk1;
773 : pari_timer ti;
774 : struct divpolmod_red d;
775 : GEN Gy;
776 371 : timer_start(&ti);
777 371 : h = FlxqX_get_red(h, T, p);
778 371 : Gy = Flxq_find_eigen_Frobenius(a4, a6, h, T, p);
779 371 : if (DEBUGLEVEL>2) err_printf(" (%ld ms)",timer_delay(&ti));
780 371 : Flxq_elldivpolmod_init(&d, a4, a6, ellk, h, T, p);
781 1685 : for (t = lambda; t < ellk; t += ellk1)
782 : {
783 1685 : GEN f = Fq_ellyn(&d, t);
784 1685 : GEN Dr = FlxqXQ_mul(Gy, gel(f,2), h, T, p);
785 1685 : if (varn(gel(f,1))!=varn(Dr)) pari_err_BUG("find_eigen_value_power");
786 1685 : if (gequal(gel(f,1), Dr)) break;
787 1441 : if (gequal(gel(f,1), FlxX_neg(Dr,p))) { t = ellk-t; break; }
788 : }
789 371 : if (DEBUGLEVEL>2) err_printf(" (%ld ms)",timer_delay(&ti));
790 371 : return gc_ulong(ltop, t);
791 : }
792 :
793 : /*Finds the eigenvalue of the Frobenius modulo ell^k given E, ell, k, h factor
794 : *of the ell-division polynomial, lambda the previous eigen value and p */
795 : static ulong
796 5586 : Fq_find_eigen_value_power(GEN a4, GEN a6, ulong ell, long k, ulong lambda, GEN h, GEN T, GEN p)
797 : {
798 5586 : pari_sp ltop = avma;
799 5586 : ulong t, ellk1 = upowuu(ell, k-1), ellk = ell*ellk1;
800 : pari_timer ti;
801 : struct divpolmod_red d;
802 : GEN Gy;
803 5586 : timer_start(&ti);
804 5586 : h = FqX_get_red(h, T, p);
805 5586 : Gy = Fq_find_eigen_Frobenius(a4, a6, h, T, p);
806 5586 : if (DEBUGLEVEL>2) err_printf(" (%ld ms)",timer_delay(&ti));
807 5586 : Fq_elldivpolmod_init(&d, a4, a6, ellk, h, T, p);
808 22262 : for (t = lambda; t < ellk; t += ellk1)
809 : {
810 22262 : GEN f = Fq_ellyn(&d, t);
811 22262 : GEN Dr = FqXQ_mul(Gy, gel(f,2), h, T, p);
812 22262 : if (varn(gel(f,1))!=varn(Dr)) pari_err_BUG("find_eigen_value_power");
813 22262 : if (gequal(gel(f,1), Dr)) break;
814 17830 : if (gequal(gel(f,1), FqX_neg(Dr,T,p))) { t = ellk-t; break; }
815 : }
816 5586 : if (DEBUGLEVEL>2) err_printf(" (%ld ms)",timer_delay(&ti));
817 5586 : return gc_ulong(ltop, t);
818 : }
819 :
820 : static ulong
821 5957 : find_eigen_value_power(GEN a4, GEN a6, ulong ell, long k, ulong lambda, GEN hq, GEN T, GEN p)
822 : {
823 5957 : ulong pp = itou_or_0(p);
824 5957 : if (pp && T)
825 : {
826 371 : GEN a4p = ZX_to_Flx(a4, pp);
827 371 : GEN a6p = ZX_to_Flx(a6, pp);
828 371 : GEN hp = ZXXT_to_FlxXT(hq, pp,varn(a4));
829 371 : GEN Tp = ZXT_to_FlxT(T, pp);
830 371 : return Flxq_find_eigen_value_power(a4p, a6p, ell, k, lambda, hp, Tp, pp);
831 : }
832 5586 : return Fq_find_eigen_value_power(a4, a6, ell, k, lambda, hq, T, p);
833 : }
834 :
835 : static GEN
836 8939 : find_kernel(GEN a4, GEN a6, long l, GEN b4, GEN b6, GEN pp1, GEN T, GEN p, GEN pp, long e)
837 : {
838 : GEN Ge, S, Sd;
839 8939 : long d = ((l+1)>>1)+1;
840 8939 : if (l == 3) return deg1pol_shallow(gen_1, Fq_neg(pp1, T, p), 0);
841 7126 : S = Zq_Weierstrass(a4, a6, b4, b6, d + 1, T, p, pp, e);
842 7126 : if (S==NULL) return NULL;
843 7126 : S = FqX_shift(S, -1);
844 7126 : Sd = FqXn_inv(S, d, T, p);
845 7126 : Ge = ZqXn_WNewton(Sd, l, a4, a6, pp1, T, p, pp, e);
846 7126 : if (!Ge) return NULL;
847 7126 : Ge = FqX_neg(Ge, T, p);
848 714 : Ge = T && lgefint(pp)==3 ? ZlXQXn_expint(Ge, d, T, p, pp[2])
849 7304 : : FqXn_expint(Ge, d, T, p);
850 7126 : Ge = RgX_recip(FqX_red(Ge, T, pp));
851 7126 : if (degpol(Ge)==(l-1)>>1) return Ge;
852 1463 : return NULL;
853 : }
854 :
855 : static GEN
856 6573 : compute_u(GEN gprime, GEN Dxxg, GEN DxJg, GEN DJJg, GEN j, GEN pJ, GEN px, ulong q, GEN E4, GEN E6, GEN T, GEN p, GEN pp, long e)
857 : {
858 6573 : pari_sp ltop = avma;
859 6573 : GEN dxxgj = FqX_eval(Dxxg, j, T, p);
860 6573 : GEN dxJgj = FqX_eval(DxJg, j, T, p);
861 6573 : GEN dJJgj = FqX_eval(DJJg, j, T, p);
862 6573 : GEN E42 = Fq_sqr(E4, T, p), E6ovE4 = Zq_div(E6, E4, T, p, pp, e);
863 6573 : GEN a = Fq_mul(gprime, dxxgj, T, p);
864 6573 : GEN b = Fq_mul(Fq_mul(Fq_mulu(j,2*q, T, p), dxJgj, T, p), E6ovE4, T, p);
865 6573 : GEN c = Fq_mul(Zq_div(Fq_sqr(E6ovE4, T, p), gprime, T, p, pp, e), j, T, p);
866 6573 : GEN d = Fq_mul(Fq_mul(c,sqru(q), T, p), Fq_add(pJ, Fq_mul(j, dJJgj, T, p), T, p), T, p);
867 6573 : GEN f = Fq_sub(Fq_div(E6ovE4,utoi(3), T, p),
868 : Zq_div(E42, Fq_mulu(E6,2,T, p), T, p, pp, e), T, p);
869 6573 : GEN g = Fq_sub(Fq_sub(b,a,T,p), d, T, p);
870 6573 : return gc_upto(ltop, Fq_add(Zq_div(g,px,T,p,pp,e), Fq_mulu(f,q,T,p), T, p));
871 : }
872 :
873 : static void
874 8890 : a4a6t(GEN *a4t, GEN *a6t, ulong l, GEN E4t, GEN E6t, GEN T, GEN p)
875 : {
876 8890 : GEN l2 = modii(sqru(l), p), l4 = Fp_sqr(l2, p), l6 = Fp_mul(l4, l2, p);
877 8890 : *a4t = Fq_mul(E4t, Fp_muls(l4, -3, p), T, p);
878 8890 : *a6t = Fq_mul(E6t, Fp_muls(l6, -2, p), T, p);
879 8890 : }
880 : static void
881 49 : a4a6t_from_J(GEN *a4t, GEN *a6t, ulong l, GEN C4t, GEN C6t, GEN T, GEN p)
882 : {
883 49 : GEN l2 = modii(sqru(l), p), l4 = Fp_sqr(l2, p), l6 = Fp_mul(l4, l2, p);
884 49 : GEN v = Fp_inv(stoi(-864), p), u = Fp_mulu(v, 18, p);
885 49 : *a4t = Fq_mul(C4t, Fp_mul(u, l4, p), T, p);
886 49 : *a6t = Fq_mul(C6t, Fp_mul(v, l6, p), T, p);
887 49 : }
888 : /* Finds the isogenous EC, and the sum of the x-coordinates of the points in
889 : * the kernel of the isogeny E -> Eb
890 : * E: elliptic curve, ell: a prime, meqn: Atkin modular equation
891 : * g: root of meqn defining isogenous curve Eb. */
892 : static GEN
893 2576 : find_isogenous_from_Atkin(GEN a4, GEN a6, ulong ell, struct meqn *MEQN, GEN g, GEN T, GEN pp, long e)
894 : {
895 2576 : pari_sp ltop = avma, btop;
896 2576 : GEN meqn = MEQN->eq, meqnx, Dmeqnx, Roots, gprime, u1;
897 2576 : long k, vJ = MEQN->vy;
898 2576 : GEN p = e==1 ? pp: powiu(pp, e);
899 2576 : GEN j = Zq_ellj(a4, a6, T, p, pp, e);
900 2576 : GEN E4 = Fq_div(a4, stoi(-3), T, p);
901 2576 : GEN E6 = Fq_neg(Fq_halve(a6, T, p), T, p);
902 2576 : GEN Dx = RgX_deriv(meqn);
903 2576 : GEN DJ = deriv(meqn, vJ);
904 2576 : GEN Dxg = FpXY_Fq_evaly(Dx, g, T, p, vJ);
905 2576 : GEN px = FqX_eval(Dxg, j, T, p), dx = Fq_mul(px, g, T, p);
906 2576 : GEN DJg = FpXY_Fq_evaly(DJ, g, T, p, vJ);
907 2576 : GEN pJ = FqX_eval(DJg, j, T, p), dJ = Fq_mul(pJ, j, T, p);
908 2576 : GEN Dxx = RgX_deriv(Dx);
909 2576 : GEN DxJg = FqX_deriv(Dxg, T, p);
910 :
911 2576 : GEN Dxxg = FpXY_Fq_evaly(Dxx, g, T, p, vJ);
912 2576 : GEN DJJg = FqX_deriv(DJg, T, p);
913 : GEN a, b;
914 2576 : if (!signe(Fq_red(dJ,T,pp)) || !signe(Fq_red(dx,T,pp)))
915 : {
916 21 : if (DEBUGLEVEL>0) err_printf("[A: d%c=0]",signe(dJ)? 'x': 'J');
917 21 : return gc_NULL(ltop);
918 : }
919 2555 : a = Fq_mul(dJ, Fq_mul(g, E6, T, p), T, p);
920 2555 : b = Fq_mul(E4, dx, T, p);
921 2555 : gprime = Zq_div(a, b, T, p, pp, e);
922 :
923 2555 : u1 = compute_u(gprime, Dxxg, DxJg, DJJg, j, pJ, px, 1, E4, E6, T, p, pp, e);
924 2555 : meqnx = FpXY_Fq_evaly(meqn, g, T, p, vJ);
925 2555 : Dmeqnx = FqX_deriv(meqnx, T, pp);
926 2555 : Roots = FqX_roots(meqnx, T, pp);
927 :
928 2555 : btop = avma;
929 4032 : for (k = lg(Roots)-1; k >= 1; k--, set_avma(btop))
930 : {
931 4032 : GEN jt = gel(Roots, k);
932 4032 : if (signe(FqX_eval(Dmeqnx, jt, T, pp))==0)
933 0 : continue;
934 4032 : if (e > 1)
935 91 : jt = ZqX_liftroot(meqnx, gel(Roots, k), T, pp, e);
936 4032 : if (signe(Fq_red(jt, T, pp)) == 0 || signe(Fq_sub(jt, utoi(1728), T, pp)) == 0)
937 : {
938 14 : if (DEBUGLEVEL>0) err_printf("[A: jt=%ld]",signe(Fq_red(jt,T,p))? 1728: 0);
939 14 : return gc_NULL(ltop);
940 : }
941 : else
942 : {
943 4018 : GEN pxstar = FqX_eval(Dxg, jt, T, p);
944 4018 : GEN dxstar = Fq_mul(pxstar, g, T, p);
945 4018 : GEN pJstar = FqX_eval(DJg, jt, T, p);
946 4018 : GEN dJstar = Fq_mul(Fq_mulu(jt, ell, T, p), pJstar, T, p);
947 4018 : GEN u = Fq_mul(Fq_mul(dxstar, dJ, T, p), E6, T, p);
948 4018 : GEN v = Fq_mul(Fq_mul(dJstar, dx, T, p), E4, T, p);
949 4018 : GEN E4t = Zq_div(Fq_mul(Fq_sqr(u, T, p), jt, T, p), Fq_mul(Fq_sqr(v, T, p), Fq_sub(jt, utoi(1728), T, p), T, p), T, p, pp, e);
950 4018 : GEN E6t = Zq_div(Fq_mul(u, E4t, T, p), v, T, p, pp, e);
951 4018 : GEN u2 = compute_u(gprime, Dxxg, DxJg, DJJg, jt, pJstar, pxstar, ell, E4t, E6t, T, p, pp, e);
952 4018 : GEN pp1 = Fq_mulu(Fq_sub(u1, u2, T, p), 3*ell, T, p);
953 : GEN a4t, a6t, h;
954 4018 : a4a6t(&a4t, &a6t, ell, E4t, E6t, T, p);
955 4018 : h = find_kernel(a4, a6, ell, a4t, a6t, pp1, T, p, pp, e);
956 4018 : if (h && signe(Fq_elldivpolmod(a4, a6, ell, h, T, pp))==0)
957 2541 : return gc_GEN(ltop, mkvec3(a4t, a6t, h));
958 : }
959 : }
960 0 : pari_err_BUG("find_isogenous_from_Atkin, kernel not found");
961 : return NULL;/*LCOV_EXCL_LINE*/
962 : }
963 :
964 : /* Finds E' ell-isogenous to E and the trace term p1 from canonical modular
965 : * equation meqn
966 : * E: elliptic curve, ell: a prime, meqn: canonical modular equation
967 : * g: root of meqn defining isogenous curve Eb. */
968 : static GEN
969 4879 : find_isogenous_from_canonical(GEN a4, GEN a6, ulong ell, struct meqn *MEQN, GEN g, GEN T, GEN pp, long e)
970 : {
971 4879 : pari_sp ltop = avma;
972 4879 : GEN meqn = MEQN->eq;
973 4879 : long vJ = MEQN->vy;
974 4879 : GEN p = e==1 ? pp: powiu(pp, e);
975 : GEN h;
976 4879 : GEN E4 = Fq_div(a4, stoi(-3), T, p);
977 4879 : GEN E6 = Fq_neg(Fq_halve(a6, T, p), T, p);
978 4879 : GEN E42 = Fq_sqr(E4, T, p);
979 4879 : GEN E43 = Fq_mul(E4, E42, T, p);
980 4879 : GEN E62 = Fq_sqr(E6, T, p);
981 4879 : GEN delta = Fq_div(Fq_sub(E43, E62, T, p), utoi(1728), T, p);
982 4879 : GEN j = Zq_div(E43, delta, T, p, pp, e);
983 4879 : GEN Dx = RgX_deriv(meqn);
984 4879 : GEN DJ = deriv(meqn, vJ);
985 4879 : GEN Dxg = FpXY_Fq_evaly(Dx, g, T, p, vJ);
986 4879 : GEN px = FqX_eval(Dxg, j, T, p), dx = Fq_mul(px, g, T, p);
987 4879 : GEN DJg = FpXY_Fq_evaly(DJ, g, T, p, vJ);
988 4879 : GEN pJ = FqX_eval(DJg, j, T, p), dJ = Fq_mul(j, pJ, T, p);
989 4879 : GEN Dxx = RgX_deriv(Dx);
990 4879 : GEN DxJg = FqX_deriv(Dxg, T, p);
991 :
992 4879 : GEN ExJ = FqX_eval(DxJg, j, T, p);
993 4879 : ulong tis = ugcd(12, ell-1), is = 12 / tis;
994 4879 : GEN itis = Fq_inv(stoi(-tis), T, p);
995 4879 : GEN deltal = Fq_div(Fq_mul(delta, Fq_powu(g, tis, T, p), T, p), powuu(ell, 12), T, p);
996 : GEN E4l, E6l, a4t, a6t, p_1;
997 4879 : if (signe(Fq_red(dx,T, pp))==0)
998 : {
999 0 : if (DEBUGLEVEL>0) err_printf("[C: dx=0]");
1000 0 : return gc_NULL(ltop);
1001 : }
1002 4879 : if (signe(Fq_red(dJ, T, pp))==0)
1003 : {
1004 : GEN jl;
1005 0 : if (DEBUGLEVEL>0) err_printf("[C: dJ=0]");
1006 0 : E4l = Fq_div(E4, sqru(ell), T, p);
1007 0 : jl = Zq_div(Fq_powu(E4l, 3, T, p), deltal, T, p, pp, e);
1008 0 : E6l = Zq_sqrt(Fq_mul(Fq_sub(jl, utoi(1728), T, p),
1009 : deltal, T, p), T, pp, e);
1010 0 : p_1 = gen_0;
1011 : }
1012 : else
1013 : {
1014 : GEN jl, f, fd, Dgs, Djs, jld;
1015 4879 : GEN E2s = Zq_div(Fq_mul(Fq_neg(Fq_mulu(E6, 12, T, p), T, p), dJ, T, p),
1016 : Fq_mul(Fq_mulu(E4, is, T, p), dx, T, p), T, p, pp, e);
1017 4879 : GEN gd = Fq_mul(Fq_mul(E2s, itis, T, p), g, T, p);
1018 4879 : GEN jd = Zq_div(Fq_mul(Fq_neg(E42, T, p), E6, T, p), delta, T, p, pp, e);
1019 4879 : GEN E0b = Zq_div(E6, Fq_mul(E4, E2s, T, p), T, p, pp, e);
1020 4879 : GEN Dxxgj = FqXY_eval(Dxx, g, j, T, p);
1021 4879 : GEN Dgd = Fq_add(Fq_mul(gd, px, T, p), Fq_mul(g, Fq_add(Fq_mul(gd, Dxxgj, T, p), Fq_mul(jd, ExJ, T, p), T, p), T, p), T, p);
1022 4879 : GEN DJgJj = FqX_eval(FqX_deriv(DJg, T, p), j, T, p);
1023 4879 : GEN Djd = Fq_add(Fq_mul(jd, pJ, T, p), Fq_mul(j, Fq_add(Fq_mul(jd, DJgJj, T, p), Fq_mul(gd, ExJ, T, p), T, p), T, p), T, p);
1024 4879 : GEN E0bd = Zq_div(Fq_sub(Fq_mul(Dgd, itis, T, p), Fq_mul(E0b, Djd, T, p), T, p), dJ, T, p, pp, e);
1025 4879 : E4l = Fq_div(Fq_sub(E4, Fq_mul(E2s, Fq_sub(Fq_sub(Fq_add(Zq_div(Fq_mulu(E0bd, 12, T, p), E0b, T, p, pp, e), Zq_div(Fq_mulu(E42, 6, T, p), E6, T, p, pp, e), T, p), Zq_div(Fq_mulu(E6, 4, T, p), E4, T, p, pp, e), T, p), E2s, T, p), T, p), T, p), sqru(ell), T, p);
1026 4879 : jl = Zq_div(Fq_powu(E4l, 3, T, p), deltal, T, p, pp, e);
1027 4879 : if (signe(Fq_red(jl,T,pp))==0)
1028 : {
1029 7 : if (DEBUGLEVEL>0) err_printf("[C: jl=0]");
1030 7 : return gc_NULL(ltop);
1031 : }
1032 4872 : f = Zq_div(powuu(ell, is), g, T, p, pp, e);
1033 4872 : fd = Fq_neg(Fq_mul(Fq_mul(E2s, f, T, p), itis, T, p), T, p);
1034 4872 : Dgs = FqXY_eval(Dx, f, jl, T, p);
1035 4872 : Djs = FqXY_eval(DJ, f, jl, T, p);
1036 4872 : jld = Zq_div(Fq_mul(Fq_neg(fd, T, p), Dgs, T, p),
1037 : Fq_mulu(Djs, ell, T, p), T, p, pp, e);
1038 4872 : E6l = Zq_div(Fq_mul(Fq_neg(E4l, T, p), jld, T, p), jl, T, p, pp, e);
1039 4872 : p_1 = Fq_neg(Fq_halve(Fq_mulu(E2s, ell, T, p), T, p),T,p);
1040 : }
1041 4872 : a4a6t(&a4t, &a6t, ell, E4l, E6l, T, p);
1042 4872 : h = find_kernel(a4, a6, ell, a4t, a6t, p_1, T, p, pp, e);
1043 4872 : if (!h) return NULL;
1044 4872 : return gc_GEN(ltop, mkvec3(a4t, a6t, h));
1045 : }
1046 :
1047 : static GEN
1048 98 : corr(GEN c4, GEN c6, GEN T, GEN p, GEN pp, long e)
1049 : {
1050 98 : GEN c46 = Zq_div(Fq_sqr(c4, T, p), c6, T, p, pp, e);
1051 98 : GEN c64 = Zq_div(c6, c4, T, p, pp, e);
1052 98 : GEN a = Fp_divu(gen_2, 3, p);
1053 98 : return Fq_add(Fq_halve(c46, T, p), Fq_mul(a, c64, T, p), T, p);
1054 : }
1055 :
1056 : static GEN
1057 168 : RgXY_deflatex(GEN H, long n, long d)
1058 : {
1059 168 : long i, l = lg(H);
1060 168 : GEN R = cgetg(l, t_POL);
1061 168 : R[1] = H[1];
1062 980 : for(i = 2; i < l; i++)
1063 : {
1064 812 : GEN Hi = gel(H, i);
1065 812 : gel(R,i) = typ(Hi)==t_POL? RgX_deflate(RgX_shift_shallow(Hi, d), n): Hi;
1066 : }
1067 168 : return RgX_renormalize_lg(R, l);
1068 : }
1069 :
1070 : static GEN
1071 70 : Fq_polmodular_eval(GEN meqn, GEN j, long N, GEN T, GEN p, long vJ)
1072 : {
1073 70 : pari_sp av = avma;
1074 : GEN R, dR, ddR;
1075 70 : long t0 = N%3 == 1 ? 2: 0;
1076 70 : long t2 = N%3 == 1 ? 0: 2;
1077 70 : if (N == 3)
1078 : {
1079 14 : GEN P = FpXX_red(meqn, p);
1080 14 : GEN dP = deriv(P, -1), ddP = deriv(dP, -1);
1081 14 : R = FpXY_Fq_evaly(P, j, T, p, vJ);
1082 14 : dR = FpXY_Fq_evaly(dP, j, T, p, vJ);
1083 14 : ddR = FpXY_Fq_evaly(ddP, j, T, p, vJ);
1084 14 : return gc_GEN(av, mkvec3(R,dR,ddR));
1085 : }
1086 : else
1087 : {
1088 56 : GEN P5 = FpXX_red(meqn, p);
1089 56 : GEN H = RgX_splitting(P5, 3);
1090 56 : GEN H0 = RgXY_deflatex(gel(H,1), 3, -t0);
1091 56 : GEN H1 = RgXY_deflatex(gel(H,2), 3, -1);
1092 56 : GEN H2 = RgXY_deflatex(gel(H,3), 3, -t2);
1093 56 : GEN h0 = FpXY_Fq_evaly(H0, j, T, p, vJ);
1094 56 : GEN h1 = FpXY_Fq_evaly(H1, j, T, p, vJ);
1095 56 : GEN h2 = FpXY_Fq_evaly(H2, j, T, p, vJ);
1096 56 : GEN dH0 = RgX_deriv(H0);
1097 56 : GEN dH1 = RgX_deriv(H1);
1098 56 : GEN dH2 = RgX_deriv(H2);
1099 56 : GEN ddH0 = RgX_deriv(dH0);
1100 56 : GEN ddH1 = RgX_deriv(dH1);
1101 56 : GEN ddH2 = RgX_deriv(dH2);
1102 56 : GEN d0 = FpXY_Fq_evaly(dH0, j, T, p, vJ);
1103 56 : GEN d1 = FpXY_Fq_evaly(dH1, j, T, p, vJ);
1104 56 : GEN d2 = FpXY_Fq_evaly(dH2, j, T, p, vJ);
1105 56 : GEN dd0 = FpXY_Fq_evaly(ddH0, j, T, p, vJ);
1106 56 : GEN dd1 = FpXY_Fq_evaly(ddH1, j, T, p, vJ);
1107 56 : GEN dd2 = FpXY_Fq_evaly(ddH2, j, T, p, vJ);
1108 : GEN h02, h12, h22, h03, h13, h23, h012, dh03, dh13, dh23, dh012;
1109 : GEN ddh03, ddh13, ddh23, ddh012;
1110 : GEN R1, dR1, ddR1, ddR2;
1111 56 : h02 = FqX_sqr(h0, T, p);
1112 56 : h12 = FqX_sqr(h1, T, p);
1113 56 : h22 = FqX_sqr(h2, T, p);
1114 56 : h03 = FqX_mul(h0, h02, T, p);
1115 56 : h13 = FqX_mul(h1, h12, T, p);
1116 56 : h23 = FqX_mul(h2, h22, T, p);
1117 56 : h012 = FqX_mul(FqX_mul(h0, h1, T, p), h2, T, p);
1118 56 : dh03 = FqX_mul(FqX_mulu(d0, 3, T, p), h02, T, p);
1119 56 : dh13 = FqX_mul(FqX_mulu(d1, 3, T, p), h12, T, p);
1120 56 : dh23 = FqX_mul(FqX_mulu(d2, 3, T, p), h22, T, p);
1121 56 : dh012 = FqX_add(FqX_add(FqX_mul(FqX_mul(d0, h1, T, p), h2, T, p), FqX_mul(FqX_mul(h0, d1, T, p), h2, T, p), T, p), FqX_mul(FqX_mul(h0, h1, T, p), d2, T, p), T, p);
1122 56 : R1 = FqX_sub(h13, FqX_mulu(h012, 3, T, p), T, p);
1123 56 : R = FqX_add(FqX_add(FqX_Fq_mul(RgX_shift_shallow(h23, t2), Fq_sqr(j, T, p), T, p), FqX_Fq_mul(RgX_shift_shallow(R1, 1), j, T, p), T, p), RgX_shift_shallow(h03, t0), T, p);
1124 56 : dR1 = FqX_sub(dh13, FqX_mulu(dh012, 3, T, p), T, p);
1125 56 : dR = FqX_add(FqX_add(RgX_shift_shallow(FqX_add(FqX_Fq_mul(dh23, Fq_sqr(j, T, p), T, p), FqX_Fq_mul(h23, Fq_mulu(j, 2, T, p), T, p), T, p), t2), RgX_shift_shallow(FqX_add(FqX_Fq_mul(dR1, j, T, p), R1, T, p), 1), T, p), RgX_shift_shallow(dh03, t0), T, p);
1126 56 : ddh03 = FqX_mulu(FqX_add(FqX_mul(dd0, h02, T, p), FqX_mul(FqX_mulu(FqX_sqr(d0, T, p), 2, T, p), h0, T, p), T, p), 3, T, p);
1127 56 : ddh13 = FqX_mulu(FqX_add(FqX_mul(dd1, h12, T, p), FqX_mul(FqX_mulu(FqX_sqr(d1, T, p), 2, T, p), h1, T, p), T, p), 3, T, p);
1128 56 : ddh23 = FqX_mulu(FqX_add(FqX_mul(dd2, h22, T, p), FqX_mul(FqX_mulu(FqX_sqr(d2, T, p), 2, T, p), h2, T, p), T, p), 3, T, p);
1129 56 : ddh012 = FqX_add(FqX_add(FqX_add(FqX_mul(FqX_mul(dd0, h1, T, p), h2, T, p), FqX_mul(FqX_mul(h0, dd1, T, p), h2, T, p), T, p), FqX_mul(FqX_mul(h0, h1, T, p), dd2, T, p), T, p), FqX_mulu(FqX_add(FqX_add(FqX_mul(FqX_mul(d0, d1, T, p), h2, T, p), FqX_mul(FqX_mul(d0, h1, T, p), d2, T, p), T, p), FqX_mul(FqX_mul(h0, d1, T, p), d2, T, p), T, p), 2, T, p), T, p);
1130 56 : ddR1 = FqX_sub(ddh13, FqX_mulu(ddh012, 3, T, p), T, p);
1131 56 : ddR2 = FqX_add(FqX_add(FqX_Fq_mul(ddh23, Fq_sqr(j, T, p), T, p), FqX_Fq_mul(dh23, Fq_mulu(j, 4, T, p), T, p), T, p), FqX_mulu(h23, 2, T, p), T, p);
1132 56 : ddR = FqX_add(FqX_add(RgX_shift_shallow(ddR2, t2), RgX_shift_shallow(FqX_add(FqX_mulu(dR1, 2, T, p), FqX_Fq_mul(ddR1, j, T, p), T, p), 1), T, p), RgX_shift_shallow(ddh03, t0), T, p);
1133 56 : return gc_GEN(av, mkvec3(R, dR, ddR));
1134 : }
1135 : }
1136 :
1137 : static GEN
1138 11606 : meqn_j(struct meqn *MEQN, GEN j, long ell, GEN T, GEN p)
1139 : {
1140 11606 : if (MEQN->type=='J')
1141 : {
1142 70 : MEQN->eval = Fq_polmodular_eval(MEQN->eq, j, ell, T, p, MEQN->vy);
1143 70 : return gel(MEQN->eval, 1);
1144 : }
1145 : else
1146 11536 : return FqXY_evalx(MEQN->eq, j, T, p);
1147 : }
1148 :
1149 : static GEN
1150 49 : find_isogenous_from_J(GEN a4, GEN a6, ulong ell, struct meqn *MEQN, GEN g, GEN T, GEN pp, long e)
1151 : {
1152 49 : pari_sp ltop = avma;
1153 49 : GEN meqn = MEQN->eval;
1154 49 : GEN p = e==1 ? pp: powiu(pp, e);
1155 : GEN h, a4t, a6t;
1156 : GEN C4, C6, C4t, C6t;
1157 : GEN j, jp, jtp, jtp2, jtp3;
1158 : GEN Py, Pxy, Pyy, Pxj, Pyj, Pxxj, Pxyj, Pyyj;
1159 : GEN s0, s1, s2, s3;
1160 : GEN den, D, co, cot, c0, p_1;
1161 49 : if (signe(g) == 0 || signe(Fq_sub(g, utoi(1728), T, p)) == 0)
1162 : {
1163 0 : if (DEBUGLEVEL>0) err_printf("[J: g=%ld]",signe(g)==0 ?0: 1728);
1164 0 : return gc_NULL(ltop);
1165 : }
1166 49 : C4 = Fq_mul(a4, stoi(-48), T, p);
1167 49 : C6 = Fq_mul(a6, stoi(-864), T, p);
1168 49 : if (signe(C4)==0 || signe(C6)==0)
1169 : {
1170 0 : if (DEBUGLEVEL>0) err_printf("[J: C%ld=0]",signe(C4)==0 ?4: 6);
1171 0 : return gc_NULL(ltop);
1172 : }
1173 49 : j = Zq_ellj(a4, a6, T, p, pp, e);
1174 49 : jp = Fq_mul(j, Zq_div(C6, C4, T, p, pp, e), T, p);
1175 49 : co = corr(C4, C6, T, p, pp, e);
1176 49 : Py = RgX_deriv(gel(meqn, 1));
1177 49 : Pxy = RgX_deriv(gel(meqn,2));
1178 49 : Pyy = RgX_deriv(Py);
1179 49 : Pxj = FqX_eval(gel(meqn, 2), g, T, p);
1180 49 : if (signe(Pxj)==0)
1181 : {
1182 0 : if (DEBUGLEVEL>0) err_printf("[J: Pxj=0]");
1183 0 : return gc_NULL(ltop);
1184 : }
1185 49 : Pyj = FqX_eval(Py, g, T, p);
1186 49 : Pxxj = FqX_eval(gel(meqn, 3), g, T, p);
1187 49 : Pxyj = FqX_eval(Pxy, g, T, p);
1188 49 : Pyyj = FqX_eval(Pyy, g, T, p);
1189 49 : jtp = Fq_div(Fq_mul(jp, Zq_div(Pxj, Pyj, T, p, pp, e), T, p),
1190 : utoineg(ell), T, p);
1191 49 : jtp2 = Fq_sqr(jtp,T,p);
1192 49 : jtp3 = Fq_mul(jtp,jtp2,T,p);
1193 49 : den = Fq_mul(Fq_sqr(g,T,p),Fq_sub(g,utoi(1728),T,p),T, p);
1194 49 : D = Zq_inv(den, T, pp, e);
1195 49 : C4t = Fq_mul(jtp2,Fq_mul(g, D, T, p), T, p);
1196 49 : C6t = Fq_mul(jtp3, D, T, p);
1197 49 : s0 = Fq_mul(Fq_sqr(jp, T, p), Pxxj, T, p);
1198 49 : s1 = Fq_mul(Fq_mulu(Fq_mul(jp,jtp,T,p),2*ell,T,p), Pxyj, T, p);
1199 49 : s2 = Fq_mul(Fq_mulu(jtp2,ell*ell,T,p), Pyyj, T, p);
1200 49 : s3 = Zq_div(Fq_add(s0, Fq_add(s1, s2, T, p), T, p),Fq_mul(jp, Pxj, T, p),T,p,pp,e);
1201 49 : cot = corr(C4t, C6t, T, p, pp, e);
1202 49 : c0 = Fq_sub(co,Fq_mulu(cot,ell,T,p),T,p);
1203 49 : p_1 = Fq_div(Fq_mulu(Fq_add(s3, c0, T, p),ell,T,p),stoi(-4),T,p);
1204 49 : a4a6t_from_J(&a4t, &a6t, ell, C4t, C6t, T, p);
1205 49 : h = find_kernel(a4, a6, ell, a4t, a6t, p_1, T, p, pp, e);
1206 49 : if (!h) return NULL;
1207 49 : return gc_GEN(ltop, mkvec3(a4t, a6t, h));
1208 : }
1209 :
1210 : static GEN
1211 7511 : find_isogenous(GEN a4,GEN a6, ulong ell, struct meqn *MEQN, GEN g, GEN T,GEN p)
1212 : {
1213 7511 : ulong pp = itou_or_0(p);
1214 7511 : long e = pp ? ulogint(((ell+1)>>1)+1, pp) + ulogint(2*ell+4, pp) + 1: 1;
1215 7511 : if (signe(a4)==0 || signe(a6)==0)
1216 : {
1217 7 : if (DEBUGLEVEL>0) err_printf("[%c: j=%ld]",MEQN->type,signe(a4)==0 ?0: 1728);
1218 7 : return NULL;
1219 : }
1220 7504 : if (e > 1)
1221 : {
1222 63 : GEN pe = powiu(p, e);
1223 63 : GEN meqnj = meqn_j(MEQN, Zq_ellj(a4, a6, T, pe, p, e), ell, T, pe);
1224 63 : g = ZqX_liftroot(meqnj, g, T, p, e);
1225 : }
1226 7504 : switch(MEQN->type)
1227 : {
1228 4879 : case 'C': return find_isogenous_from_canonical(a4,a6,ell, MEQN, g, T,p,e);
1229 2576 : case 'A': return find_isogenous_from_Atkin(a4,a6,ell, MEQN, g, T,p,e);
1230 49 : default: return find_isogenous_from_J(a4,a6,ell, MEQN, g, T,p,e);
1231 : }
1232 : }
1233 :
1234 : static GEN
1235 6181 : FqX_homogenous_eval(GEN P, GEN A, GEN B, GEN T, GEN p)
1236 : {
1237 6181 : long d = degpol(P), i, v = varn(A);
1238 6181 : GEN s = scalar_ZX_shallow(gel(P, d+2), v), Bn = pol_1(v);
1239 20454 : for (i = d-1; i >= 0; i--)
1240 : {
1241 14273 : Bn = FqX_mul(Bn, B, T, p);
1242 14273 : s = FqX_add(FqX_mul(s, A, T, p), FqX_Fq_mul(Bn, gel(P,i+2), T, p), T, p);
1243 : }
1244 6181 : return s;
1245 : }
1246 :
1247 : static GEN
1248 1295 : FqX_homogenous_div(GEN P, GEN Q, GEN A, GEN B, GEN T, GEN p)
1249 : {
1250 1295 : GEN z = cgetg(3, t_RFRAC);
1251 1295 : long d = degpol(Q)-degpol(P);
1252 1295 : gel(z, 1) = FqX_homogenous_eval(P, A, B, T, p);
1253 1295 : gel(z, 2) = FqX_homogenous_eval(Q, A, B, T, p);
1254 1295 : if (d > 0)
1255 0 : gel(z, 1) = FqX_mul(gel(z, 1), FqX_powu(B, d, T, p), T, p);
1256 1295 : else if (d < 0)
1257 1295 : gel(z, 2) = FqX_mul(gel(z, 2), FqX_powu(B, -d, T, p), T, p);
1258 1295 : return z;
1259 : }
1260 :
1261 : static GEN
1262 1540 : find_kernel_power(GEN Eba4, GEN Eba6, GEN Eca4, GEN Eca6, ulong ell, struct meqn *MEQN, GEN kpoly, GEN Ib, GEN T, GEN p)
1263 : {
1264 1540 : pari_sp ltop = avma, btop;
1265 : GEN a4t, a6t, gtmp;
1266 1540 : GEN num_iso = FqX_numer_isog_abscissa(kpoly, Eba4, Eba6, T, p, 0);
1267 1540 : GEN mpoly = meqn_j(MEQN, Fq_ellj(Eca4, Eca6, T, p), ell, T, p);
1268 1540 : GEN mroots = FqX_roots(mpoly, T, p);
1269 1540 : GEN kpoly2 = FqX_sqr(kpoly, T, p);
1270 1540 : long i, l1 = lg(mroots);
1271 1540 : btop = avma;
1272 2541 : for (i = 1; i < l1; i++)
1273 : {
1274 : GEN h;
1275 2303 : GEN tmp = find_isogenous(Eca4, Eca6, ell, MEQN, gel(mroots, i), T, p);
1276 2303 : if (!tmp) return gc_NULL(ltop);
1277 2296 : a4t = gel(tmp, 1);
1278 2296 : a6t = gel(tmp, 2);
1279 2296 : gtmp = gel(tmp, 3);
1280 :
1281 : /*check that the kernel kpoly is the good one */
1282 2296 : h = FqX_homogenous_eval(gtmp, num_iso, kpoly2, T, p);
1283 2296 : if (signe(Fq_elldivpolmod(Eba4, Eba6, ell, h, T, p)))
1284 : {
1285 1295 : GEN Ic = FqX_homogenous_div(num_iso,kpoly2, numer_i(Ib),denom_i(Ib), T,p);
1286 1295 : GEN kpoly_new = FqX_homogenous_eval(gtmp, numer_i(Ic),denom_i(Ic), T,p);
1287 1295 : return gc_GEN(ltop, mkvecn(5, a4t, a6t, kpoly_new, gtmp, Ic));
1288 : }
1289 1001 : set_avma(btop);
1290 : }
1291 238 : return gc_NULL(ltop);
1292 : }
1293 :
1294 : /****************************************************************************/
1295 : /* TRACE */
1296 : /****************************************************************************/
1297 : enum mod_type {MTcm, MTpathological, MTAtkin, MTElkies, MTone_root, MTroots};
1298 :
1299 : static GEN
1300 678 : Flxq_study_eqn(GEN mpoly, GEN T, ulong p, long *pt_dG, long *pt_r)
1301 : {
1302 678 : GEN Xq = FlxqX_Frobenius(mpoly, T, p);
1303 678 : GEN G = FlxqX_gcd(FlxX_sub(Xq, pol_x(0), p), mpoly, T, p);
1304 678 : *pt_dG = degpol(G);
1305 678 : if (!*pt_dG) { *pt_r = FlxqX_ddf_degree(mpoly, Xq, T, p); return NULL; }
1306 410 : return gel(FlxqX_roots(G, T, p), 1);
1307 : }
1308 :
1309 : static GEN
1310 8988 : Fp_study_eqn(GEN mpoly, GEN p, long *pt_dG, long *pt_r)
1311 : {
1312 8988 : GEN T = FpX_get_red(mpoly, p);
1313 8988 : GEN XP = FpX_Frobenius(T, p);
1314 8988 : GEN G = FpX_gcd(FpX_sub(XP, pol_x(0), p), mpoly, p);
1315 8988 : *pt_dG = degpol(G);
1316 8988 : if (!*pt_dG) { *pt_r = FpX_ddf_degree(T, XP, p); return NULL; }
1317 4732 : return FpX_oneroot(G, p);
1318 : }
1319 :
1320 : static GEN
1321 9989 : Fq_study_eqn(GEN mpoly, GEN T, GEN p, long *pt_dG, long *pt_r)
1322 : {
1323 : GEN G;
1324 9989 : if (!T) return Fp_study_eqn(mpoly, p, pt_dG, pt_r);
1325 1001 : if (lgefint(p)==3)
1326 : {
1327 678 : ulong pp = p[2];
1328 678 : GEN Tp = ZXT_to_FlxT(T,pp);
1329 678 : GEN mpolyp = ZXX_to_FlxX(mpoly,pp,get_FpX_var(T));
1330 678 : G = Flxq_study_eqn(mpolyp, Tp, pp, pt_dG, pt_r);
1331 678 : return G ? Flx_to_ZX(G): NULL;
1332 : }
1333 : else
1334 : {
1335 323 : GEN Xq = FpXQX_Frobenius(mpoly, T, p);
1336 323 : G = FpXQX_gcd(FpXX_sub(Xq, pol_x(0), p), mpoly, T, p);
1337 323 : *pt_dG = degpol(G);
1338 323 : if (!*pt_dG) { *pt_r = FpXQX_ddf_degree(mpoly, Xq, T, p); return NULL; }
1339 136 : return gel(FpXQX_roots(G, T, p), 1);
1340 : }
1341 : }
1342 :
1343 : /* Berlekamp variant */
1344 : static GEN
1345 10003 : study_modular_eqn(long ell, GEN mpoly, GEN T, GEN p, enum mod_type *mt, long *ptr_r)
1346 : {
1347 10003 : pari_sp ltop = avma;
1348 10003 : GEN g = gen_0;
1349 10003 : *ptr_r = 0; /*gcc -Wall*/
1350 10003 : if (!FqX_is_squarefree(mpoly, T, p)) *mt = MTcm;
1351 : else
1352 : {
1353 : long dG;
1354 9989 : g = Fq_study_eqn(mpoly, T, p, &dG, ptr_r);
1355 9989 : switch(dG)
1356 : {
1357 4711 : case 0: *mt = MTAtkin; break;
1358 539 : case 1: *mt = MTone_root; break;
1359 4669 : case 2: *mt = MTElkies; break;
1360 70 : default: *mt = (dG == ell + 1)? MTroots: MTpathological;
1361 : }
1362 : }
1363 10003 : if (DEBUGLEVEL) switch(*mt)
1364 : {
1365 0 : case MTone_root: err_printf("One root\t"); break;
1366 0 : case MTElkies: err_printf("Elkies\t"); break;
1367 0 : case MTroots: err_printf("l+1 roots\t"); break;
1368 0 : case MTAtkin: err_printf("Atkin\t"); break;
1369 0 : case MTpathological: err_printf("Pathological\n"); break;
1370 0 : case MTcm: err_printf("CM\t"); break;
1371 : }
1372 10003 : return g ? gc_GEN(ltop, g): NULL;
1373 : }
1374 :
1375 : /*Returns the trace modulo ell^k when ell is an Elkies prime */
1376 : static GEN
1377 5208 : find_trace_Elkies_power(GEN a4, GEN a6, ulong ell, long *pt_k, struct meqn *MEQN, GEN g, GEN tr, GEN q, GEN T, GEN p, long smallfact, pari_timer *ti)
1378 : {
1379 5208 : pari_sp ltop = avma, btop;
1380 : GEN tmp, Eba4, Eba6, Eca4, Eca6, Ib, kpoly;
1381 5208 : long k = *pt_k;
1382 5208 : ulong lambda, ellk = upowuu(ell, k), pellk = umodiu(q, ellk);
1383 : long cnt;
1384 :
1385 5208 : if (DEBUGLEVEL) { err_printf("mod %ld", ell); }
1386 5208 : Eba4 = a4;
1387 5208 : Eba6 = a6;
1388 5208 : tmp = find_isogenous(a4,a6, ell, MEQN, g, T, p);
1389 5208 : if (!tmp) return gc_NULL(ltop);
1390 5166 : Eca4 = gel(tmp, 1);
1391 5166 : Eca6 = gel(tmp, 2);
1392 5166 : kpoly = gel(tmp, 3);
1393 5166 : Ib = pol_x(0);
1394 5166 : lambda = tr ? find_eigen_value_oneroot(a4, a6, ell, tr, kpoly, T, p):
1395 4662 : find_eigen_value_power(a4, a6, ell, 1, 1, kpoly, T, p);
1396 5166 : if (DEBUGLEVEL>1) err_printf(" [%ld ms]", timer_delay(ti));
1397 5166 : if (smallfact && smallfact%(long)ell!=0)
1398 : {
1399 378 : ulong pell = pellk%ell;
1400 378 : ulong ap = Fl_add(lambda, Fl_div(pell, lambda, ell), ell);
1401 378 : if (Fl_sub(pell, ap, ell)==ell-1) { set_avma(ltop); return mkvecsmall(ap); }
1402 364 : if (smallfact < 0 && Fl_add(pell, ap, ell)==ell-1) { set_avma(ltop); return mkvecsmall(ap); }
1403 : }
1404 5138 : btop = avma;
1405 6433 : for (cnt = 2; cnt <= k; cnt++)
1406 : {
1407 1540 : GEN tmp = find_kernel_power(Eba4, Eba6, Eca4, Eca6, ell, MEQN, kpoly, Ib, T, p);
1408 1540 : if (!tmp) { k = cnt-1; break; }
1409 1295 : if (DEBUGLEVEL) err_printf(", %Ps", powuu(ell, cnt));
1410 1295 : lambda = find_eigen_value_power(a4, a6, ell, cnt, lambda, gel(tmp,3), T, p);
1411 1295 : Eba4 = Eca4;
1412 1295 : Eba6 = Eca6;
1413 1295 : Eca4 = gel(tmp,1);
1414 1295 : Eca6 = gel(tmp,2);
1415 1295 : kpoly = gel(tmp,4);
1416 1295 : Ib = gel(tmp, 5);
1417 1295 : if (gc_needed(btop, 1))
1418 : {
1419 0 : if(DEBUGMEM>1) pari_warn(warnmem,"find_trace_Elkies_power");
1420 0 : (void)gc_all(btop, 6, &Eba4, &Eba6, &Eca4, &Eca6, &kpoly, &Ib);
1421 : }
1422 1295 : if (DEBUGLEVEL>1) err_printf(" [%ld ms]", timer_delay(ti));
1423 : }
1424 5138 : set_avma(ltop);
1425 5138 : ellk = upowuu(ell, k);
1426 5138 : pellk = umodiu(q, ellk);
1427 5138 : *pt_k = k;
1428 5138 : return mkvecsmall(Fl_add(lambda, Fl_div(pellk, lambda, ellk), ellk));
1429 : }
1430 :
1431 : /*Returns the possible values of the trace when ell is an Atkin prime, */
1432 : /*given r the splitting degree of the modular equation at J = E.j */
1433 : static GEN
1434 4711 : find_trace_Atkin(ulong ell, long r, GEN q)
1435 : {
1436 4711 : pari_sp ltop = avma;
1437 4711 : long nval = 0;
1438 4711 : ulong teta, pell = umodiu(q, ell), invp = Fl_inv(pell, ell);
1439 4711 : GEN val_pos = cgetg(1+ell, t_VECSMALL), P = gel(factoru(r), 1);
1440 4711 : GEN S = mkvecsmall4(0, pell, 0, 1);
1441 4711 : GEN U = mkvecsmall3(0, ell-1, 0);
1442 4711 : pari_sp btop = avma;
1443 4711 : if (r==2 && krouu(ell-pell, ell) < 0)
1444 707 : val_pos[++nval] = 0;
1445 92099 : for (teta = 1; teta < ell; teta++, set_avma(btop))
1446 : {
1447 87388 : ulong disc = Fl_sub(Fl_sqr(teta,ell), Fl_mul(4UL,pell,ell), ell);
1448 : GEN a;
1449 87388 : if (krouu(disc, ell) >= 0) continue;
1450 43162 : S[3] = Fl_neg(teta, ell);
1451 43162 : U[3] = Fl_mul(invp, teta, ell);
1452 43162 : a = Flxq_powu(U, r/P[1], S, ell);
1453 43162 : if (!Flx_equal1(a) && Flx_equal1(Flxq_powu(a, P[1], S, ell)))
1454 : {
1455 29260 : pari_sp av = avma;
1456 29260 : long i, l=lg(P);
1457 49924 : for (i = 2; i < l; i++, set_avma(av))
1458 26250 : if (Flx_equal1(Flxq_powu(U, r/P[i], S, ell))) break;
1459 29260 : if (i==l) val_pos[++nval] = teta;
1460 : }
1461 : }
1462 4711 : return gc_upto(ltop, vecsmall_shorten(val_pos, nval));
1463 : }
1464 :
1465 : /*Returns the possible traces when there is only one root */
1466 : static GEN
1467 539 : find_trace_one_root(ulong ell, GEN q)
1468 : {
1469 539 : ulong a = Fl_double(Fl_sqrt(umodiu(q,ell), ell), ell);
1470 539 : return mkvecsmall2(a, ell - a);
1471 : }
1472 :
1473 : static GEN
1474 70 : find_trace_lp1_roots(long ell, GEN q)
1475 : {
1476 70 : ulong ell2 = ell * ell, pell = umodiu(q, ell2);
1477 70 : ulong a = Fl_sqrt(pell%ell, ell);
1478 70 : ulong pa = Fl_add(Fl_div(pell, a, ell2), a, ell2);
1479 70 : return mkvecsmall2(pa, ell2 - pa);
1480 : }
1481 :
1482 : /*ell odd prime; trace modulo ell^k: [], [t] or [t1,...,td] */
1483 : static GEN
1484 10003 : find_trace(GEN a4, GEN a6, GEN j, ulong ell, GEN q, GEN T, GEN p, long *ptr_kt,
1485 : long smallfact, long vx, long vy)
1486 : {
1487 10003 : pari_sp ltop = avma;
1488 : GEN g, meqnj, tr, tr2;
1489 : long kt, r;
1490 : enum mod_type mt;
1491 : struct meqn MEQN;
1492 : pari_timer ti;
1493 :
1494 10003 : kt = maxss((long)(log(expi(q)*M_LN2)/log((double)ell)), 1);
1495 10003 : if (DEBUGLEVEL)
1496 0 : { err_printf("SEA: Prime %5ld ", ell); timer_start(&ti); }
1497 10003 : get_modular_eqn(&MEQN, ell, vx, vy);
1498 10003 : meqnj = meqn_j(&MEQN, j, ell, T, p);
1499 10003 : g = study_modular_eqn(ell, meqnj, T, p, &mt, &r);
1500 : /* If l is an Elkies prime, search for a factor of the l-division polynomial.
1501 : * Then deduce the trace by looking for eigenvalues of the Frobenius by
1502 : * computing modulo this factor */
1503 10003 : switch (mt)
1504 : {
1505 539 : case MTone_root:
1506 539 : tr2 = find_trace_one_root(ell, q);
1507 539 : tr = find_trace_Elkies_power(a4,a6,ell, &kt, &MEQN, g, tr2, q, T, p, smallfact, &ti);
1508 539 : if (!tr) { tr = tr2; kt = 1; }
1509 539 : break;
1510 4669 : case MTElkies:
1511 : /* Contrary to MTone_root, may look mod higher powers of ell */
1512 4669 : if (abscmpiu(p, 2*ell+3) <= 0)
1513 49 : kt = 1; /* Not implemented in this case */
1514 4669 : tr = find_trace_Elkies_power(a4,a6,ell, &kt, &MEQN, g, NULL, q, T, p, smallfact, &ti);
1515 4669 : if (!tr)
1516 : {
1517 7 : if (DEBUGLEVEL) err_printf("[fail]\n");
1518 7 : return gc_NULL(ltop);
1519 : }
1520 4662 : break;
1521 70 : case MTroots:
1522 70 : tr = find_trace_lp1_roots(ell, q);
1523 70 : kt = 2;
1524 70 : break;
1525 4711 : case MTAtkin:
1526 4711 : tr = find_trace_Atkin(ell, r, q);
1527 4711 : if (lg(tr)==1) pari_err_PRIME("ellap",p);
1528 4711 : kt = 1;
1529 4711 : break;
1530 14 : case MTcm:
1531 : {
1532 14 : long D = find_CM(ell, j, T, p);
1533 14 : GEN C = Fq_ellcard_CM(D, a4, a6, T, p);
1534 14 : if (DEBUGLEVEL>1) err_printf(" D=%ld [%ld ms]\n", D, timer_delay(&ti));
1535 14 : return gc_const(ltop, C);
1536 : }
1537 0 : default: /* case MTpathological: */
1538 0 : return gc_NULL(ltop);
1539 : }
1540 9982 : if (DEBUGLEVEL) {
1541 0 : long n = lg(tr)-1;
1542 0 : if (n > 1 || mt == MTAtkin)
1543 : {
1544 0 : err_printf("%3ld trace(s)",n);
1545 0 : if (DEBUGLEVEL>1) err_printf(" [%ld ms]", timer_delay(&ti));
1546 : }
1547 0 : if (n > 1) err_printf("\n");
1548 : }
1549 9982 : *ptr_kt = kt;
1550 9982 : return gc_upto(ltop, tr);
1551 : }
1552 :
1553 : /* A partition of compile_atkin in baby and giant is represented as the binary
1554 : developpement of an integer; if the i-th bit is 1, the i-th prime in
1555 : compile-atkin is a baby. The optimum is obtained when the ratio between
1556 : the number of possibilities for traces modulo giants (p_g) and babies (p_b)
1557 : is near 3/4. */
1558 : static long
1559 910 : separation(GEN cnt)
1560 : {
1561 : pari_sp btop;
1562 910 : long k = lg(cnt)-1, l = (1L<<k)-1, best_i, i, j;
1563 : GEN best_r, P, P3, r;
1564 :
1565 910 : P = gen_1;
1566 4550 : for (j = 1; j <= k; ++j) P = mulis(P, cnt[j]);
1567 : /* p_b * p_g = P is constant */
1568 910 : P3 = mulsi(3, P);
1569 910 : btop = avma;
1570 910 : best_i = 0;
1571 910 : best_r = P3;
1572 44282 : for (i = 1; i < l; i++)
1573 : {
1574 : /* scan all possibilities */
1575 43463 : GEN p_b = gen_1;
1576 415947 : for (j = 0; j < k; j++)
1577 372484 : if (i & (1L<<j)) p_b = mulis(p_b, cnt[1+j]);
1578 43463 : r = subii(shifti(sqri(p_b), 2), P3); /* (p_b/p_g - 3/4)*4*P */
1579 43463 : if (!signe(r)) { best_i = i; break; }
1580 43372 : if (abscmpii(r, best_r) < 0) { best_i = i; best_r = r; }
1581 43372 : if (gc_needed(btop, 1))
1582 0 : best_r = gc_INT(btop, best_r);
1583 : }
1584 910 : return best_i;
1585 : }
1586 :
1587 : /* x VEC defined modulo P (= *P), y VECSMALL modulo q, (q,P) = 1. */
1588 : /* Update in place:
1589 : * x to vector mod q P congruent to x mod P (resp. y mod q). */
1590 : /* P ( <-- qP ) */
1591 : static void
1592 1820 : multiple_crt(GEN x, GEN y, GEN q, GEN P)
1593 : {
1594 1820 : pari_sp ltop = avma, av;
1595 1820 : long i, j, k, lx = lg(x)-1, ly = lg(y)-1;
1596 : GEN a1, a2, u, v, A2X;
1597 1820 : (void)bezout(P,q,&u,&v);
1598 1820 : a1 = mulii(P,u);
1599 1820 : a2 = mulii(q,v); A2X = ZC_Z_mul(x, a2);
1600 1820 : av = avma; affii(mulii(P,q), P);
1601 73010 : for (i = 1, k = 1; i <= lx; i++, set_avma(av))
1602 : {
1603 71190 : GEN a2x = gel(A2X,i);
1604 1194718 : for (j = 1; j <= ly; ++j)
1605 : {
1606 1123528 : GEN t = Fp_add(Fp_mulu(a1, y[j], P), a2x, P);
1607 1123528 : affii(t, gel(x, k++));
1608 : }
1609 : }
1610 1820 : setlg(x, k); set_avma(ltop);
1611 1820 : }
1612 :
1613 : /****************************************************************************/
1614 : /* MATCH AND SORT */
1615 : /****************************************************************************/
1616 :
1617 : static GEN
1618 1820 : possible_traces(GEN compile, GEN mask, GEN *P, int larger)
1619 : {
1620 1820 : GEN V, Pfinal = gen_1, C = shallowextract(compile, mask);
1621 1820 : long i, lfinal = 1, lC = lg(C), lP;
1622 1820 : pari_sp av = avma;
1623 :
1624 5460 : for (i = 1; i < lC; i++)
1625 : {
1626 3640 : GEN c = gel(C,i), t;
1627 3640 : Pfinal = mulii(Pfinal, gel(c,1));
1628 3640 : t = muluu(lfinal, lg(gel(c,2))-1);
1629 3640 : lfinal = itou(t);
1630 : }
1631 1820 : Pfinal = gc_INT(av, Pfinal);
1632 1820 : if (larger)
1633 910 : lP = lgefint(shifti(Pfinal,1));
1634 : else
1635 910 : lP = lgefint(Pfinal);
1636 1820 : lfinal++;
1637 : /* allocate room for final result */
1638 1820 : V = cgetg(lfinal, t_VEC);
1639 1061382 : for (i = 1; i < lfinal; i++) gel(V,i) = cgeti(lP);
1640 :
1641 : {
1642 1820 : GEN c = gel(C,1), v = gel(c,2);
1643 1820 : long l = lg(v);
1644 9044 : for (i = 1; i < l; i++) affsi(v[i], gel(V,i));
1645 1820 : setlg(V, l); affii(gel(c,1), Pfinal); /* reset Pfinal */
1646 : }
1647 3640 : for (i = 2; i < lC; i++)
1648 : {
1649 1820 : GEN c = gel(C,i);
1650 1820 : multiple_crt(V, gel(c,2), gel(c,1), Pfinal); /* Pfinal updated! */
1651 : }
1652 1820 : *P = Pfinal; return V;
1653 : }
1654 :
1655 : static GEN
1656 459375 : cost(long mask, GEN cost_vec)
1657 : {
1658 459375 : pari_sp ltop = avma;
1659 : long i;
1660 459375 : GEN c = gen_1;
1661 7173831 : for (i = 1; i < lg(cost_vec); i++)
1662 6714456 : if (mask&(1L<<(i-1)))
1663 2976967 : c = mulis(c, cost_vec[i]);
1664 459375 : return gc_INT(ltop, c);
1665 : }
1666 :
1667 : static GEN
1668 369894 : value(long mask, GEN atkin, long k)
1669 : {
1670 369894 : pari_sp ltop = avma;
1671 : long i;
1672 369894 : GEN c = gen_1;
1673 5777625 : for (i = 1; i <= k; i++)
1674 5407731 : if (mask&(1L<<(i-1)))
1675 2386237 : c = mulii(c, gmael(atkin, i, 1));
1676 369894 : return gc_INT(ltop, c);
1677 : }
1678 :
1679 : static void
1680 182616 : set_cost(GEN B, long b, GEN cost_vec, long *pi)
1681 : {
1682 182616 : pari_sp av = avma;
1683 182616 : GEN costb = cost(b, cost_vec);
1684 182616 : long i = *pi;
1685 250474 : while (cmpii(costb, cost(B[i], cost_vec)) < 0) --i;
1686 182616 : B[++i] = b;
1687 182616 : *pi = i; set_avma(av);
1688 182616 : }
1689 :
1690 : static GEN
1691 1925 : get_lgatkin(GEN compile_atkin, long k)
1692 : {
1693 1925 : GEN v = cgetg(k+1, t_VECSMALL);
1694 : long j;
1695 10248 : for (j = 1; j <= k; ++j) v[j] = lg(gmael(compile_atkin, j, 2))-1;
1696 1925 : return v;
1697 : }
1698 :
1699 : static GEN
1700 1015 : champion(GEN atkin, long k, GEN bound_champ)
1701 : {
1702 1015 : const long two_k = 1L<<k;
1703 1015 : pari_sp ltop = avma;
1704 : long i, j, n, i1, i2;
1705 1015 : GEN B, Bp, cost_vec, res = NULL;
1706 :
1707 1015 : cost_vec = get_lgatkin(atkin, k);
1708 1015 : if (k == 1) return mkvec2(gen_1, utoipos(cost_vec[1]));
1709 :
1710 1001 : B = zero_zv(two_k);
1711 1001 : Bp = zero_zv(two_k);
1712 1001 : Bp[2] = 1;
1713 4669 : for (n = 2, j = 2; j <= k; j++)
1714 : {
1715 : long b;
1716 3668 : i = 1;
1717 173418 : for (i1 = 2, i2 = 1; i1 <= n; )
1718 : {
1719 169750 : pari_sp av = avma;
1720 169750 : long b1 = Bp[i1], b2 = Bp[i2]|(1L<<(j-1));
1721 169750 : if (cmpii(value(b1, atkin, k), value(b2, atkin, k)) < 0)
1722 169750 : { b = b1; i1++; } else { b = b2; i2++; }
1723 169750 : set_avma(av);
1724 169750 : set_cost(B, b, cost_vec, &i);
1725 : }
1726 16534 : for ( ; i2 <= n; i2++)
1727 : {
1728 12866 : b = Bp[i2]|(1L<<(j-1));
1729 12866 : set_cost(B, b, cost_vec, &i);
1730 : }
1731 3668 : n = i;
1732 122094 : for (i = 1; i <= n; i++)
1733 118426 : Bp[i] = B[i];
1734 : }
1735 9631069 : for (i = 1; i <= two_k; i++)
1736 9630068 : if (B[i])
1737 : {
1738 26285 : GEN b = cost (B[i], cost_vec);
1739 26285 : GEN v = value(B[i], atkin, k);
1740 26285 : if (cmpii(v, bound_champ) <=0) continue;
1741 5005 : if (res && gcmp(b, gel(res, 2)) >=0) continue;
1742 1001 : res = mkvec2(utoi(B[i]), b);
1743 : }
1744 1001 : return gc_GEN(ltop, res);
1745 : }
1746 :
1747 : static GEN
1748 1820 : compute_diff(GEN v)
1749 : {
1750 1820 : long i, l = lg(v) - 1;
1751 1820 : GEN diff = cgetg(l, t_VEC);
1752 1059562 : for (i = 1; i < l; i++) gel(diff, i) = subii(gel(v, i+1), gel(v, i));
1753 1820 : return ZV_sort_uniq_shallow(diff);
1754 : }
1755 :
1756 : static int
1757 17276 : cmp_atkin(void*E, GEN a, GEN b)
1758 : {
1759 17276 : long ta=typ(a)==t_INT, tb=typ(b)==t_INT, c;
1760 : (void) E;
1761 17276 : if (ta || tb) return ta-tb;
1762 5670 : c = lg(gel(a,2)) - lg(gel(b,2));
1763 5670 : if (c) return c;
1764 847 : return cmpii(gel(b,1), gel(a,1));
1765 : }
1766 :
1767 : static void
1768 4109 : add_atkin(GEN atkin, GEN trace, long *nb)
1769 : {
1770 4109 : long l = lg(atkin)-1;
1771 4109 : long i, k = gen_search(atkin, trace, NULL, cmp_atkin);
1772 4109 : if (k > 0 || (k = -k) > l) return;
1773 79926 : for (i = l; i > k; i--) gel(atkin,i) = gel(atkin,i-1);
1774 4109 : if (typ(gel(atkin,l))==t_INT) (*nb)++;
1775 4109 : gel(atkin,k) = trace;
1776 : }
1777 :
1778 : /* V = baby / giant, P = Pb / Pg */
1779 : static GEN
1780 1820 : BSGS_pre(GEN *pdiff, GEN V, GEN P, void *E, const struct bb_group *grp)
1781 : {
1782 1820 : GEN diff = compute_diff(V);
1783 1820 : GEN pre = cgetg(lg(diff), t_VEC);
1784 1820 : long i, l = lg(diff);
1785 1820 : gel(pre, 1) = grp->pow(E, P, gel(diff, 1));
1786 : /* what we'd _really_ want here is a hashtable diff[i] -> pre[i] */
1787 39018 : for (i = 2; i < l; i++)
1788 : {
1789 37198 : pari_sp av = avma;
1790 37198 : GEN d = subii(gel(diff, i), gel(diff, i-1));
1791 37198 : GEN Q = grp->mul(E, gel(pre, i-1), grp->pow(E, P, d));
1792 37198 : gel(pre, i) = gc_GEN(av, Q);
1793 : }
1794 1820 : *pdiff = diff; return pre;
1795 : }
1796 :
1797 : /* u = trace_elkies, Mu = prod_elkies. Let caller collect garbage */
1798 : /* Match & sort: variant from Lercier's thesis, section 11.2.3 */
1799 : /* baby/giant/table updated in place: this routines uses
1800 : * size(baby)+size(giant)+size(table)+size(table_ind) + O(log p)
1801 : * bits of stack */
1802 : static GEN
1803 966 : match_and_sort(GEN compile_atkin, GEN Mu, GEN u, GEN q, void *E, const struct bb_group *grp)
1804 : {
1805 : pari_sp av1, av2;
1806 966 : GEN baby, giant, SgMb, Mb, Mg, den, Sg, dec_inf, div, pp1 = addiu(q,1);
1807 : GEN P, Pb, Pg, point, diff, pre, table, table_ind;
1808 966 : long best_i, i, lbaby, lgiant, k = lg(compile_atkin)-1;
1809 966 : GEN bound = sqrti(shifti(q, 2)), card;
1810 966 : const long lcard = 100;
1811 966 : long lq = lgefint(q), nbcard;
1812 : pari_timer ti;
1813 :
1814 966 : if (k == 1)
1815 : { /*only one Atkin prime, check the cardinality with random points */
1816 56 : GEN r = gel(compile_atkin, 1), r1 = gel(r,1), r2 = gel(r,2);
1817 56 : long l = lg(r2), j;
1818 56 : GEN card = cgetg(l, t_VEC), Cs2, C, U;
1819 56 : Z_chinese_pre(Mu, r1, &C,&U, NULL);
1820 56 : Cs2 = shifti(C, -1);
1821 378 : for (j = 1, i = 1; i < l; i++)
1822 : {
1823 322 : GEN t = Z_chinese_post(u, stoi(r2[i]), C, U, NULL);
1824 322 : t = Fp_center_i(t, C, Cs2);
1825 322 : if (abscmpii(t, bound) <= 0) gel(card, j++) = subii(pp1, t);
1826 : }
1827 56 : setlg(card, j);
1828 56 : return gen_select_order(card, E, grp);
1829 : }
1830 910 : if (DEBUGLEVEL>=2) timer_start(&ti);
1831 910 : av1 = avma;
1832 910 : best_i = separation( get_lgatkin(compile_atkin, k) );
1833 910 : set_avma(av1);
1834 :
1835 910 : baby = possible_traces(compile_atkin, utoi(best_i), &Mb, 1);
1836 910 : giant = possible_traces(compile_atkin, subiu(int2n(k), best_i+1), &Mg, 0);
1837 910 : lbaby = lg(baby);
1838 910 : lgiant = lg(giant);
1839 910 : den = Fp_inv(Fp_mul(Mu, Mb, Mg), Mg);
1840 910 : av2 = avma;
1841 622790 : for (i = 1; i < lgiant; i++, set_avma(av2))
1842 621880 : affii(Fp_mul(gel(giant,i), den, Mg), gel(giant,i));
1843 910 : ZV_sort_inplace(giant);
1844 910 : Sg = Fp_mul(negi(u), den, Mg);
1845 910 : den = Fp_inv(Fp_mul(Mu, Mg, Mb), Mb);
1846 910 : dec_inf = divii(mulii(Mb,addii(Mg,shifti(Sg,1))), shifti(Mg,1));
1847 910 : togglesign(dec_inf); /* now, dec_inf = ceil(- (Mb/2 + Sg Mb/Mg) ) */
1848 910 : div = mulii(truedivii(dec_inf, Mb), Mb);
1849 910 : av2 = avma;
1850 438592 : for (i = 1; i < lbaby; i++, set_avma(av2))
1851 : {
1852 437682 : GEN b = addii(Fp_mul(Fp_sub(gel(baby,i), u, Mb), den, Mb), div);
1853 437682 : if (cmpii(b, dec_inf) < 0) b = addii(b, Mb);
1854 437682 : affii(b, gel(baby,i));
1855 : }
1856 910 : ZV_sort_inplace(baby);
1857 :
1858 910 : SgMb = mulii(Sg, Mb);
1859 910 : card = cgetg(lcard+1,t_VEC);
1860 91910 : for (i = 1; i <= lcard; i++) gel(card,i) = cgetipos(lq+1);
1861 :
1862 910 : av2 = avma;
1863 910 : MATCH_RESTART:
1864 910 : set_avma(av2);
1865 910 : nbcard = 0;
1866 910 : P = grp->rand(E);
1867 910 : point = grp->pow(E,P, Mu);
1868 910 : Pb = grp->pow(E,point, Mg);
1869 910 : Pg = grp->pow(E,point, Mb);
1870 : /* Precomputation for babies */
1871 910 : pre = BSGS_pre(&diff, baby, Pb, E, grp);
1872 :
1873 : /*Now we compute the table of babies, this table contains only the */
1874 : /*lifted x-coordinate of the points in order to use less memory */
1875 910 : table = cgetg(lbaby, t_VECSMALL);
1876 910 : av1 = avma;
1877 : /* (p+1 - u - Mu*Mb*Sg) P - (baby[1]) Pb */
1878 910 : point = grp->pow(E,P, subii(subii(pp1, u), mulii(Mu, addii(SgMb, mulii(Mg, gel(baby,1))))));
1879 910 : table[1] = grp->hash(gel(point,1));
1880 437682 : for (i = 2; i < lbaby; i++)
1881 : {
1882 436772 : GEN d = subii(gel(baby, i), gel(baby, i-1));
1883 436772 : point = grp->mul(E, point, grp->pow(E, gel(pre, ZV_search(diff, d)), gen_m1));
1884 436772 : table[i] = grp->hash(gel(point,1));
1885 436772 : if (gc_needed(av1,3))
1886 : {
1887 19 : if(DEBUGMEM>1) pari_warn(warnmem,"match_and_sort, baby = %ld", i);
1888 19 : point = gc_upto(av1, point);
1889 : }
1890 : }
1891 910 : set_avma(av1);
1892 : /* Precomputations for giants */
1893 910 : pre = BSGS_pre(&diff, giant, Pg, E, grp);
1894 :
1895 : /* Look for a collision among the x-coordinates */
1896 910 : table_ind = vecsmall_indexsort(table);
1897 910 : table = perm_mul(table,table_ind);
1898 :
1899 910 : av1 = avma;
1900 910 : point = grp->pow(E, Pg, gel(giant, 1));
1901 910 : for (i = 1; ; i++)
1902 620970 : {
1903 : GEN d;
1904 621880 : long h = grp->hash(gel(point, 1));
1905 621880 : long s = zv_search(table, h);
1906 621880 : if (s) {
1907 1820 : while (table[s] == h && s) s--;
1908 1820 : for (s++; s < lbaby && table[s] == h; s++)
1909 : {
1910 910 : GEN B = gel(baby,table_ind[s]), G = gel(giant,i);
1911 910 : GEN GMb = mulii(G, Mb), BMg = mulii(B, Mg);
1912 910 : GEN Be = subii(subii(pp1, u), mulii(Mu, addii(SgMb, BMg)));
1913 910 : GEN Bp = grp->pow(E,P, Be);
1914 : /* p+1 - u - Mu (Sg Mb + GIANT Mb + BABY Mg) */
1915 910 : if (gequal(gel(Bp,1),gel(point,1)))
1916 : {
1917 910 : GEN card1 = subii(Be, mulii(Mu, GMb));
1918 910 : GEN card2 = addii(card1, mulii(mulsi(2,Mu), GMb));
1919 910 : if (abscmpii(subii(pp1, card1), bound) <= 0)
1920 798 : affii(card1, gel(card, ++nbcard));
1921 910 : if (nbcard >= lcard) goto MATCH_RESTART;
1922 910 : if (abscmpii(subii(pp1, card2), bound) <= 0)
1923 490 : affii(card2, gel(card, ++nbcard));
1924 910 : if (nbcard >= lcard) goto MATCH_RESTART;
1925 : }
1926 : }
1927 : }
1928 621880 : if (i==lgiant-1) break;
1929 620970 : d = subii(gel(giant, i+1), gel(giant, i));
1930 620970 : point = grp->mul(E,point, gel(pre, ZV_search(diff, d)));
1931 620970 : if (gc_needed(av1,3))
1932 : {
1933 26 : if(DEBUGMEM>1) pari_warn(warnmem,"match_and_sort, giant = %ld", i);
1934 26 : point = gc_upto(av1, point);
1935 : }
1936 : }
1937 910 : setlg(card, nbcard+1);
1938 910 : if (DEBUGLEVEL>=2) timer_printf(&ti,"match_and_sort");
1939 910 : return gen_select_order(card, E, grp);
1940 : }
1941 :
1942 : static GEN
1943 1015 : get_bound_bsgs(long lp)
1944 : {
1945 : GEN B;
1946 1015 : if (lp <= 160)
1947 980 : B = divru(powru(dbltor(1.048), lp), 9);
1948 35 : else if (lp <= 192)
1949 28 : B = divrr(powru(dbltor(1.052), lp), dbltor(16.65));
1950 : else
1951 7 : B = mulrr(powru(dbltor(1.035), minss(lp,307)), dbltor(1.35));
1952 1015 : return mulru(B, 1000000);
1953 : }
1954 :
1955 : /* E is an elliptic curve defined over Z or over Fp in ellinit format, defined
1956 : * by the equation E: y^2 + a1*x*y + a2*y = x^3 + a2*x^2 + a4*x + a6
1957 : * p is a prime number
1958 : * set smallfact to stop whenever a small factor of the order, not dividing smallfact,
1959 : * is detected. Useful when searching for a good curve for cryptographic
1960 : * applications */
1961 : GEN
1962 1064 : Fq_ellcard_SEA(GEN a4, GEN a6, GEN q, GEN T, GEN p, long smallfact)
1963 : {
1964 1064 : const long MAX_ATKIN = 21;
1965 1064 : pari_sp ltop = avma, btop;
1966 : long ell, i, nb_atkin, vx,vy;
1967 : GEN TR, TR_mod, compile_atkin, bound, bound_bsgs, champ;
1968 1064 : GEN prod_atkin = gen_1, max_traces = gen_0;
1969 : GEN j;
1970 1064 : double bound_gr = 1.;
1971 1064 : const double growth_factor = 1.26;
1972 : forprime_t TT;
1973 :
1974 1064 : j = Fq_ellj(a4, a6, T, p);
1975 1064 : if (signe(j) == 0 || signe(Fq_sub(j, utoi(1728), T, p)) == 0)
1976 0 : return T ? FpXQ_ellcard(Fq_to_FpXQ(a4, T, p), Fq_to_FpXQ(a6, T, p), T, p)
1977 14 : : Fp_ellcard(a4, a6, p);
1978 1050 : if (Fq_elljissupersingular(j, T, p))
1979 21 : return Fq_ellcard_supersingular(a4, a6, T, p);
1980 : /*First compute the trace modulo 2 */
1981 1029 : switch(FqX_nbroots(rhs(a4, a6, 0), T, p))
1982 : {
1983 70 : case 3: /* bonus time: 4 | #E(Fq) = q+1 - t */
1984 70 : i = mod4(q)+1; if (i > 2) i -= 4;
1985 70 : TR_mod = utoipos(4);
1986 70 : TR = stoi(i); break;
1987 511 : case 1:
1988 511 : TR_mod = gen_2;
1989 511 : TR = gen_0; break;
1990 448 : default : /* 0 */
1991 448 : TR_mod = gen_2;
1992 448 : TR = gen_1; break;
1993 : }
1994 1029 : if (odd(smallfact) && !mpodd(TR))
1995 : {
1996 14 : if (DEBUGLEVEL) err_printf("Aborting: #E(Fq) divisible by 2\n");
1997 14 : set_avma(ltop); return gen_0;
1998 : }
1999 1015 : vy = fetch_var();
2000 1015 : vx = fetch_var_higher();
2001 :
2002 : /* compile_atkin is a vector containing informations about Atkin primes,
2003 : * informations about Elkies primes lie in Mod(TR, TR_mod). */
2004 1015 : u_forprime_init(&TT, 3, ULONG_MAX);
2005 1015 : bound = sqrti(shifti(q, 4));
2006 1015 : bound_bsgs = get_bound_bsgs(expi(q));
2007 1015 : compile_atkin = zerovec(MAX_ATKIN); nb_atkin = 0;
2008 1015 : btop = avma;
2009 10010 : while ( (ell = u_forprime_next(&TT)) )
2010 : {
2011 10010 : long ellkt, kt = 1, nbtrace;
2012 : GEN trace_mod;
2013 10017 : if (absequalui(ell, p)) continue;
2014 10003 : trace_mod = find_trace(a4, a6, j, ell, q, T, p, &kt, smallfact, vx,vy);
2015 10003 : if (!trace_mod) continue;
2016 9996 : if (typ(trace_mod)==t_INT)
2017 : {
2018 14 : delete_var(); delete_var();
2019 1015 : return gc_INT(ltop, trace_mod);
2020 : }
2021 9982 : nbtrace = lg(trace_mod) - 1;
2022 9982 : ellkt = (long)upowuu(ell, kt);
2023 9982 : if (nbtrace == 1)
2024 : {
2025 5873 : long t_mod_ellkt = trace_mod[1];
2026 5873 : if (smallfact && smallfact%ell!=0)
2027 : { /* does ell divide q + 1 - t ? */
2028 385 : long q_mod_ell_plus_one = umodiu(q,ell) + 1;
2029 385 : ulong card_mod_ell = umodsu(q_mod_ell_plus_one - t_mod_ellkt, ell);
2030 385 : ulong tcard_mod_ell = 1;
2031 385 : if (card_mod_ell && smallfact < 0)
2032 133 : tcard_mod_ell = umodsu(q_mod_ell_plus_one + t_mod_ellkt, ell);
2033 385 : if (!card_mod_ell || !tcard_mod_ell)
2034 : {
2035 28 : if (DEBUGLEVEL)
2036 0 : err_printf("\nAborting: #E%s(Fq) divisible by %ld\n",
2037 : tcard_mod_ell ? "" : "_twist", ell);
2038 28 : delete_var(); delete_var();
2039 28 : return gc_const(ltop, gen_0);
2040 : }
2041 : }
2042 5845 : (void)Z_incremental_CRT(&TR, t_mod_ellkt, &TR_mod, ellkt);
2043 5845 : if (DEBUGLEVEL)
2044 0 : err_printf(", missing %ld bits\n",expi(bound)-expi(TR_mod));
2045 : }
2046 : else
2047 : {
2048 4109 : add_atkin(compile_atkin, mkvec2(utoipos(ellkt), trace_mod), &nb_atkin);
2049 4109 : prod_atkin = value(-1, compile_atkin, nb_atkin);
2050 : }
2051 9954 : if (cmpii(mulii(TR_mod, prod_atkin), bound) > 0)
2052 : {
2053 : GEN bound_tr;
2054 1057 : if (!nb_atkin)
2055 : {
2056 7 : delete_var(); delete_var();
2057 7 : return gc_INT(ltop, subii(addiu(q, 1), TR));
2058 : }
2059 1050 : bound_tr = mulrr(bound_bsgs, dbltor(bound_gr));
2060 1050 : bound_gr *= growth_factor;
2061 1050 : if (signe(max_traces))
2062 : {
2063 84 : max_traces = divis(muliu(max_traces,nbtrace), ellkt);
2064 84 : if (DEBUGLEVEL>=3)
2065 0 : err_printf("At least %Ps remaining possibilities.\n",max_traces);
2066 : }
2067 1050 : if (cmpir(max_traces, bound_tr) < 0)
2068 : {
2069 1015 : GEN bound_atkin = truedivii(bound, TR_mod);
2070 1015 : champ = champion(compile_atkin, nb_atkin, bound_atkin);
2071 1015 : max_traces = gel(champ,2);
2072 1015 : if (DEBUGLEVEL>=2)
2073 0 : err_printf("%Ps remaining possibilities.\n", max_traces);
2074 1015 : if (cmpir(max_traces, bound_tr) < 0)
2075 : {
2076 966 : GEN res, cat = shallowextract(compile_atkin, gel(champ,1));
2077 : const struct bb_group *grp;
2078 : void *E;
2079 966 : if (DEBUGLEVEL)
2080 0 : err_printf("Match and sort for %Ps possibilities.\n", max_traces);
2081 966 : delete_var();
2082 966 : delete_var();
2083 966 : grp = get_FqE_group(&E,a4,a6,T,p);
2084 966 : res = match_and_sort(cat, TR_mod, TR, q, E, grp);
2085 966 : return gc_INT(ltop, res);
2086 : }
2087 : }
2088 : }
2089 8981 : if (gc_needed(btop, 1))
2090 0 : (void)gc_all(btop,5, &TR,&TR_mod, &compile_atkin, &max_traces, &prod_atkin);
2091 : }
2092 : return NULL;/*LCOV_EXCL_LINE*/
2093 : }
2094 :
2095 : GEN
2096 973 : Fp_ellcard_SEA(GEN a4, GEN a6, GEN p, long smallfact)
2097 973 : { return Fq_ellcard_SEA(a4, a6, p, NULL, p, smallfact); }
|
