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