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