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 1628 : pari_set_seadata(GEN mod) { modular_eqn = mod; }
43 : GEN
44 0 : pari_get_seadata(void) { return modular_eqn; }
45 :
46 : static char *
47 78 : seadata_filename(ulong ell)
48 78 : { return stack_sprintf("%s/seadata/sea%ld", pari_datadir, ell); }
49 :
50 : static GEN
51 78 : get_seadata(ulong ell)
52 : {
53 78 : pari_sp av = avma;
54 : GEN eqn;
55 78 : char *s = seadata_filename(ell);
56 78 : pariFILE *F = pari_fopengz(s);
57 78 : if (!F) return NULL;
58 30 : if (ell) /* large single polynomial */
59 6 : eqn = gp_read_stream(F->file);
60 : else
61 : { /* table of polynomials of small level */
62 24 : eqn = gp_readvec_stream(F->file);
63 24 : modular_eqn = eqn = gclone(eqn);
64 24 : set_avma(av);
65 : }
66 30 : pari_fclose(F);
67 30 : return eqn;
68 : }
69 :
70 : /*Builds the modular equation corresponding to the vector list. Shallow */
71 : static GEN
72 8544 : list_to_pol(GEN list, long vx, long vy)
73 : {
74 8544 : long i, l = lg(list);
75 8544 : GEN P = cgetg(l, t_VEC);
76 175818 : for (i = 1; i < l; i++)
77 : {
78 167274 : GEN L = gel(list,i);
79 167274 : if (typ(L) == t_VEC) L = RgV_to_RgX_reverse(L, vy);
80 167274 : gel(P, i) = L;
81 : }
82 8544 : 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 8592 : seadata_cache(ulong ell)
93 : {
94 8592 : long n = uprimepi(ell)-1;
95 : GEN C;
96 8592 : if (!modular_eqn && !get_seadata(0))
97 48 : C = NULL;
98 8544 : else if (n && n < lg(modular_eqn))
99 8538 : C = gel(modular_eqn, n);
100 : else
101 6 : C = get_seadata(ell);
102 8592 : return C;
103 : }
104 : /* C = [prime level, type "A" or "C", pol. coeffs] */
105 : static void
106 8544 : seadata_parse(struct meqn *M, GEN C, long vx, long vy)
107 : {
108 8544 : M->type = *GSTR(gel(C,2));
109 8544 : M->eq = list_to_pol(gel(C,3), vx, vy);
110 8544 : }
111 : static void
112 8574 : get_modular_eqn(struct meqn *M, ulong ell, long vx, long vy)
113 : {
114 8574 : GEN C = seadata_cache(ell);
115 8574 : M->vx = vx;
116 8574 : M->vy = vy;
117 8574 : M->eval = gen_0;
118 8574 : if (C) seadata_parse(M, C, vx, vy);
119 : else
120 : {
121 48 : M->type = 'J'; /* j^(1/3) for ell != 3, j for 3 */
122 48 : M->eq = polmodular_ZXX(ell, ell==3? 0: 5, vx, vy);
123 : }
124 8574 : }
125 :
126 : GEN
127 30 : ellmodulareqn(long ell, long vx, long vy)
128 : {
129 30 : pari_sp av = avma;
130 : struct meqn meqn;
131 : GEN C;
132 30 : if (vx < 0) vx = 0;
133 30 : if (vy < 0) vy = 1;
134 30 : if (varncmp(vx,vy) >= 0)
135 6 : pari_err_PRIORITY("ellmodulareqn", pol_x(vx), ">=", vy);
136 24 : if (ell < 2 || !uisprime(ell))
137 6 : pari_err_PRIME("ellmodulareqn (level)", stoi(ell));
138 18 : C = seadata_cache(ell);
139 18 : if (!C) pari_err_FILE("seadata file", seadata_filename(ell));
140 18 : seadata_parse(&meqn, C, vx, vy);
141 18 : 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 104 : Fq_to_Flx(GEN a4, GEN T, ulong p)
152 104 : { 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 840 : get_FqE_group(void ** pt_E, GEN a4, GEN a6, GEN T, GEN p)
157 : {
158 840 : if (!T) return get_FpE_group(pt_E,a4,a6,p);
159 66 : else if (lgefint(p)==3)
160 : {
161 52 : ulong pp = uel(p,2);
162 52 : GEN Tp = ZXT_to_FlxT(T,pp);
163 52 : return get_FlxqE_group(pt_E, Fq_to_Flx(a4, Tp, pp), Fq_to_Flx(a6, Tp, pp),
164 : Tp, pp);
165 : }
166 14 : 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 12 : list_singular_discs(long l)
179 : {
180 12 : const long _4l2 = 4*l*l;
181 : long v;
182 12 : GEN V = zero_F2v(_4l2);
183 : /* special cased for efficiency + remove factor l^2 from conductor */
184 12 : F2v_set(V, 4); /* v = 0 */
185 12 : F2v_set(V, 3); /* v = l */
186 1056 : for (v = 1; v < 2*l; v++)
187 1044 : if (v != l)
188 : { /* l does not divide _4l2 - v*v */
189 1032 : GEN F = factoru(_4l2 - v*v), P, E, c;
190 1032 : ulong d = coredisc2u_fact(F, -1, &P, &E);
191 : long i, lc;
192 1032 : c = divisorsu_fact(mkvec2(P,E));
193 1032 : lc = lg(c);
194 3024 : for (i = 1; i < lc; i++)
195 1992 : F2v_set(V, d * uel(c,i)*uel(c,i));
196 : }
197 12 : 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 12 : find_CM(long l, GEN j, GEN T, GEN p)
204 : {
205 12 : const long inv = 0;
206 12 : GEN v = list_singular_discs(l);
207 12 : long i, n = v[1];
208 12 : GEN db = polmodular_db_init(inv);
209 738 : for (i = 1; i < n; i++)
210 738 : if (F2v_coeff(v,i))
211 : {
212 138 : GEN C = polclass0(-i, inv, 0, &db);
213 138 : GEN F = FqX_eval(C, j, T, p);
214 138 : if (signe(F)==0) break;
215 : }
216 12 : gunclone_deep(db); return i < n ? -i: 0;
217 : }
218 :
219 : static GEN
220 12 : vecpoints_to_vecx(GEN x, GEN q1)
221 : {
222 36 : 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 12 : Fq_ellcard_CM(long disc, GEN a4, GEN a6, GEN T, GEN p)
228 : {
229 : const struct bb_group *grp;
230 : void *E;
231 12 : long d = T ? degpol(T): 1;
232 12 : GEN q = powiu(p, d), q1 = addiu(q, 1), Q, S;
233 12 : Q = qfbsolve(Qfb0(gen_1,gen_0,stoi(-disc)), mkmat22(gen_2, gen_2, p, utoi(d)), 3);
234 12 : if (lg(Q)==1) return q1;
235 12 : S = vecpoints_to_vecx(Q, q1);
236 12 : grp = get_FqE_group(&E, a4, a6, T, p);
237 12 : 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 124464 : divpol_f2(GEN t, GEN r2, long n, void *E, const struct bb_algebra *ff)
251 : {
252 124464 : if (n==0) return ff->zero(E);
253 124464 : if (n<=2) return ff->one(E);
254 103224 : if (gmael(t,2,n)) return gmael(t,2,n);
255 38064 : gmael(t,2,n) = gclone(ff->sqr(E,divpol(t,r2,n,E,ff)));
256 38064 : return gmael(t,2,n);
257 : }
258 :
259 : /* f_n f_{n-2}, return ff->zero or a clone */
260 : static GEN
261 75612 : divpol_ff(GEN t, GEN r2, long n, void *E, const struct bb_algebra *ff)
262 : {
263 75612 : if (n<=2) return ff->zero(E);
264 75612 : if (gmael(t,3,n)) return gmael(t,3,n);
265 48672 : if (n<=4) return divpol(t,r2,n,E,ff);
266 21438 : gmael(t,3,n) = gclone(ff->mul(E,divpol(t,r2,n,E,ff), divpol(t,r2,n-2,E,ff)));
267 21438 : return gmael(t,3,n);
268 : }
269 :
270 : /* f_n, return ff->zero or a clone */
271 : static GEN
272 161184 : divpol(GEN t, GEN r2, long n, void *E, const struct bb_algebra *ff)
273 : {
274 161184 : long m = n/2;
275 161184 : pari_sp av = avma;
276 : GEN f;
277 161184 : if (n==0) return ff->zero(E);
278 157920 : if (gmael(t,1,n)) return gmael(t,1,n);
279 43998 : switch(n)
280 : {
281 6192 : case 1:
282 : case 2:
283 6192 : f = ff->one(E);
284 6192 : break;
285 37806 : default:
286 37806 : if (odd(n))
287 21966 : if (odd(m))
288 9756 : f = ff->sub(E, ff->mul(E, divpol_ff(t,r2,m+2,E,ff),
289 : divpol_f2(t,r2,m,E,ff)),
290 9756 : ff->mul(E, r2,
291 9756 : ff->mul(E,divpol_ff(t,r2,m+1,E,ff),
292 : divpol_f2(t,r2,m+1,E,ff))));
293 : else
294 24420 : f = ff->sub(E, ff->mul(E, r2,
295 12210 : ff->mul(E, divpol_ff(t,r2,m+2,E,ff),
296 : divpol_f2(t,r2,m,E,ff))),
297 12210 : ff->mul(E, divpol_ff(t,r2,m+1,E,ff),
298 : divpol_f2(t,r2,m+1,E,ff)));
299 : else
300 15840 : 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 15840 : ff->mul(E, divpol_ff(t,r2,m,E,ff),
303 : divpol_f2(t,r2,m+1,E,ff)));
304 : }
305 43998 : gmael(t,1,n) = f = gclone( ff->red(E, f) );
306 43998 : set_avma(av); return f;
307 : }
308 :
309 : static void
310 9696 : divpol_free(GEN t)
311 : {
312 9696 : long i, l = lg(gel(t,1));
313 179730 : for (i=1; i<l; i++)
314 : {
315 170034 : guncloneNULL(gmael(t,1,i));
316 170034 : guncloneNULL(gmael(t,2,i));
317 170034 : guncloneNULL(gmael(t,3,i));
318 : }
319 9696 : }
320 :
321 : static GEN
322 1296 : Flxq_elldivpol34(long n, GEN a4, GEN a6, GEN S, GEN T, ulong p)
323 : {
324 : GEN res;
325 1296 : long vs = T[1];
326 1296 : switch(n)
327 : {
328 648 : case 3:
329 648 : 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 648 : break;
332 648 : case 4:
333 : {
334 648 : GEN a42 = Flxq_sqr(a4, T, p);
335 1296 : 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 648 : Flx_sub(Flx_mulu(Flxq_sqr(a6, T, p), p-8%p, p),
339 : Flxq_mul(a4, a42, T, p), p));
340 648 : res = FlxX_double(res, p);
341 : }
342 648 : break;
343 0 : default:
344 0 : pari_err_BUG("Flxq_elldivpol34");
345 : return NULL;/*LCOV_EXCL_LINE*/
346 : }
347 1296 : setvarn(res, get_FlxqX_var(S));
348 1296 : return FlxqX_rem(res, S, T, p);
349 : }
350 :
351 : static GEN
352 18096 : Fq_elldivpol34(long n, GEN a4, GEN a6, GEN S, GEN T, GEN p)
353 : {
354 : GEN res;
355 18096 : switch(n)
356 : {
357 9048 : case 3:
358 9048 : 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 9048 : break;
361 9048 : case 4:
362 : {
363 9048 : GEN a42 = Fq_sqr(a4, T, p);
364 9048 : 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 9048 : res = FqX_mulu(res, 2, T, p);
370 : }
371 9048 : break;
372 0 : default:
373 0 : pari_err_BUG("Fq_elldivpol34");
374 : return NULL;/*LCOV_EXCL_LINE*/
375 : }
376 18096 : if (S)
377 : {
378 18096 : setvarn(res, get_FpXQX_var(S));
379 18096 : res = FqX_rem(res, S, T, p);
380 : }
381 18096 : return res;
382 : }
383 :
384 : static GEN
385 15152 : rhs(GEN a4, GEN a6, long v)
386 : {
387 15152 : GEN RHS = mkpoln(4, gen_1, gen_0, a4, a6);
388 15152 : setvarn(RHS, v); return RHS;
389 : }
390 :
391 : static GEN
392 964 : Flxq_rhs(GEN a4, GEN a6, long v, long vs)
393 : {
394 964 : GEN RHS = mkpoln(4, pol1_Flx(vs), pol0_Flx(vs), a4, a6);
395 964 : 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 9696 : divpolmod_init(struct divpolmod_red *d, GEN D3, GEN D4, GEN RHS, long n,
407 : void *E, const struct bb_algebra *ff)
408 : {
409 9696 : long k = n+2;
410 9696 : d->ff = ff; d->E = E;
411 9696 : d->t = mkvec3(const_vec(k, NULL),const_vec(k, NULL),const_vec(k, NULL));
412 9696 : if (k>=3) gmael(d->t,1,3) = gclone(D3);
413 9696 : if (k>=4) gmael(d->t,1,4) = gclone(D4);
414 9696 : d->r2 = ff->sqr(E, RHS);
415 9696 : }
416 :
417 : static void
418 9048 : 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 9048 : GEN RHS, D3 = NULL, D4 = NULL;
423 9048 : long v = h ? get_FpXQX_var(h): 0;
424 9048 : D3 = n>=0 ? Fq_elldivpol34(3, a4, a6, h, T, p): NULL;
425 9048 : D4 = n>=1 ? Fq_elldivpol34(4, a4, a6, h, T, p): NULL;
426 9048 : RHS = rhs(a4, a6, v);
427 9048 : RHS = h ? FqX_rem(RHS, h, T, p): RHS;
428 9048 : RHS = FqX_mulu(RHS, 4, T, p);
429 9048 : 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 9048 : divpolmod_init(d, D3, D4, RHS, n, E, ff);
432 9048 : }
433 :
434 : static void
435 648 : 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 648 : GEN RHS, D3 = NULL, D4 = NULL;
440 648 : long v = get_FlxqX_var(h), vT = get_Flx_var(T);
441 648 : D3 = n>=0 ? Flxq_elldivpol34(3, a4, a6, h, T, p): NULL;
442 648 : D4 = n>=1 ? Flxq_elldivpol34(4, a4, a6, h, T, p): NULL;
443 648 : RHS = FlxX_Fl_mul(FlxqX_rem(Flxq_rhs(a4, a6, v, vT), h, T, p), 4, p);
444 648 : ff = get_FlxqXQ_algebra(&E, h, T, p);
445 648 : divpolmod_init(d, D3, D4, RHS, n, E, ff);
446 648 : }
447 :
448 : /*Computes the n-division polynomial modulo the polynomial h \in Fq[x] */
449 : GEN
450 332 : Flxq_elldivpolmod(GEN a4, GEN a6, long n, GEN h, GEN T, ulong p)
451 : {
452 : struct divpolmod_red d;
453 332 : pari_sp ltop = avma;
454 : GEN res;
455 332 : Flxq_elldivpolmod_init(&d, a4, a6, n, h, T, p);
456 332 : res = gcopy(divpol(d.t,d.r2,n,d.E,d.ff));
457 332 : divpol_free(d.t);
458 332 : return gerepileupto(ltop, res);
459 : }
460 :
461 : /*Computes the n-division polynomial modulo the polynomial h \in Fq[x] */
462 : GEN
463 4158 : Fq_elldivpolmod(GEN a4, GEN a6, long n, GEN h, GEN T, GEN p)
464 : {
465 : struct divpolmod_red d;
466 4158 : pari_sp ltop = avma;
467 : GEN res;
468 4158 : if (lgefint(p)==3 && T)
469 : {
470 332 : ulong pp = p[2];
471 332 : GEN a4p = ZX_to_Flx(a4,pp), a6p = ZX_to_Flx(a6,pp);
472 332 : GEN hp = ZXX_to_FlxX(h, pp, get_FpX_var(T)), Tp = ZXT_to_FlxT(T , pp);
473 332 : res = Flxq_elldivpolmod(a4p, a6p, n, hp, Tp, pp);
474 332 : return gerepileupto(ltop, FlxX_to_ZXX(res));
475 : }
476 3826 : Fq_elldivpolmod_init(&d, a4, a6, n, h, T, p);
477 3826 : res = gcopy(divpol(d.t,d.r2,n,d.E,d.ff));
478 3826 : divpol_free(d.t);
479 3826 : 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 20958 : Fq_ellyn(struct divpolmod_red *d, long k)
492 : {
493 20958 : pari_sp av = avma;
494 20958 : void *E = d->E;
495 20958 : const struct bb_algebra *ff = d->ff;
496 20958 : if (k==1) return mkvec2(ff->one(E), ff->one(E));
497 : else
498 : {
499 16284 : GEN t = d->t, r2 = d->r2;
500 16284 : GEN pn2 = divpol(t,r2,k-2,E,ff);
501 16284 : GEN pp2 = divpol(t,r2,k+2,E,ff);
502 16284 : GEN pn12 = divpol_f2(t,r2,k-1,E,ff);
503 16284 : GEN pp12 = divpol_f2(t,r2,k+1,E,ff);
504 16284 : GEN on = ff->red(E,ff->sub(E, ff->mul(E,pp2,pn12), ff->mul(E,pn2,pp12)));
505 16284 : GEN f = divpol(t,r2,k,E,ff);
506 16284 : GEN f2 = divpol_f2(t,r2,k,E,ff);
507 16284 : GEN f3 = ff->mul(E,f,f2);
508 16284 : if (!odd(k)) f3 = ff->mul(E,f3,r2);
509 16284 : return gerepilecopy(av,mkvec2(on, f3));
510 : }
511 : }
512 :
513 : static void
514 5538 : Fq_elldivpolmod_close(struct divpolmod_red *d)
515 5538 : { divpol_free(d->t); }
516 : static GEN
517 1320 : Fq_elldivpol2(GEN a4, GEN a6, GEN T, GEN p)
518 1320 : { return mkpoln(4, utoi(4), gen_0, Fq_mulu(a4, 4, T, p), Fq_mulu(a6, 4, T, p)); }
519 :
520 : static GEN
521 1320 : Fq_elldivpol2d(GEN a4, GEN T, GEN p)
522 1320 : { return mkpoln(3, utoi(6), gen_0, Fq_mulu(a4, 2, T, p)); }
523 :
524 : static GEN
525 1320 : 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 1320 : long m = degpol(h);
529 1320 : mp1 = gel(h, m + 1); /* negative of first power sum */
530 1320 : dh = FqX_deriv(h, T, p);
531 1320 : ddh = FqX_deriv(dh, T, p);
532 1320 : t = Fq_elldivpol2(a4, a6, T, p);
533 1320 : u = Fq_elldivpol2d(a4, T, p);
534 1320 : t1 = FqX_sub(FqX_sqr(dh, T, p), FqX_mul(ddh, h, T, p), T, p);
535 1320 : t2 = FqX_mul(u, FqX_mul(h, dh, T, p), T, p);
536 1320 : t3 = FqX_mul(FqX_sqr(h, T, p),
537 : deg1pol_shallow(stoi(2*m), Fq_mulu(mp1, 2, T, p), vx), T, p);
538 1320 : f0 = FqX_add(FqX_sub(FqX_mul(t, t1, T, p), t2, T, p), t3, T, p);
539 1320 : t4 = FqX_mul(pol_x(vx), FqX_sqr(h, T, p), T, p);
540 1320 : return FqX_add(t4, f0, T, p);
541 : }
542 :
543 : static GEN
544 888 : Zq_inv(GEN b, GEN T, GEN p, long e)
545 : {
546 1734 : return e==1 ? Fq_inv(b, T, p):
547 846 : typ(b)==t_INT ? Zp_inv(b, p, e): ZpXQ_inv(b, T, p, e);
548 : }
549 :
550 : static GEN
551 84378 : Zq_div(GEN a, GEN b, GEN T, GEN q, GEN p, long e)
552 : {
553 84378 : if (e==1) return Fq_div(a, b, T, p);
554 846 : 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 12 : Zq_divexact(GEN a, GEN b)
566 12 : { return typ(a)==t_INT ? diviiexact(a, b): ZX_Z_divexact(a, b); }
567 :
568 : static long
569 12 : Zq_pval(GEN a, GEN p)
570 12 : { return typ(a)==t_INT ? Z_pval(a, p): ZX_pval(a, p); }
571 :
572 : static GEN
573 103032 : Zq_divu_safe(GEN a, ulong b, GEN T, GEN q, GEN p, long e)
574 : {
575 : long v, w;
576 103032 : if (e==1) return Fq_div(a, utoi(b), T, q);
577 2238 : v = u_pvalrem(b, p, &b);
578 2238 : if (v > 0)
579 : {
580 12 : if (signe(a)==0) return gen_0;
581 12 : w = Zq_pval(a, p);
582 12 : if (v > w) return NULL;
583 12 : a = Zq_divexact(a, powiu(p,v));
584 : }
585 2238 : return Fq_Fp_mul(a, Fp_inv(utoi(b), q), T, q);
586 : }
587 :
588 : static GEN
589 140898 : FqX_shift(GEN P,long n)
590 140898 : { return RgX_shift_shallow(P, n); }
591 :
592 : static GEN
593 33276 : FqX_mulhigh_i(GEN f, GEN g, long n, GEN T, GEN p)
594 33276 : { return FqX_shift(FqX_mul(f,g,T, p),-n); }
595 :
596 : static GEN
597 33276 : FqX_mulhigh(GEN f, GEN g, long n2, long n, GEN T, GEN p)
598 : {
599 33276 : GEN F = RgX_blocks(f, n2, 2), fl = gel(F,1), fh = gel(F,2);
600 33276 : 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 16638 : FqX_invlift1(GEN Q, GEN P, long t1, long t2, GEN T, GEN p)
605 : {
606 16638 : GEN H = FqXn_mul(FqX_mulhigh(Q, P, t1, t2, T, p), Q, t2-t1, T, p);
607 16638 : return FqX_sub(Q, FqX_shift(H, t1), T, p);
608 : }
609 :
610 : static GEN
611 16638 : FqX_invsqrtlift1(GEN Q, GEN P, long t1, long t2, GEN T, GEN p)
612 : {
613 16638 : GEN D = FqX_mulhigh(P, FqX_sqr(Q, T, p), t1, t2, T, p);
614 16638 : GEN H = FqXn_mul(Q, FqX_halve(D, T, p), t2-t1, T, p);
615 16638 : 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 22746 : ZqX_integ2Xn(GEN P, long N, GEN T, GEN p, GEN pp, long e)
621 : {
622 22746 : long d = degpol(P), v = varn(P);
623 : long k;
624 : GEN Q;
625 22746 : if(d==-1) return pol_0(v);
626 16638 : Q = cgetg(d+3,t_POL);
627 16638 : Q[1] = evalsigne(1) | evalvarn(v);
628 71208 : for (k = 0; k <= d; k++)
629 : {
630 54570 : GEN q = Zq_divu_safe(gel(P,2+k), 2*(k+N)+1, T, p, pp, e);
631 54570 : if (!q) return NULL;
632 54570 : gel(Q, 2+k) = q;
633 : }
634 16638 : return ZXX_renormalize(Q,d+3);
635 : }
636 :
637 : /* solution of G*(S'^2)=(S/x)*(HoS) mod x^m */
638 : static GEN
639 6108 : Zq_Weierstrass(GEN a4, GEN a6, GEN b4, GEN b6, long m, GEN T, GEN p, GEN pp, long n)
640 : {
641 6108 : pari_sp av = avma;
642 6108 : long v = 0;
643 6108 : ulong mask = quadratic_prec_mask(m);
644 6108 : GEN iGdS2 = pol_1(v);
645 6108 : GEN G = mkpoln(4, a6, a4, gen_0, gen_1);
646 6108 : GEN GdS2 = G, S = pol_x(v), sG = pol_1(v), isG = sG, dS = sG;
647 6108 : long N = 1;
648 22746 : for (;mask>1;)
649 : {
650 : GEN S2, HS, K, dK, E;
651 22746 : long N2 = N, d;
652 22746 : N<<=1; if (mask & 1) N--;
653 22746 : mask >>= 1;
654 22746 : d = N-N2;
655 22746 : S2 = FqX_sqr(S, T, p);
656 22746 : HS = FqX_Fq_add(FqX_Fq_mul(S, b6, T, p), b4, T, p);
657 22746 : HS = FqX_Fq_add(FqXn_mul(S2, HS, N, T, p), gen_1, T, p);
658 22746 : HS = FqXn_mul(HS, FqX_shift(S,-1), N, T, p);
659 22746 : sG = FqXn_mul(G, isG, N2, T, p);
660 : /* (HS-Gds2)/(Gds2*sG) */
661 22746 : dK = FqXn_mul(FqX_shift(FqX_sub(HS, GdS2, T, p), -N2),
662 : FqXn_mul(iGdS2, isG, d, T, p), d, T, p);
663 22746 : K = ZqX_integ2Xn(dK, N2, T, p, pp, n);
664 22746 : if (!K) return gc_NULL(av);
665 22746 : E = FqXn_mul(FqXn_mul(K, sG, d, T, p), dS, d, T, p);
666 22746 : S = FqX_add(S, FqX_shift(E, N2+1), T, p);
667 22746 : if (mask <= 1) break;
668 16638 : isG = FqX_invsqrtlift1(isG, G, N2, N, T, p);
669 16638 : dS = FqX_deriv(S, T, p);
670 16638 : GdS2 = FqX_mul(G, FqX_sqr(dS, T, p), T, p);
671 16638 : iGdS2 = FqX_invlift1(iGdS2, GdS2, N2, N, T, p);
672 : }
673 6108 : return gerepileupto(av, S);
674 : }
675 :
676 : static GEN
677 6108 : ZqXn_WNewton(GEN S, long l, GEN a4, GEN a6, GEN pp1, GEN T, GEN p, GEN pp, long e)
678 : {
679 6108 : long d = degpol(S);
680 : long k;
681 6108 : GEN Ge = cgetg(2+d,t_POL);
682 6108 : Ge[1] = evalsigne(1);
683 6108 : gel(Ge,2) = pp1;
684 6108 : if (d >= 2)
685 : {
686 6108 : 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 6108 : if (!g) return NULL;
688 6108 : gel(Ge, 3) = g;
689 : }
690 6108 : if (d >= 3)
691 : {
692 6108 : 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 6108 : Fq_mulu(a6,(l-1)*2,T,p),T,p),10,T,p,pp,e);
695 6108 : if (!g) return NULL;
696 6108 : gel(Ge, 4) = g;
697 : }
698 42354 : for (k = 4; k <= d; k++)
699 : {
700 72492 : GEN g = Zq_divu_safe(Fq_sub(Fq_sub(gel(S,4+k-2),
701 36246 : Fq_mul(a4,Fq_mulu(gel(Ge,k-1),4*k-6,T,p),T,p),T,p),
702 36246 : Fq_mul(a6,Fq_mulu(gel(Ge,k-2),4*k-8,T,p),T,p),T,p),
703 36246 : 4*k-2, T, p, pp, e);
704 36246 : if (!g) return NULL;
705 36246 : gel(Ge, k+1) = g;
706 : }
707 6108 : return ZXX_renormalize(Ge, 2+d);
708 : }
709 :
710 : /****************************************************************************/
711 : /* SIMPLE ELLIPTIC CURVE OVER Fq */
712 : /****************************************************************************/
713 :
714 : static GEN
715 2232 : Fq_ellj(GEN a4, GEN a6, GEN T, GEN p)
716 : {
717 2232 : pari_sp ltop=avma;
718 2232 : GEN a43 = Fq_mulu(Fq_powu(a4, 3, T, p), 4, T, p);
719 2232 : 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 2232 : return gerepileupto(ltop, j);
722 : }
723 :
724 : static GEN
725 2304 : Zq_ellj(GEN a4, GEN a6, GEN T, GEN p, GEN pp, long e)
726 : {
727 2304 : pari_sp ltop=avma;
728 2304 : GEN a43 = Fq_mulu(Fq_powu(a4, 3, T, p), 4, T, p);
729 2304 : 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 2304 : return gerepileupto(ltop, j);
732 : }
733 : /****************************************************************************/
734 : /* EIGENVALUE */
735 : /****************************************************************************/
736 :
737 : static GEN
738 316 : Flxq_find_eigen_Frobenius(GEN a4, GEN a6, GEN h, GEN T, ulong p)
739 : {
740 316 : long v = get_FlxqX_var(h), vT = get_Flx_var(T);
741 316 : GEN RHS = FlxqX_rem(Flxq_rhs(a4, a6, v, vT), h, T, p);
742 316 : return FlxqXQ_halfFrobenius(RHS, h, T, p);
743 : }
744 :
745 : static GEN
746 5222 : Fq_find_eigen_Frobenius(GEN a4, GEN a6, GEN h, GEN T, GEN p)
747 : {
748 5222 : long v = T ? get_FpXQX_var(h): get_FpX_var(h);
749 5222 : GEN RHS = FqX_rem(rhs(a4, a6, v), h, T, p);
750 10238 : return T ? FpXQXQ_halfFrobenius(RHS, h, T, p):
751 5016 : 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 432 : find_eigen_value_oneroot(GEN a4, GEN a6, ulong ell, GEN tr, GEN h, GEN T, GEN p)
758 : {
759 432 : pari_sp ltop = avma;
760 : ulong t;
761 : struct divpolmod_red d;
762 : GEN f, Dy, Gy;
763 432 : h = FqX_get_red(h, T, p);
764 432 : Gy = Fq_find_eigen_Frobenius(a4, a6, h, T, p);
765 432 : t = Fl_div(tr[1], 2, ell);
766 432 : if (t < (ell>>1)) t = ell - t;
767 432 : Fq_elldivpolmod_init(&d, a4, a6, t, h, T, p);
768 432 : f = Fq_ellyn(&d, t);
769 432 : Dy = FqXQ_mul(Gy, gel(f,2), h, T, p);
770 432 : if (!gequal(gel(f,1), Dy)) t = ell-t;
771 432 : Fq_elldivpolmod_close(&d);
772 432 : return gc_ulong(ltop, t);
773 : }
774 :
775 : static ulong
776 316 : Flxq_find_eigen_value_power(GEN a4, GEN a6, ulong ell, long k, ulong lambda,
777 : GEN h, GEN T, ulong p)
778 : {
779 316 : pari_sp ltop = avma;
780 316 : ulong t, ellk1 = upowuu(ell, k-1), ellk = ell*ellk1;
781 : pari_timer ti;
782 : struct divpolmod_red d;
783 : GEN Gy;
784 316 : timer_start(&ti);
785 316 : h = FlxqX_get_red(h, T, p);
786 316 : Gy = Flxq_find_eigen_Frobenius(a4, a6, h, T, p);
787 316 : if (DEBUGLEVEL>2) err_printf(" (%ld ms)",timer_delay(&ti));
788 316 : Flxq_elldivpolmod_init(&d, a4, a6, ellk, h, T, p);
789 1437 : for (t = lambda; t < ellk; t += ellk1)
790 : {
791 1437 : GEN f = Fq_ellyn(&d, t);
792 1437 : GEN Dr = FlxqXQ_mul(Gy, gel(f,2), h, T, p);
793 1437 : if (varn(gel(f,1))!=varn(Dr)) pari_err_BUG("find_eigen_value_power");
794 1437 : if (gequal(gel(f,1), Dr)) break;
795 1229 : if (gequal(gel(f,1), FlxX_neg(Dr,p))) { t = ellk-t; break; }
796 : }
797 316 : if (DEBUGLEVEL>2) err_printf(" (%ld ms)",timer_delay(&ti));
798 316 : Fq_elldivpolmod_close(&d);
799 316 : 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 4790 : Fq_find_eigen_value_power(GEN a4, GEN a6, ulong ell, long k, ulong lambda, GEN h, GEN T, GEN p)
806 : {
807 4790 : pari_sp ltop = avma;
808 4790 : ulong t, ellk1 = upowuu(ell, k-1), ellk = ell*ellk1;
809 : pari_timer ti;
810 : struct divpolmod_red d;
811 : GEN Gy;
812 4790 : timer_start(&ti);
813 4790 : h = FqX_get_red(h, T, p);
814 4790 : Gy = Fq_find_eigen_Frobenius(a4, a6, h, T, p);
815 4790 : if (DEBUGLEVEL>2) err_printf(" (%ld ms)",timer_delay(&ti));
816 4790 : Fq_elldivpolmod_init(&d, a4, a6, ellk, h, T, p);
817 19089 : for (t = lambda; t < ellk; t += ellk1)
818 : {
819 19089 : GEN f = Fq_ellyn(&d, t);
820 19089 : GEN Dr = FqXQ_mul(Gy, gel(f,2), h, T, p);
821 19089 : if (varn(gel(f,1))!=varn(Dr)) pari_err_BUG("find_eigen_value_power");
822 19089 : if (gequal(gel(f,1), Dr)) break;
823 15289 : if (gequal(gel(f,1), FqX_neg(Dr,T,p))) { t = ellk-t; break; }
824 : }
825 4790 : if (DEBUGLEVEL>2) err_printf(" (%ld ms)",timer_delay(&ti));
826 4790 : Fq_elldivpolmod_close(&d);
827 4790 : return gc_ulong(ltop, t);
828 : }
829 :
830 : static ulong
831 5106 : find_eigen_value_power(GEN a4, GEN a6, ulong ell, long k, ulong lambda, GEN hq, GEN T, GEN p)
832 : {
833 5106 : ulong pp = itou_or_0(p);
834 5106 : if (pp && T)
835 : {
836 316 : GEN a4p = ZX_to_Flx(a4, pp);
837 316 : GEN a6p = ZX_to_Flx(a6, pp);
838 316 : GEN hp = ZXXT_to_FlxXT(hq, pp,varn(a4));
839 316 : GEN Tp = ZXT_to_FlxT(T, pp);
840 316 : return Flxq_find_eigen_value_power(a4p, a6p, ell, k, lambda, hp, Tp, pp);
841 : }
842 4790 : return Fq_find_eigen_value_power(a4, a6, ell, k, lambda, hq, T, p);
843 : }
844 :
845 : static GEN
846 7662 : 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 7662 : long d = ((l+1)>>1)+1;
850 7662 : if(l==3)
851 1554 : return deg1pol(gen_1, Fq_neg(pp1, T, p), 0);
852 6108 : S = Zq_Weierstrass(a4, a6, b4, b6, d + 1, T, p, pp, e);
853 6108 : if (S==NULL) return NULL;
854 6108 : S = FqX_shift(S, -1);
855 6108 : Sd = FqXn_inv(S, d, T, p);
856 6108 : Ge = ZqXn_WNewton(Sd, l, a4, a6, pp1, T, p, pp, e);
857 6108 : if (!Ge) return NULL;
858 6108 : Ge = FqX_neg(Ge, T, p);
859 612 : Ge = T && lgefint(pp)==3 ? ZlXQXn_expint(Ge, d, T, p, pp[2])
860 6263 : : FqXn_expint(Ge, d, T, p);
861 6108 : Ge = RgX_recip(FqX_red(Ge, T, pp));
862 6108 : if (degpol(Ge)==(l-1)>>1) return Ge;
863 1254 : return NULL;
864 : }
865 :
866 : static GEN
867 5634 : 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 5634 : pari_sp ltop = avma;
870 5634 : GEN dxxgj = FqX_eval(Dxxg, j, T, p);
871 5634 : GEN dxJgj = FqX_eval(DxJg, j, T, p);
872 5634 : GEN dJJgj = FqX_eval(DJJg, j, T, p);
873 5634 : GEN E42 = Fq_sqr(E4, T, p), E6ovE4 = Zq_div(E6, E4, T, p, pp, e);
874 5634 : GEN a = Fq_mul(gprime, dxxgj, T, p);
875 5634 : GEN b = Fq_mul(Fq_mul(Fq_mulu(j,2*q, T, p), dxJgj, T, p), E6ovE4, T, p);
876 5634 : GEN c = Fq_mul(Zq_div(Fq_sqr(E6ovE4, T, p), gprime, T, p, pp, e), j, T, p);
877 5634 : 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 5634 : 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 5634 : GEN g = Fq_sub(Fq_sub(b,a,T,p), d, T, p);
881 5634 : 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 7620 : a4a6t(GEN *a4t, GEN *a6t, ulong l, GEN E4t, GEN E6t, GEN T, GEN p)
886 : {
887 7620 : GEN l2 = modii(sqru(l), p), l4 = Fp_sqr(l2, p), l6 = Fp_mul(l4, l2, p);
888 7620 : *a4t = Fq_mul(E4t, Fp_muls(l4, -3, p), T, p);
889 7620 : *a6t = Fq_mul(E6t, Fp_muls(l6, -2, p), T, p);
890 7620 : }
891 : static void
892 42 : a4a6t_from_J(GEN *a4t, GEN *a6t, ulong l, GEN C4t, GEN C6t, GEN T, GEN p)
893 : {
894 42 : GEN l2 = modii(sqru(l), p), l4 = Fp_sqr(l2, p), l6 = Fp_mul(l4, l2, p);
895 42 : GEN v = Fp_inv(stoi(-864), p), u = Fp_mulu(v, 18, p);
896 42 : *a4t = Fq_mul(C4t, Fp_mul(u, l4, p), T, p);
897 42 : *a6t = Fq_mul(C6t, Fp_mul(v, l6, p), T, p);
898 42 : }
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 2208 : find_isogenous_from_Atkin(GEN a4, GEN a6, ulong ell, struct meqn *MEQN, GEN g, GEN T, GEN pp, long e)
905 : {
906 2208 : pari_sp ltop = avma, btop;
907 2208 : GEN meqn = MEQN->eq, meqnx, Dmeqnx, Roots, gprime, u1;
908 2208 : long k, vJ = MEQN->vy;
909 2208 : GEN p = e==1 ? pp: powiu(pp, e);
910 2208 : GEN j = Zq_ellj(a4, a6, T, p, pp, e);
911 2208 : GEN E4 = Fq_div(a4, stoi(-3), T, p);
912 2208 : GEN E6 = Fq_neg(Fq_halve(a6, T, p), T, p);
913 2208 : GEN Dx = RgX_deriv(meqn);
914 2208 : GEN DJ = deriv(meqn, vJ);
915 2208 : GEN Dxg = FpXY_Fq_evaly(Dx, g, T, p, vJ);
916 2208 : GEN px = FqX_eval(Dxg, j, T, p), dx = Fq_mul(px, g, T, p);
917 2208 : GEN DJg = FpXY_Fq_evaly(DJ, g, T, p, vJ);
918 2208 : GEN pJ = FqX_eval(DJg, j, T, p), dJ = Fq_mul(pJ, j, T, p);
919 2208 : GEN Dxx = RgX_deriv(Dx);
920 2208 : GEN DxJg = FqX_deriv(Dxg, T, p);
921 :
922 2208 : GEN Dxxg = FpXY_Fq_evaly(Dxx, g, T, p, vJ);
923 2208 : GEN DJJg = FqX_deriv(DJg, T, p);
924 : GEN a, b;
925 2208 : if (!signe(Fq_red(dJ,T,pp)) || !signe(Fq_red(dx,T,pp)))
926 : {
927 18 : if (DEBUGLEVEL>0) err_printf("[A: d%c=0]",signe(dJ)? 'x': 'J');
928 18 : return gc_NULL(ltop);
929 : }
930 2190 : a = Fq_mul(dJ, Fq_mul(g, E6, T, p), T, p);
931 2190 : b = Fq_mul(E4, dx, T, p);
932 2190 : gprime = Zq_div(a, b, T, p, pp, e);
933 :
934 2190 : u1 = compute_u(gprime, Dxxg, DxJg, DJJg, j, pJ, px, 1, E4, E6, T, p, pp, e);
935 2190 : meqnx = FpXY_Fq_evaly(meqn, g, T, p, vJ);
936 2190 : Dmeqnx = FqX_deriv(meqnx, T, pp);
937 2190 : Roots = FqX_roots(meqnx, T, pp);
938 :
939 2190 : btop = avma;
940 3456 : for (k = lg(Roots)-1; k >= 1; k--, set_avma(btop))
941 : {
942 3456 : GEN jt = gel(Roots, k);
943 3456 : if (signe(FqX_eval(Dmeqnx, jt, T, pp))==0)
944 0 : continue;
945 3456 : if (e > 1)
946 78 : jt = ZqX_liftroot(meqnx, gel(Roots, k), T, pp, e);
947 3456 : if (signe(Fq_red(jt, T, pp)) == 0 || signe(Fq_sub(jt, utoi(1728), T, pp)) == 0)
948 : {
949 12 : if (DEBUGLEVEL>0) err_printf("[A: jt=%ld]",signe(Fq_red(jt,T,p))? 1728: 0);
950 12 : return gc_NULL(ltop);
951 : }
952 : else
953 : {
954 3444 : GEN pxstar = FqX_eval(Dxg, jt, T, p);
955 3444 : GEN dxstar = Fq_mul(pxstar, g, T, p);
956 3444 : GEN pJstar = FqX_eval(DJg, jt, T, p);
957 3444 : GEN dJstar = Fq_mul(Fq_mulu(jt, ell, T, p), pJstar, T, p);
958 3444 : GEN u = Fq_mul(Fq_mul(dxstar, dJ, T, p), E6, T, p);
959 3444 : GEN v = Fq_mul(Fq_mul(dJstar, dx, T, p), E4, T, p);
960 3444 : 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 3444 : GEN E6t = Zq_div(Fq_mul(u, E4t, T, p), v, T, p, pp, e);
962 3444 : GEN u2 = compute_u(gprime, Dxxg, DxJg, DJJg, jt, pJstar, pxstar, ell, E4t, E6t, T, p, pp, e);
963 3444 : GEN pp1 = Fq_mulu(Fq_sub(u1, u2, T, p), 3*ell, T, p);
964 : GEN a4t, a6t, h;
965 3444 : a4a6t(&a4t, &a6t, ell, E4t, E6t, T, p);
966 3444 : h = find_kernel(a4, a6, ell, a4t, a6t, pp1, T, p, pp, e);
967 3444 : if (h && signe(Fq_elldivpolmod(a4, a6, ell, h, T, pp))==0)
968 2178 : 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 4182 : find_isogenous_from_canonical(GEN a4, GEN a6, ulong ell, struct meqn *MEQN, GEN g, GEN T, GEN pp, long e)
981 : {
982 4182 : pari_sp ltop = avma;
983 4182 : GEN meqn = MEQN->eq;
984 4182 : long vJ = MEQN->vy;
985 4182 : GEN p = e==1 ? pp: powiu(pp, e);
986 : GEN h;
987 4182 : GEN E4 = Fq_div(a4, stoi(-3), T, p);
988 4182 : GEN E6 = Fq_neg(Fq_halve(a6, T, p), T, p);
989 4182 : GEN E42 = Fq_sqr(E4, T, p);
990 4182 : GEN E43 = Fq_mul(E4, E42, T, p);
991 4182 : GEN E62 = Fq_sqr(E6, T, p);
992 4182 : GEN delta = Fq_div(Fq_sub(E43, E62, T, p), utoi(1728), T, p);
993 4182 : GEN j = Zq_div(E43, delta, T, p, pp, e);
994 4182 : GEN Dx = RgX_deriv(meqn);
995 4182 : GEN DJ = deriv(meqn, vJ);
996 4182 : GEN Dxg = FpXY_Fq_evaly(Dx, g, T, p, vJ);
997 4182 : GEN px = FqX_eval(Dxg, j, T, p), dx = Fq_mul(px, g, T, p);
998 4182 : GEN DJg = FpXY_Fq_evaly(DJ, g, T, p, vJ);
999 4182 : GEN pJ = FqX_eval(DJg, j, T, p), dJ = Fq_mul(j, pJ, T, p);
1000 4182 : GEN Dxx = RgX_deriv(Dx);
1001 4182 : GEN DxJg = FqX_deriv(Dxg, T, p);
1002 :
1003 4182 : GEN ExJ = FqX_eval(DxJg, j, T, p);
1004 4182 : ulong tis = ugcd(12, ell-1), is = 12 / tis;
1005 4182 : GEN itis = Fq_inv(stoi(-tis), T, p);
1006 4182 : 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 4182 : 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 4182 : 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 4182 : 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 4182 : GEN gd = Fq_mul(Fq_mul(E2s, itis, T, p), g, T, p);
1029 4182 : GEN jd = Zq_div(Fq_mul(Fq_neg(E42, T, p), E6, T, p), delta, T, p, pp, e);
1030 4182 : GEN E0b = Zq_div(E6, Fq_mul(E4, E2s, T, p), T, p, pp, e);
1031 4182 : GEN Dxxgj = FqXY_eval(Dxx, g, j, T, p);
1032 4182 : 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 4182 : GEN DJgJj = FqX_eval(FqX_deriv(DJg, T, p), j, T, p);
1034 4182 : 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 4182 : 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 4182 : 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 4182 : jl = Zq_div(Fq_powu(E4l, 3, T, p), deltal, T, p, pp, e);
1038 4182 : if (signe(Fq_red(jl,T,pp))==0)
1039 : {
1040 6 : if (DEBUGLEVEL>0) err_printf("[C: jl=0]");
1041 6 : return gc_NULL(ltop);
1042 : }
1043 4176 : f = Zq_div(powuu(ell, is), g, T, p, pp, e);
1044 4176 : fd = Fq_neg(Fq_mul(Fq_mul(E2s, f, T, p), itis, T, p), T, p);
1045 4176 : Dgs = FqXY_eval(Dx, f, jl, T, p);
1046 4176 : Djs = FqXY_eval(DJ, f, jl, T, p);
1047 4176 : 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 4176 : E6l = Zq_div(Fq_mul(Fq_neg(E4l, T, p), jld, T, p), jl, T, p, pp, e);
1050 4176 : p_1 = Fq_neg(Fq_halve(Fq_mulu(E2s, ell, T, p), T, p),T,p);
1051 : }
1052 4176 : a4a6t(&a4t, &a6t, ell, E4l, E6l, T, p);
1053 4176 : h = find_kernel(a4, a6, ell, a4t, a6t, p_1, T, p, pp, e);
1054 4176 : if (!h) return NULL;
1055 4176 : return gerepilecopy(ltop, mkvec3(a4t, a6t, h));
1056 : }
1057 :
1058 : static GEN
1059 84 : corr(GEN c4, GEN c6, GEN T, GEN p, GEN pp, long e)
1060 : {
1061 84 : GEN c46 = Zq_div(Fq_sqr(c4, T, p), c6, T, p, pp, e);
1062 84 : GEN c64 = Zq_div(c6, c4, T, p, pp, e);
1063 84 : GEN a = Fp_divu(gen_2, 3, p);
1064 84 : return Fq_add(Fq_halve(c46, T, p), Fq_mul(a, c64, T, p), T, p);
1065 : }
1066 :
1067 : static GEN
1068 144 : RgXY_deflatex(GEN H, long n, long d)
1069 : {
1070 144 : long i, l = lg(H);
1071 144 : GEN R = cgetg(l, t_POL);
1072 144 : R[1] = H[1];
1073 840 : for(i = 2; i < l; i++)
1074 : {
1075 696 : GEN Hi = gel(H, i);
1076 696 : gel(R,i) = typ(Hi)==t_POL? RgX_deflate(RgX_shift_shallow(Hi, d), n): Hi;
1077 : }
1078 144 : return RgX_renormalize_lg(R, l);
1079 : }
1080 :
1081 : static GEN
1082 60 : Fq_polmodular_eval(GEN meqn, GEN j, long N, GEN T, GEN p, long vJ)
1083 : {
1084 60 : pari_sp av = avma;
1085 : GEN R, dR, ddR;
1086 60 : long t0 = N%3 == 1 ? 2: 0;
1087 60 : long t2 = N%3 == 1 ? 0: 2;
1088 60 : if (N == 3)
1089 : {
1090 12 : GEN P = FpXX_red(meqn, p);
1091 12 : GEN dP = deriv(P, -1), ddP = deriv(dP, -1);
1092 12 : R = FpXY_Fq_evaly(P, j, T, p, vJ);
1093 12 : dR = FpXY_Fq_evaly(dP, j, T, p, vJ);
1094 12 : ddR = FpXY_Fq_evaly(ddP, j, T, p, vJ);
1095 12 : return gerepilecopy(av, mkvec3(R,dR,ddR));
1096 : }
1097 : else
1098 : {
1099 48 : GEN P5 = FpXX_red(meqn, p);
1100 48 : GEN H = RgX_splitting(P5, 3);
1101 48 : GEN H0 = RgXY_deflatex(gel(H,1), 3, -t0);
1102 48 : GEN H1 = RgXY_deflatex(gel(H,2), 3, -1);
1103 48 : GEN H2 = RgXY_deflatex(gel(H,3), 3, -t2);
1104 48 : GEN h0 = FpXY_Fq_evaly(H0, j, T, p, vJ);
1105 48 : GEN h1 = FpXY_Fq_evaly(H1, j, T, p, vJ);
1106 48 : GEN h2 = FpXY_Fq_evaly(H2, j, T, p, vJ);
1107 48 : GEN dH0 = RgX_deriv(H0);
1108 48 : GEN dH1 = RgX_deriv(H1);
1109 48 : GEN dH2 = RgX_deriv(H2);
1110 48 : GEN ddH0 = RgX_deriv(dH0);
1111 48 : GEN ddH1 = RgX_deriv(dH1);
1112 48 : GEN ddH2 = RgX_deriv(dH2);
1113 48 : GEN d0 = FpXY_Fq_evaly(dH0, j, T, p, vJ);
1114 48 : GEN d1 = FpXY_Fq_evaly(dH1, j, T, p, vJ);
1115 48 : GEN d2 = FpXY_Fq_evaly(dH2, j, T, p, vJ);
1116 48 : GEN dd0 = FpXY_Fq_evaly(ddH0, j, T, p, vJ);
1117 48 : GEN dd1 = FpXY_Fq_evaly(ddH1, j, T, p, vJ);
1118 48 : 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 48 : h02 = FqX_sqr(h0, T, p);
1123 48 : h12 = FqX_sqr(h1, T, p);
1124 48 : h22 = FqX_sqr(h2, T, p);
1125 48 : h03 = FqX_mul(h0, h02, T, p);
1126 48 : h13 = FqX_mul(h1, h12, T, p);
1127 48 : h23 = FqX_mul(h2, h22, T, p);
1128 48 : h012 = FqX_mul(FqX_mul(h0, h1, T, p), h2, T, p);
1129 48 : dh03 = FqX_mul(FqX_mulu(d0, 3, T, p), h02, T, p);
1130 48 : dh13 = FqX_mul(FqX_mulu(d1, 3, T, p), h12, T, p);
1131 48 : dh23 = FqX_mul(FqX_mulu(d2, 3, T, p), h22, T, p);
1132 48 : 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 48 : R1 = FqX_sub(h13, FqX_mulu(h012, 3, T, p), T, p);
1134 48 : 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 48 : dR1 = FqX_sub(dh13, FqX_mulu(dh012, 3, T, p), T, p);
1136 48 : 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 48 : 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 48 : 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 48 : 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 48 : 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 48 : ddR1 = FqX_sub(ddh13, FqX_mulu(ddh012, 3, T, p), T, p);
1142 48 : 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 48 : 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 48 : return gerepilecopy(av, mkvec3(R ,dR ,ddR));
1145 : }
1146 : }
1147 :
1148 : static GEN
1149 9948 : meqn_j(struct meqn *MEQN, GEN j, long ell, GEN T, GEN p)
1150 : {
1151 9948 : if (MEQN->type=='J')
1152 : {
1153 60 : MEQN->eval = Fq_polmodular_eval(MEQN->eq, j, ell, T, p, MEQN->vy);
1154 60 : return gel(MEQN->eval, 1);
1155 : }
1156 : else
1157 9888 : return FqXY_evalx(MEQN->eq, j, T, p);
1158 : }
1159 :
1160 : static GEN
1161 42 : find_isogenous_from_J(GEN a4, GEN a6, ulong ell, struct meqn *MEQN, GEN g, GEN T, GEN pp, long e)
1162 : {
1163 42 : pari_sp ltop = avma;
1164 42 : GEN meqn = MEQN->eval;
1165 42 : 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 42 : 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 42 : C4 = Fq_mul(a4, stoi(-48), T, p);
1178 42 : C6 = Fq_mul(a6, stoi(-864), T, p);
1179 42 : 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 42 : j = Zq_ellj(a4, a6, T, p, pp, e);
1185 42 : jp = Fq_mul(j, Zq_div(C6, C4, T, p, pp, e), T, p);
1186 42 : co = corr(C4, C6, T, p, pp, e);
1187 42 : Py = RgX_deriv(gel(meqn, 1));
1188 42 : Pxy = RgX_deriv(gel(meqn,2));
1189 42 : Pyy = RgX_deriv(Py);
1190 42 : Pxj = FqX_eval(gel(meqn, 2), g, T, p);
1191 42 : if (signe(Pxj)==0)
1192 : {
1193 0 : if (DEBUGLEVEL>0) err_printf("[J: Pxj=0]");
1194 0 : return gc_NULL(ltop);
1195 : }
1196 42 : Pyj = FqX_eval(Py, g, T, p);
1197 42 : Pxxj = FqX_eval(gel(meqn, 3), g, T, p);
1198 42 : Pxyj = FqX_eval(Pxy, g, T, p);
1199 42 : Pyyj = FqX_eval(Pyy, g, T, p);
1200 42 : jtp = Fq_div(Fq_mul(jp, Zq_div(Pxj, Pyj, T, p, pp, e), T, p),
1201 : utoineg(ell), T, p);
1202 42 : jtp2 = Fq_sqr(jtp,T,p);
1203 42 : jtp3 = Fq_mul(jtp,jtp2,T,p);
1204 42 : den = Fq_mul(Fq_sqr(g,T,p),Fq_sub(g,utoi(1728),T,p),T, p);
1205 42 : D = Zq_inv(den, T, pp, e);
1206 42 : C4t = Fq_mul(jtp2,Fq_mul(g, D, T, p), T, p);
1207 42 : C6t = Fq_mul(jtp3, D, T, p);
1208 42 : s0 = Fq_mul(Fq_sqr(jp, T, p), Pxxj, T, p);
1209 42 : s1 = Fq_mul(Fq_mulu(Fq_mul(jp,jtp,T,p),2*ell,T,p), Pxyj, T, p);
1210 42 : s2 = Fq_mul(Fq_mulu(jtp2,ell*ell,T,p), Pyyj, T, p);
1211 42 : 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 42 : cot = corr(C4t, C6t, T, p, pp, e);
1213 42 : c0 = Fq_sub(co,Fq_mulu(cot,ell,T,p),T,p);
1214 42 : p_1 = Fq_div(Fq_mulu(Fq_add(s3, c0, T, p),ell,T,p),stoi(-4),T,p);
1215 42 : a4a6t_from_J(&a4t, &a6t, ell, C4t, C6t, T, p);
1216 42 : h = find_kernel(a4, a6, ell, a4t, a6t, p_1, T, p, pp, e);
1217 42 : if (!h) return NULL;
1218 42 : return gerepilecopy(ltop, mkvec3(a4t, a6t, h));
1219 : }
1220 :
1221 : static GEN
1222 6438 : find_isogenous(GEN a4,GEN a6, ulong ell, struct meqn *MEQN, GEN g, GEN T,GEN p)
1223 : {
1224 6438 : ulong pp = itou_or_0(p);
1225 6438 : long e = pp ? ulogint(((ell+1)>>1)+1, pp) + ulogint(2*ell+4, pp) + 1: 1;
1226 6438 : if (signe(a4)==0 || signe(a6)==0)
1227 : {
1228 6 : if (DEBUGLEVEL>0) err_printf("[%c: j=%ld]",MEQN->type,signe(a4)==0 ?0: 1728);
1229 6 : return NULL;
1230 : }
1231 6432 : if (e > 1)
1232 : {
1233 54 : GEN pe = powiu(p, e);
1234 54 : GEN meqnj = meqn_j(MEQN, Zq_ellj(a4, a6, T, pe, p, e), ell, T, pe);
1235 54 : g = ZqX_liftroot(meqnj, g, T, p, e);
1236 : }
1237 6432 : switch(MEQN->type)
1238 : {
1239 4182 : case 'C': return find_isogenous_from_canonical(a4,a6,ell, MEQN, g, T,p,e);
1240 2208 : case 'A': return find_isogenous_from_Atkin(a4,a6,ell, MEQN, g, T,p,e);
1241 42 : default: return find_isogenous_from_J(a4,a6,ell, MEQN, g, T,p,e);
1242 : }
1243 : }
1244 :
1245 : static GEN
1246 5298 : FqX_homogenous_eval(GEN P, GEN A, GEN B, GEN T, GEN p)
1247 : {
1248 5298 : long d = degpol(P), i, v = varn(A);
1249 5298 : GEN s = scalar_ZX_shallow(gel(P, d+2), v), Bn = pol_1(v);
1250 17532 : for (i = d-1; i >= 0; i--)
1251 : {
1252 12234 : Bn = FqX_mul(Bn, B, T, p);
1253 12234 : s = FqX_add(FqX_mul(s, A, T, p), FqX_Fq_mul(Bn, gel(P,i+2), T, p), T, p);
1254 : }
1255 5298 : return s;
1256 : }
1257 :
1258 : static GEN
1259 1110 : FqX_homogenous_div(GEN P, GEN Q, GEN A, GEN B, GEN T, GEN p)
1260 : {
1261 1110 : GEN z = cgetg(3, t_RFRAC);
1262 1110 : long d = degpol(Q)-degpol(P);
1263 1110 : gel(z, 1) = FqX_homogenous_eval(P, A, B, T, p);
1264 1110 : gel(z, 2) = FqX_homogenous_eval(Q, A, B, T, p);
1265 1110 : if (d > 0)
1266 0 : gel(z, 1) = FqX_mul(gel(z, 1), FqX_powu(B, d, T, p), T, p);
1267 1110 : else if (d < 0)
1268 1110 : gel(z, 2) = FqX_mul(gel(z, 2), FqX_powu(B, -d, T, p), T, p);
1269 1110 : return z;
1270 : }
1271 :
1272 : static GEN
1273 1320 : 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 1320 : pari_sp ltop = avma, btop;
1276 : GEN a4t, a6t, gtmp;
1277 1320 : GEN num_iso = FqX_numer_isog_abscissa(kpoly, Eba4, Eba6, T, p, 0);
1278 1320 : GEN mpoly = meqn_j(MEQN, Fq_ellj(Eca4, Eca6, T, p), ell, T, p);
1279 1320 : GEN mroots = FqX_roots(mpoly, T, p);
1280 1320 : GEN kpoly2 = FqX_sqr(kpoly, T, p);
1281 1320 : long i, l1 = lg(mroots);
1282 1320 : btop = avma;
1283 2178 : for (i = 1; i < l1; i++)
1284 : {
1285 : GEN h;
1286 1974 : GEN tmp = find_isogenous(Eca4, Eca6, ell, MEQN, gel(mroots, i), T, p);
1287 1974 : if (!tmp) return gc_NULL(ltop);
1288 1968 : a4t = gel(tmp, 1);
1289 1968 : a6t = gel(tmp, 2);
1290 1968 : gtmp = gel(tmp, 3);
1291 :
1292 : /*check that the kernel kpoly is the good one */
1293 1968 : h = FqX_homogenous_eval(gtmp, num_iso, kpoly2, T, p);
1294 1968 : if (signe(Fq_elldivpolmod(Eba4, Eba6, ell, h, T, p)))
1295 : {
1296 1110 : GEN Ic = FqX_homogenous_div(num_iso,kpoly2, numer_i(Ib),denom_i(Ib), T,p);
1297 1110 : GEN kpoly_new = FqX_homogenous_eval(gtmp, numer_i(Ic),denom_i(Ic), T,p);
1298 1110 : return gerepilecopy(ltop, mkvecn(5, a4t, a6t, kpoly_new, gtmp, Ic));
1299 : }
1300 858 : set_avma(btop);
1301 : }
1302 204 : 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 578 : Flxq_study_eqn(GEN mpoly, GEN T, ulong p, long *pt_dG, long *pt_r)
1312 : {
1313 578 : GEN Xq = FlxqX_Frobenius(mpoly, T, p);
1314 578 : GEN G = FlxqX_gcd(FlxX_sub(Xq, pol_x(0), p), mpoly, T, p);
1315 578 : *pt_dG = degpol(G);
1316 578 : if (!*pt_dG) { *pt_r = FlxqX_ddf_degree(mpoly, Xq, T, p); return NULL; }
1317 349 : return gel(FlxqX_roots(G, T, p), 1);
1318 : }
1319 :
1320 : static GEN
1321 7704 : Fp_study_eqn(GEN mpoly, GEN p, long *pt_dG, long *pt_r)
1322 : {
1323 7704 : GEN T = FpX_get_red(mpoly, p);
1324 7704 : GEN XP = FpX_Frobenius(T, p);
1325 7704 : GEN G = FpX_gcd(FpX_sub(XP, pol_x(0), p), mpoly, p);
1326 7704 : *pt_dG = degpol(G);
1327 7704 : if (!*pt_dG) { *pt_r = FpX_ddf_degree(T, XP, p); return NULL; }
1328 4056 : return FpX_oneroot(G, p);
1329 : }
1330 :
1331 : static GEN
1332 8562 : Fq_study_eqn(GEN mpoly, GEN T, GEN p, long *pt_dG, long *pt_r)
1333 : {
1334 : GEN G;
1335 8562 : if (!T) return Fp_study_eqn(mpoly, p, pt_dG, pt_r);
1336 858 : if (lgefint(p)==3)
1337 : {
1338 578 : ulong pp = p[2];
1339 578 : GEN Tp = ZXT_to_FlxT(T,pp);
1340 578 : GEN mpolyp = ZXX_to_FlxX(mpoly,pp,get_FpX_var(T));
1341 578 : G = Flxq_study_eqn(mpolyp, Tp, pp, pt_dG, pt_r);
1342 578 : return G ? Flx_to_ZX(G): NULL;
1343 : }
1344 : else
1345 : {
1346 280 : GEN Xq = FpXQX_Frobenius(mpoly, T, p);
1347 280 : G = FpXQX_gcd(FpXX_sub(Xq, pol_x(0), p), mpoly, T, p);
1348 280 : *pt_dG = degpol(G);
1349 280 : if (!*pt_dG) { *pt_r = FpXQX_ddf_degree(mpoly, Xq, T, p); return NULL; }
1350 119 : return gel(FpXQX_roots(G, T, p), 1);
1351 : }
1352 : }
1353 :
1354 : /* Berlekamp variant */
1355 : static GEN
1356 8574 : study_modular_eqn(long ell, GEN mpoly, GEN T, GEN p, enum mod_type *mt, long *ptr_r)
1357 : {
1358 8574 : pari_sp ltop = avma;
1359 8574 : GEN g = gen_0;
1360 8574 : *ptr_r = 0; /*gcc -Wall*/
1361 8574 : if (!FqX_is_squarefree(mpoly, T, p)) *mt = MTcm;
1362 : else
1363 : {
1364 : long dG;
1365 8562 : g = Fq_study_eqn(mpoly, T, p, &dG, ptr_r);
1366 8562 : switch(dG)
1367 : {
1368 4038 : case 0: *mt = MTAtkin; break;
1369 462 : case 1: *mt = MTone_root; break;
1370 4002 : case 2: *mt = MTElkies; break;
1371 60 : default: *mt = (dG == ell + 1)? MTroots: MTpathological;
1372 : }
1373 : }
1374 8574 : 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 8574 : 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 4464 : 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 4464 : pari_sp ltop = avma, btop;
1391 : GEN tmp, Eba4, Eba6, Eca4, Eca6, Ib, kpoly;
1392 4464 : long k = *pt_k;
1393 4464 : ulong lambda, ellk = upowuu(ell, k), pellk = umodiu(q, ellk);
1394 : long cnt;
1395 :
1396 4464 : if (DEBUGLEVEL) { err_printf("mod %ld", ell); }
1397 4464 : Eba4 = a4;
1398 4464 : Eba6 = a6;
1399 4464 : tmp = find_isogenous(a4,a6, ell, MEQN, g, T, p);
1400 4464 : if (!tmp) return gc_NULL(ltop);
1401 4428 : Eca4 = gel(tmp, 1);
1402 4428 : Eca6 = gel(tmp, 2);
1403 4428 : kpoly = gel(tmp, 3);
1404 4428 : Ib = pol_x(0);
1405 4428 : lambda = tr ? find_eigen_value_oneroot(a4, a6, ell, tr, kpoly, T, p):
1406 3996 : find_eigen_value_power(a4, a6, ell, 1, 1, kpoly, T, p);
1407 4428 : if (DEBUGLEVEL>1) err_printf(" [%ld ms]", timer_delay(ti));
1408 4428 : if (smallfact && smallfact%(long)ell!=0)
1409 : {
1410 324 : ulong pell = pellk%ell;
1411 324 : ulong ap = Fl_add(lambda, Fl_div(pell, lambda, ell), ell);
1412 324 : if (Fl_sub(pell, ap, ell)==ell-1) { set_avma(ltop); return mkvecsmall(ap); }
1413 312 : if (smallfact < 0 && Fl_add(pell, ap, ell)==ell-1) { set_avma(ltop); return mkvecsmall(ap); }
1414 : }
1415 4404 : btop = avma;
1416 5514 : for (cnt = 2; cnt <= k; cnt++)
1417 : {
1418 1320 : GEN tmp = find_kernel_power(Eba4, Eba6, Eca4, Eca6, ell, MEQN, kpoly, Ib, T, p);
1419 1320 : if (!tmp) { k = cnt-1; break; }
1420 1110 : if (DEBUGLEVEL) err_printf(", %Ps", powuu(ell, cnt));
1421 1110 : lambda = find_eigen_value_power(a4, a6, ell, cnt, lambda, gel(tmp,3), T, p);
1422 1110 : Eba4 = Eca4;
1423 1110 : Eba6 = Eca6;
1424 1110 : Eca4 = gel(tmp,1);
1425 1110 : Eca6 = gel(tmp,2);
1426 1110 : kpoly = gel(tmp,4);
1427 1110 : Ib = gel(tmp, 5);
1428 1110 : 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 1110 : if (DEBUGLEVEL>1) err_printf(" [%ld ms]", timer_delay(ti));
1434 : }
1435 4404 : set_avma(ltop);
1436 4404 : ellk = upowuu(ell, k);
1437 4404 : pellk = umodiu(q, ellk);
1438 4404 : *pt_k = k;
1439 4404 : 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 4038 : find_trace_Atkin(ulong ell, long r, GEN q)
1446 : {
1447 4038 : pari_sp ltop = avma;
1448 4038 : long nval = 0;
1449 4038 : ulong teta, pell = umodiu(q, ell), invp = Fl_inv(pell, ell);
1450 4038 : GEN val_pos = cgetg(1+ell, t_VECSMALL), P = gel(factoru(r), 1);
1451 4038 : GEN S = mkvecsmall4(0, pell, 0, 1);
1452 4038 : GEN U = mkvecsmall3(0, ell-1, 0);
1453 4038 : pari_sp btop = avma;
1454 4038 : if (r==2 && krouu(ell-pell, ell) < 0)
1455 606 : val_pos[++nval] = 0;
1456 78942 : for (teta = 1; teta < ell; teta++, set_avma(btop))
1457 : {
1458 74904 : ulong disc = Fl_sub(Fl_sqr(teta,ell), Fl_mul(4UL,pell,ell), ell);
1459 : GEN a;
1460 74904 : if (krouu(disc, ell) >= 0) continue;
1461 36996 : S[3] = Fl_neg(teta, ell);
1462 36996 : U[3] = Fl_mul(invp, teta, ell);
1463 36996 : a = Flxq_powu(U, r/P[1], S, ell);
1464 36996 : if (!Flx_equal1(a) && Flx_equal1(Flxq_powu(a, P[1], S, ell)))
1465 : {
1466 25080 : pari_sp av = avma;
1467 25080 : long i, l=lg(P);
1468 42792 : for (i = 2; i < l; i++, set_avma(av))
1469 22500 : if (Flx_equal1(Flxq_powu(U, r/P[i], S, ell))) break;
1470 25080 : if (i==l) val_pos[++nval] = teta;
1471 : }
1472 : }
1473 4038 : 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 462 : find_trace_one_root(ulong ell, GEN q)
1479 : {
1480 462 : ulong a = Fl_double(Fl_sqrt(umodiu(q,ell), ell), ell);
1481 462 : return mkvecsmall2(a, ell - a);
1482 : }
1483 :
1484 : static GEN
1485 60 : find_trace_lp1_roots(long ell, GEN q)
1486 : {
1487 60 : ulong ell2 = ell * ell, pell = umodiu(q, ell2);
1488 60 : ulong a = Fl_sqrt(pell%ell, ell);
1489 60 : ulong pa = Fl_add(Fl_div(pell, a, ell2), a, ell2);
1490 60 : return mkvecsmall2(pa, ell2 - pa);
1491 : }
1492 :
1493 : /*ell odd prime; trace modulo ell^k: [], [t] or [t1,...,td] */
1494 : static GEN
1495 8574 : 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 8574 : 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 8574 : kt = maxss((long)(log(expi(q)*M_LN2)/log((double)ell)), 1);
1506 8574 : if (DEBUGLEVEL)
1507 0 : { err_printf("SEA: Prime %5ld ", ell); timer_start(&ti); }
1508 8574 : get_modular_eqn(&MEQN, ell, vx, vy);
1509 8574 : meqnj = meqn_j(&MEQN, j, ell, T, p);
1510 8574 : 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 8574 : switch (mt)
1515 : {
1516 462 : case MTone_root:
1517 462 : tr2 = find_trace_one_root(ell, q);
1518 462 : tr = find_trace_Elkies_power(a4,a6,ell, &kt, &MEQN, g, tr2, q, T, p, smallfact, &ti);
1519 462 : if (!tr) { tr = tr2; kt = 1; }
1520 462 : break;
1521 4002 : case MTElkies:
1522 : /* Contrary to MTone_root, may look mod higher powers of ell */
1523 4002 : if (abscmpiu(p, 2*ell+3) <= 0)
1524 42 : kt = 1; /* Not implemented in this case */
1525 4002 : tr = find_trace_Elkies_power(a4,a6,ell, &kt, &MEQN, g, NULL, q, T, p, smallfact, &ti);
1526 4002 : if (!tr)
1527 : {
1528 6 : if (DEBUGLEVEL) err_printf("[fail]\n");
1529 6 : return gc_NULL(ltop);
1530 : }
1531 3996 : break;
1532 60 : case MTroots:
1533 60 : tr = find_trace_lp1_roots(ell, q);
1534 60 : kt = 2;
1535 60 : break;
1536 4038 : case MTAtkin:
1537 4038 : tr = find_trace_Atkin(ell, r, q);
1538 4038 : if (lg(tr)==1) pari_err_PRIME("ellap",p);
1539 4038 : kt = 1;
1540 4038 : break;
1541 12 : case MTcm:
1542 : {
1543 12 : long D = find_CM(ell, j, T, p);
1544 12 : GEN C = Fq_ellcard_CM(D, a4, a6, T, p);
1545 12 : if (DEBUGLEVEL>1) err_printf(" D=%ld [%ld ms]\n", D, timer_delay(&ti));
1546 12 : return gc_const(ltop, C);
1547 : }
1548 0 : default: /* case MTpathological: */
1549 0 : return gc_NULL(ltop);
1550 : }
1551 8556 : 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 8556 : *ptr_kt = kt;
1561 8556 : 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 780 : separation(GEN cnt)
1571 : {
1572 : pari_sp btop;
1573 780 : long k = lg(cnt)-1, l = (1L<<k)-1, best_i, i, j;
1574 : GEN best_r, P, P3, r;
1575 :
1576 780 : P = gen_1;
1577 3900 : for (j = 1; j <= k; ++j) P = mulis(P, cnt[j]);
1578 : /* p_b * p_g = P is constant */
1579 780 : P3 = mulsi(3, P);
1580 780 : btop = avma;
1581 780 : best_i = 0;
1582 780 : best_r = P3;
1583 37956 : for (i = 1; i < l; i++)
1584 : {
1585 : /* scan all possibilities */
1586 37254 : GEN p_b = gen_1;
1587 356526 : for (j = 0; j < k; j++)
1588 319272 : if (i & (1L<<j)) p_b = mulis(p_b, cnt[1+j]);
1589 37254 : r = subii(shifti(sqri(p_b), 2), P3); /* (p_b/p_g - 3/4)*4*P */
1590 37254 : if (!signe(r)) { best_i = i; break; }
1591 37176 : if (abscmpii(r, best_r) < 0) { best_i = i; best_r = r; }
1592 37176 : if (gc_needed(btop, 1))
1593 0 : best_r = gerepileuptoint(btop, best_r);
1594 : }
1595 780 : 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 1560 : multiple_crt(GEN x, GEN y, GEN q, GEN P)
1604 : {
1605 1560 : pari_sp ltop = avma, av;
1606 1560 : long i, j, k, lx = lg(x)-1, ly = lg(y)-1;
1607 : GEN a1, a2, u, v, A2X;
1608 1560 : (void)bezout(P,q,&u,&v);
1609 1560 : a1 = mulii(P,u);
1610 1560 : a2 = mulii(q,v); A2X = ZC_Z_mul(x, a2);
1611 1560 : av = avma; affii(mulii(P,q), P);
1612 62580 : for (i = 1, k = 1; i <= lx; i++, set_avma(av))
1613 : {
1614 61020 : GEN a2x = gel(A2X,i);
1615 1024044 : for (j = 1; j <= ly; ++j)
1616 : {
1617 963024 : GEN t = Fp_add(Fp_mulu(a1, y[j], P), a2x, P);
1618 963024 : affii(t, gel(x, k++));
1619 : }
1620 : }
1621 1560 : setlg(x, k); set_avma(ltop);
1622 1560 : }
1623 :
1624 : /****************************************************************************/
1625 : /* MATCH AND SORT */
1626 : /****************************************************************************/
1627 :
1628 : static GEN
1629 1560 : possible_traces(GEN compile, GEN mask, GEN *P, int larger)
1630 : {
1631 1560 : GEN V, Pfinal = gen_1, C = shallowextract(compile, mask);
1632 1560 : long i, lfinal = 1, lC = lg(C), lP;
1633 1560 : pari_sp av = avma;
1634 :
1635 4680 : for (i = 1; i < lC; i++)
1636 : {
1637 3120 : GEN c = gel(C,i), t;
1638 3120 : Pfinal = mulii(Pfinal, gel(c,1));
1639 3120 : t = muluu(lfinal, lg(gel(c,2))-1);
1640 3120 : lfinal = itou(t);
1641 : }
1642 1560 : Pfinal = gerepileuptoint(av, Pfinal);
1643 1560 : if (larger)
1644 780 : lP = lgefint(shifti(Pfinal,1));
1645 : else
1646 780 : lP = lgefint(Pfinal);
1647 1560 : lfinal++;
1648 : /* allocate room for final result */
1649 1560 : V = cgetg(lfinal, t_VEC);
1650 909756 : for (i = 1; i < lfinal; i++) gel(V,i) = cgeti(lP);
1651 :
1652 : {
1653 1560 : GEN c = gel(C,1), v = gel(c,2);
1654 1560 : long l = lg(v);
1655 7752 : for (i = 1; i < l; i++) affsi(v[i], gel(V,i));
1656 1560 : setlg(V, l); affii(gel(c,1), Pfinal); /* reset Pfinal */
1657 : }
1658 3120 : for (i = 2; i < lC; i++)
1659 : {
1660 1560 : GEN c = gel(C,i);
1661 1560 : multiple_crt(V, gel(c,2), gel(c,1), Pfinal); /* Pfinal updated! */
1662 : }
1663 1560 : *P = Pfinal; return V;
1664 : }
1665 :
1666 : static GEN
1667 393750 : cost(long mask, GEN cost_vec)
1668 : {
1669 393750 : pari_sp ltop = avma;
1670 : long i;
1671 393750 : GEN c = gen_1;
1672 6148998 : for (i = 1; i < lg(cost_vec); i++)
1673 5755248 : if (mask&(1L<<(i-1)))
1674 2551686 : c = mulis(c, cost_vec[i]);
1675 393750 : return gerepileuptoint(ltop, c);
1676 : }
1677 :
1678 : static GEN
1679 317052 : value(long mask, GEN atkin, long k)
1680 : {
1681 317052 : pari_sp ltop = avma;
1682 : long i;
1683 317052 : GEN c = gen_1;
1684 4952250 : for (i = 1; i <= k; i++)
1685 4635198 : if (mask&(1L<<(i-1)))
1686 2045346 : c = mulii(c, gmael(atkin, i, 1));
1687 317052 : return gerepileuptoint(ltop, c);
1688 : }
1689 :
1690 : static void
1691 156528 : set_cost(GEN B, long b, GEN cost_vec, long *pi)
1692 : {
1693 156528 : pari_sp av = avma;
1694 156528 : GEN costb = cost(b, cost_vec);
1695 156528 : long i = *pi;
1696 214692 : while (cmpii(costb, cost(B[i], cost_vec)) < 0) --i;
1697 156528 : B[++i] = b;
1698 156528 : *pi = i; set_avma(av);
1699 156528 : }
1700 :
1701 : static GEN
1702 1650 : get_lgatkin(GEN compile_atkin, long k)
1703 : {
1704 1650 : GEN v = cgetg(k+1, t_VECSMALL);
1705 : long j;
1706 8784 : for (j = 1; j <= k; ++j) v[j] = lg(gmael(compile_atkin, j, 2))-1;
1707 1650 : return v;
1708 : }
1709 :
1710 : static GEN
1711 870 : champion(GEN atkin, long k, GEN bound_champ)
1712 : {
1713 870 : const long two_k = 1L<<k;
1714 870 : pari_sp ltop = avma;
1715 : long i, j, n, i1, i2;
1716 870 : GEN B, Bp, cost_vec, res = NULL;
1717 :
1718 870 : cost_vec = get_lgatkin(atkin, k);
1719 870 : if (k == 1) return mkvec2(gen_1, utoipos(cost_vec[1]));
1720 :
1721 858 : B = zero_zv(two_k);
1722 858 : Bp = zero_zv(two_k);
1723 858 : Bp[2] = 1;
1724 4002 : for (n = 2, j = 2; j <= k; j++)
1725 : {
1726 : long b;
1727 3144 : i = 1;
1728 148644 : for (i1 = 2, i2 = 1; i1 <= n; )
1729 : {
1730 145500 : pari_sp av = avma;
1731 145500 : long b1 = Bp[i1], b2 = Bp[i2]|(1L<<(j-1));
1732 145500 : if (cmpii(value(b1, atkin, k), value(b2, atkin, k)) < 0)
1733 145500 : { b = b1; i1++; } else { b = b2; i2++; }
1734 145500 : set_avma(av);
1735 145500 : set_cost(B, b, cost_vec, &i);
1736 : }
1737 14172 : for ( ; i2 <= n; i2++)
1738 : {
1739 11028 : b = Bp[i2]|(1L<<(j-1));
1740 11028 : set_cost(B, b, cost_vec, &i);
1741 : }
1742 3144 : n = i;
1743 104652 : for (i = 1; i <= n; i++)
1744 101508 : Bp[i] = B[i];
1745 : }
1746 8255202 : for (i = 1; i <= two_k; i++)
1747 8254344 : if (B[i])
1748 : {
1749 22530 : GEN b = cost (B[i], cost_vec);
1750 22530 : GEN v = value(B[i], atkin, k);
1751 22530 : if (cmpii(v, bound_champ) <=0) continue;
1752 4290 : if (res && gcmp(b, gel(res, 2)) >=0) continue;
1753 858 : res = mkvec2(utoi(B[i]), b);
1754 : }
1755 858 : return gerepilecopy(ltop, res);
1756 : }
1757 :
1758 : static GEN
1759 1560 : compute_diff(GEN v)
1760 : {
1761 1560 : long i, l = lg(v) - 1;
1762 1560 : GEN diff = cgetg(l, t_VEC);
1763 908196 : for (i = 1; i < l; i++) gel(diff, i) = subii(gel(v, i+1), gel(v, i));
1764 1560 : return ZV_sort_uniq_shallow(diff);
1765 : }
1766 :
1767 : static int
1768 14808 : cmp_atkin(void*E, GEN a, GEN b)
1769 : {
1770 14808 : long ta=typ(a)==t_INT, tb=typ(b)==t_INT, c;
1771 : (void) E;
1772 14808 : if (ta || tb) return ta-tb;
1773 4860 : c = lg(gel(a,2)) - lg(gel(b,2));
1774 4860 : if (c) return c;
1775 726 : return cmpii(gel(b,1), gel(a,1));
1776 : }
1777 :
1778 : static void
1779 3522 : add_atkin(GEN atkin, GEN trace, long *nb)
1780 : {
1781 3522 : long l = lg(atkin)-1;
1782 3522 : long i, k = gen_search(atkin, trace, NULL, cmp_atkin);
1783 3522 : if (k > 0 || (k = -k) > l) return;
1784 68508 : for (i = l; i > k; i--) gel(atkin,i) = gel(atkin,i-1);
1785 3522 : if (typ(gel(atkin,l))==t_INT) (*nb)++;
1786 3522 : gel(atkin,k) = trace;
1787 : }
1788 :
1789 : /* V = baby / giant, P = Pb / Pg */
1790 : static GEN
1791 1560 : BSGS_pre(GEN *pdiff, GEN V, GEN P, void *E, const struct bb_group *grp)
1792 : {
1793 1560 : GEN diff = compute_diff(V);
1794 1560 : GEN pre = cgetg(lg(diff), t_VEC);
1795 1560 : long i, l = lg(diff);
1796 1560 : 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 33444 : for (i = 2; i < l; i++)
1799 : {
1800 31884 : pari_sp av = avma;
1801 31884 : GEN d = subii(gel(diff, i), gel(diff, i-1));
1802 31884 : GEN Q = grp->mul(E, gel(pre, i-1), grp->pow(E, P, d));
1803 31884 : gel(pre, i) = gerepilecopy(av, Q);
1804 : }
1805 1560 : *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 828 : 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 828 : 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 828 : long best_i, i, lbaby, lgiant, k = lg(compile_atkin)-1;
1820 828 : GEN bound = sqrti(shifti(q, 2)), card;
1821 828 : const long lcard = 100;
1822 828 : long lq = lgefint(q), nbcard;
1823 : pari_timer ti;
1824 :
1825 828 : if (k == 1)
1826 : { /*only one Atkin prime, check the cardinality with random points */
1827 48 : GEN r = gel(compile_atkin, 1), r1 = gel(r,1), r2 = gel(r,2);
1828 48 : long l = lg(r2), j;
1829 48 : GEN card = cgetg(l, t_VEC), Cs2, C, U;
1830 48 : Z_chinese_pre(Mu, r1, &C,&U, NULL);
1831 48 : Cs2 = shifti(C, -1);
1832 324 : for (j = 1, i = 1; i < l; i++)
1833 : {
1834 276 : GEN t = Z_chinese_post(u, stoi(r2[i]), C, U, NULL);
1835 276 : t = Fp_center_i(t, C, Cs2);
1836 276 : if (abscmpii(t, bound) <= 0) gel(card, j++) = subii(pp1, t);
1837 : }
1838 48 : setlg(card, j);
1839 48 : return gen_select_order(card, E, grp);
1840 : }
1841 780 : if (DEBUGLEVEL>=2) timer_start(&ti);
1842 780 : av1 = avma;
1843 780 : best_i = separation( get_lgatkin(compile_atkin, k) );
1844 780 : set_avma(av1);
1845 :
1846 780 : baby = possible_traces(compile_atkin, utoi(best_i), &Mb, 1);
1847 780 : giant = possible_traces(compile_atkin, subiu(int2n(k), best_i+1), &Mg, 0);
1848 780 : lbaby = lg(baby);
1849 780 : lgiant = lg(giant);
1850 780 : den = Fp_inv(Fp_mul(Mu, Mb, Mg), Mg);
1851 780 : av2 = avma;
1852 533820 : for (i = 1; i < lgiant; i++, set_avma(av2))
1853 533040 : affii(Fp_mul(gel(giant,i), den, Mg), gel(giant,i));
1854 780 : ZV_sort_inplace(giant);
1855 780 : Sg = Fp_mul(negi(u), den, Mg);
1856 780 : den = Fp_inv(Fp_mul(Mu, Mg, Mb), Mb);
1857 780 : dec_inf = divii(mulii(Mb,addii(Mg,shifti(Sg,1))), shifti(Mg,1));
1858 780 : togglesign(dec_inf); /* now, dec_inf = ceil(- (Mb/2 + Sg Mb/Mg) ) */
1859 780 : div = mulii(truedivii(dec_inf, Mb), Mb);
1860 780 : av2 = avma;
1861 375936 : for (i = 1; i < lbaby; i++, set_avma(av2))
1862 : {
1863 375156 : GEN b = addii(Fp_mul(Fp_sub(gel(baby,i), u, Mb), den, Mb), div);
1864 375156 : if (cmpii(b, dec_inf) < 0) b = addii(b, Mb);
1865 375156 : affii(b, gel(baby,i));
1866 : }
1867 780 : ZV_sort_inplace(baby);
1868 :
1869 780 : SgMb = mulii(Sg, Mb);
1870 780 : card = cgetg(lcard+1,t_VEC);
1871 78780 : for (i = 1; i <= lcard; i++) gel(card,i) = cgetipos(lq+1);
1872 :
1873 780 : av2 = avma;
1874 780 : MATCH_RESTART:
1875 780 : set_avma(av2);
1876 780 : nbcard = 0;
1877 780 : P = grp->rand(E);
1878 780 : point = grp->pow(E,P, Mu);
1879 780 : Pb = grp->pow(E,point, Mg);
1880 780 : Pg = grp->pow(E,point, Mb);
1881 : /* Precomputation for babies */
1882 780 : 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 780 : table = cgetg(lbaby, t_VECSMALL);
1887 780 : av1 = avma;
1888 : /* (p+1 - u - Mu*Mb*Sg) P - (baby[1]) Pb */
1889 780 : point = grp->pow(E,P, subii(subii(pp1, u), mulii(Mu, addii(SgMb, mulii(Mg, gel(baby,1))))));
1890 780 : table[1] = grp->hash(gel(point,1));
1891 375156 : for (i = 2; i < lbaby; i++)
1892 : {
1893 374376 : GEN d = subii(gel(baby, i), gel(baby, i-1));
1894 374376 : point = grp->mul(E, point, grp->pow(E, gel(pre, ZV_search(diff, d)), gen_m1));
1895 374376 : table[i] = grp->hash(gel(point,1));
1896 374376 : if (gc_needed(av1,3))
1897 : {
1898 16 : if(DEBUGMEM>1) pari_warn(warnmem,"match_and_sort, baby = %ld", i);
1899 16 : point = gerepileupto(av1, point);
1900 : }
1901 : }
1902 780 : set_avma(av1);
1903 : /* Precomputations for giants */
1904 780 : pre = BSGS_pre(&diff, giant, Pg, E, grp);
1905 :
1906 : /* Look for a collision among the x-coordinates */
1907 780 : table_ind = vecsmall_indexsort(table);
1908 780 : table = perm_mul(table,table_ind);
1909 :
1910 780 : av1 = avma;
1911 780 : point = grp->pow(E, Pg, gel(giant, 1));
1912 780 : for (i = 1; ; i++)
1913 532260 : {
1914 : GEN d;
1915 533040 : long h = grp->hash(gel(point, 1));
1916 533040 : long s = zv_search(table, h);
1917 533040 : if (s) {
1918 1560 : while (table[s] == h && s) s--;
1919 1560 : for (s++; s < lbaby && table[s] == h; s++)
1920 : {
1921 780 : GEN B = gel(baby,table_ind[s]), G = gel(giant,i);
1922 780 : GEN GMb = mulii(G, Mb), BMg = mulii(B, Mg);
1923 780 : GEN Be = subii(subii(pp1, u), mulii(Mu, addii(SgMb, BMg)));
1924 780 : GEN Bp = grp->pow(E,P, Be);
1925 : /* p+1 - u - Mu (Sg Mb + GIANT Mb + BABY Mg) */
1926 780 : if (gequal(gel(Bp,1),gel(point,1)))
1927 : {
1928 780 : GEN card1 = subii(Be, mulii(Mu, GMb));
1929 780 : GEN card2 = addii(card1, mulii(mulsi(2,Mu), GMb));
1930 780 : if (abscmpii(subii(pp1, card1), bound) <= 0)
1931 684 : affii(card1, gel(card, ++nbcard));
1932 780 : if (nbcard >= lcard) goto MATCH_RESTART;
1933 780 : if (abscmpii(subii(pp1, card2), bound) <= 0)
1934 420 : affii(card2, gel(card, ++nbcard));
1935 780 : if (nbcard >= lcard) goto MATCH_RESTART;
1936 : }
1937 : }
1938 : }
1939 533040 : if (i==lgiant-1) break;
1940 532260 : d = subii(gel(giant, i+1), gel(giant, i));
1941 532260 : point = grp->mul(E,point, gel(pre, ZV_search(diff, d)));
1942 532260 : if (gc_needed(av1,3))
1943 : {
1944 22 : if(DEBUGMEM>1) pari_warn(warnmem,"match_and_sort, giant = %ld", i);
1945 22 : point = gerepileupto(av1, point);
1946 : }
1947 : }
1948 780 : setlg(card, nbcard+1);
1949 780 : if (DEBUGLEVEL>=2) timer_printf(&ti,"match_and_sort");
1950 780 : return gen_select_order(card, E, grp);
1951 : }
1952 :
1953 : static GEN
1954 870 : get_bound_bsgs(long lp)
1955 : {
1956 : GEN B;
1957 870 : if (lp <= 160)
1958 840 : B = divru(powru(dbltor(1.048), lp), 9);
1959 30 : else if (lp <= 192)
1960 24 : B = divrr(powru(dbltor(1.052), lp), dbltor(16.65));
1961 : else
1962 6 : B = mulrr(powru(dbltor(1.035), minss(lp,307)), dbltor(1.35));
1963 870 : 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 912 : Fq_ellcard_SEA(GEN a4, GEN a6, GEN q, GEN T, GEN p, long smallfact)
1974 : {
1975 912 : const long MAX_ATKIN = 21;
1976 912 : 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 912 : GEN prod_atkin = gen_1, max_traces = gen_0;
1980 : GEN j;
1981 912 : double bound_gr = 1.;
1982 912 : const double growth_factor = 1.26;
1983 : forprime_t TT;
1984 :
1985 912 : j = Fq_ellj(a4, a6, T, p);
1986 912 : 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 12 : : Fp_ellcard(a4, a6, p);
1989 900 : if (Fq_elljissupersingular(j, T, p))
1990 18 : return Fq_ellcard_supersingular(a4, a6, T, p);
1991 : /*First compute the trace modulo 2 */
1992 882 : switch(FqX_nbroots(rhs(a4, a6, 0), T, p))
1993 : {
1994 60 : case 3: /* bonus time: 4 | #E(Fq) = q+1 - t */
1995 60 : i = mod4(q)+1; if (i > 2) i -= 4;
1996 60 : TR_mod = utoipos(4);
1997 60 : TR = stoi(i); break;
1998 438 : case 1:
1999 438 : TR_mod = gen_2;
2000 438 : TR = gen_0; break;
2001 384 : default : /* 0 */
2002 384 : TR_mod = gen_2;
2003 384 : TR = gen_1; break;
2004 : }
2005 882 : if (odd(smallfact) && !mpodd(TR))
2006 : {
2007 12 : if (DEBUGLEVEL) err_printf("Aborting: #E(Fq) divisible by 2\n");
2008 12 : set_avma(ltop); return gen_0;
2009 : }
2010 870 : vy = fetch_var();
2011 870 : 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 870 : u_forprime_init(&TT, 3, ULONG_MAX);
2016 870 : bound = sqrti(shifti(q, 4));
2017 870 : bound_bsgs = get_bound_bsgs(expi(q));
2018 870 : compile_atkin = zerovec(MAX_ATKIN); nb_atkin = 0;
2019 870 : btop = avma;
2020 8580 : while ( (ell = u_forprime_next(&TT)) )
2021 : {
2022 8580 : long ellkt, kt = 1, nbtrace;
2023 : GEN trace_mod;
2024 8586 : if (absequalui(ell, p)) continue;
2025 8574 : trace_mod = find_trace(a4, a6, j, ell, q, T, p, &kt, smallfact, vx,vy);
2026 8574 : if (!trace_mod) continue;
2027 8568 : if (typ(trace_mod)==t_INT)
2028 : {
2029 12 : delete_var(); delete_var();
2030 870 : return gerepileuptoint(ltop, trace_mod);
2031 : }
2032 8556 : nbtrace = lg(trace_mod) - 1;
2033 8556 : ellkt = (long)upowuu(ell, kt);
2034 8556 : if (nbtrace == 1)
2035 : {
2036 5034 : long t_mod_ellkt = trace_mod[1];
2037 5034 : if (smallfact && smallfact%ell!=0)
2038 : { /* does ell divide q + 1 - t ? */
2039 330 : long q_mod_ell_plus_one = umodiu(q,ell) + 1;
2040 330 : ulong card_mod_ell = umodsu(q_mod_ell_plus_one - t_mod_ellkt, ell);
2041 330 : ulong tcard_mod_ell = 1;
2042 330 : if (card_mod_ell && smallfact < 0)
2043 114 : tcard_mod_ell = umodsu(q_mod_ell_plus_one + t_mod_ellkt, ell);
2044 330 : if (!card_mod_ell || !tcard_mod_ell)
2045 : {
2046 24 : if (DEBUGLEVEL)
2047 0 : err_printf("\nAborting: #E%s(Fq) divisible by %ld\n",
2048 : tcard_mod_ell ? "" : "_twist", ell);
2049 24 : delete_var(); delete_var();
2050 24 : return gc_const(ltop, gen_0);
2051 : }
2052 : }
2053 5010 : (void)Z_incremental_CRT(&TR, t_mod_ellkt, &TR_mod, ellkt);
2054 5010 : if (DEBUGLEVEL)
2055 0 : err_printf(", missing %ld bits\n",expi(bound)-expi(TR_mod));
2056 : }
2057 : else
2058 : {
2059 3522 : add_atkin(compile_atkin, mkvec2(utoipos(ellkt), trace_mod), &nb_atkin);
2060 3522 : prod_atkin = value(-1, compile_atkin, nb_atkin);
2061 : }
2062 8532 : if (cmpii(mulii(TR_mod, prod_atkin), bound) > 0)
2063 : {
2064 : GEN bound_tr;
2065 906 : if (!nb_atkin)
2066 : {
2067 6 : delete_var(); delete_var();
2068 6 : return gerepileuptoint(ltop, subii(addiu(q, 1), TR));
2069 : }
2070 900 : bound_tr = mulrr(bound_bsgs, dbltor(bound_gr));
2071 900 : bound_gr *= growth_factor;
2072 900 : if (signe(max_traces))
2073 : {
2074 72 : max_traces = divis(muliu(max_traces,nbtrace), ellkt);
2075 72 : if (DEBUGLEVEL>=3)
2076 0 : err_printf("At least %Ps remaining possibilities.\n",max_traces);
2077 : }
2078 900 : if (cmpir(max_traces, bound_tr) < 0)
2079 : {
2080 870 : GEN bound_atkin = truedivii(bound, TR_mod);
2081 870 : champ = champion(compile_atkin, nb_atkin, bound_atkin);
2082 870 : max_traces = gel(champ,2);
2083 870 : if (DEBUGLEVEL>=2)
2084 0 : err_printf("%Ps remaining possibilities.\n", max_traces);
2085 870 : if (cmpir(max_traces, bound_tr) < 0)
2086 : {
2087 828 : GEN res, cat = shallowextract(compile_atkin, gel(champ,1));
2088 : const struct bb_group *grp;
2089 : void *E;
2090 828 : if (DEBUGLEVEL)
2091 0 : err_printf("Match and sort for %Ps possibilities.\n", max_traces);
2092 828 : delete_var();
2093 828 : delete_var();
2094 828 : grp = get_FqE_group(&E,a4,a6,T,p);
2095 828 : res = match_and_sort(cat, TR_mod, TR, q, E, grp);
2096 828 : return gerepileuptoint(ltop, res);
2097 : }
2098 : }
2099 : }
2100 7698 : 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 834 : Fp_ellcard_SEA(GEN a4, GEN a6, GEN p, long smallfact)
2108 834 : { return Fq_ellcard_SEA(a4, a6, p, NULL, p, smallfact); }
|
