Line data Source code
1 : /* Copyright (C) 2014 The PARI group.
2 :
3 : This file is part of the PARI/GP package.
4 :
5 : PARI/GP is free software; you can redistribute it and/or modify it under the
6 : terms of the GNU General Public License as published by the Free Software
7 : Foundation; either version 2 of the License, or (at your option) any later
8 : version. It is distributed in the hope that it will be useful, but WITHOUT
9 : ANY WARRANTY WHATSOEVER.
10 :
11 : Check the License for details. You should have received a copy of it, along
12 : with the package; see the file 'COPYING'. If not, write to the Free Software
13 : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
14 :
15 : /********************************************************************/
16 : /** **/
17 : /** HYPERELLIPTIC CURVES **/
18 : /** **/
19 : /********************************************************************/
20 : #include "pari.h"
21 : #include "paripriv.h"
22 :
23 : #define DEBUGLEVEL DEBUGLEVEL_hyperell
24 :
25 : /*****************************************************************************
26 : ******* *******
27 : ******* Naive algorithms for genus2 *******
28 : ******* *******
29 : *****************************************************************************/
30 :
31 : static GEN
32 392 : F2x_genus2charpoly_naive(GEN P, GEN Q)
33 : {
34 392 : long a, b = 1, c = 0;
35 392 : GEN T = mkvecsmall2(P[1], 7);
36 392 : GEN PT = F2x_rem(P, T), QT = F2x_rem(Q, T);
37 392 : long q0 = F2x_eval(Q, 0), q1 = F2x_eval(Q, 1);
38 392 : long dP = F2x_degree(P), dQ = F2x_degree(Q);
39 392 : a= dQ<3 ? 0: dP<=5 ? 1: -1;
40 392 : a += (q0? F2x_eval(P, 0)? -1: 1: 0) + (q1? F2x_eval(P, 1)? -1: 1: 0);
41 392 : b += q0 + q1;
42 392 : if (lgpol(QT))
43 308 : c = (F2xq_trace(F2xq_div(PT, F2xq_sqr(QT, T), T), T)==0 ? 1: -1);
44 392 : return mkvecsmalln(6, 0UL, 4UL, 2*a, (b+2*c+a*a)>>1, a, 1UL);
45 : }
46 :
47 : static GEN
48 19544 : Flx_difftable(GEN P, ulong p)
49 : {
50 19544 : long i, n = degpol(P);
51 19544 : GEN V = cgetg(n+2, t_VEC);
52 19544 : gel(V, n+1) = P;
53 132832 : for(i = n; i >= 1; i--)
54 113288 : gel(V, i) = Flx_diff1(gel(V, i+1), p);
55 19544 : return V;
56 : }
57 :
58 : static GEN
59 94234 : Flx_difftable_constant(GEN P, ulong p)
60 : {
61 94234 : long i, n = degpol(P);
62 94234 : GEN V = cgetg(n+2, t_VECSMALL);
63 94234 : uel(V, n+1) = Flx_constant(P);
64 641704 : for(i = n; i >= 1; i--)
65 : {
66 547470 : P = Flx_diff1(P, p);
67 547470 : uel(V, i) = Flx_constant(P);
68 : }
69 94234 : return V;
70 : }
71 : static GEN
72 766955 : FlxV_Fl2_eval_pre(GEN V, GEN x, ulong D, ulong p, ulong pi)
73 : {
74 766955 : long i, n = lg(V)-1;
75 766955 : GEN r = cgetg(n+1, t_VEC);
76 5981661 : for (i = 1; i <= n; i++)
77 5214706 : gel(r, i) = Flx_Fl2_eval_pre(gel(V, i), x, D, p, pi);
78 766955 : return r;
79 : }
80 :
81 : static GEN
82 97965182 : Fl2V_next(GEN V, ulong p)
83 : {
84 97965182 : long i, n = lg(V)-1;
85 97965182 : GEN r = cgetg(n+1, t_VEC);
86 97965182 : gel(r, 1) = gel(V, 1);
87 666067388 : for (i = 2; i <= n; i++)
88 568102206 : gel(r, i) = Flv_add(gel(V, i), gel(V, i-1), p);
89 97965182 : return r;
90 : }
91 :
92 : static GEN
93 19544 : FlxV_constant(GEN x)
94 152376 : { pari_APPLY_long(Flx_constant(gel(x,i))) }
95 :
96 : static GEN
97 19544 : Flx_genus2charpoly_naive(GEN H, ulong p)
98 : {
99 19544 : pari_sp av = avma, av2;
100 19544 : ulong pi = get_Fl_red(p);
101 19544 : ulong i, j, p2 = p>>1, D = 2, e = ((p&2UL) == 0) ? -1 : 1;
102 19544 : long a, b, c = 0, n = degpol(H);
103 19544 : GEN t, d, k = const_vecsmall(p, -1);
104 19544 : k[1] = 0;
105 786499 : for (i=1, j=1; i < p; i += 2, j = Fl_add(j, i, p)) k[j+1] = 1;
106 35679 : while (k[1+D] >= 0) D++;
107 19544 : b = n == 5 ? 0 : 1;
108 19544 : a = b ? k[1+Flx_lead(H)]: 0;
109 19544 : t = Flx_difftable(H, p);
110 19544 : d = FlxV_constant(t);
111 19544 : av2 = avma;
112 1572998 : for (i=0; i < p; i++)
113 : {
114 1553454 : ulong v = uel(d,n+1);
115 1553454 : a += k[1+v];
116 1553454 : b += !!v;
117 1553454 : if (n==6)
118 1241520 : uel(d,7) = Fl_add(uel(d,7), uel(d,6), p);
119 1553454 : uel(d,6) = Fl_add(uel(d,6), uel(d,5), p);
120 1553454 : uel(d,5) = Fl_add(uel(d,5), uel(d,4), p);
121 1553454 : uel(d,4) = Fl_add(uel(d,4), uel(d,3), p);
122 1553454 : uel(d,3) = Fl_add(uel(d,3), uel(d,2), p);
123 1553454 : uel(d,2) = Fl_add(uel(d,2), uel(d,1), p);
124 : }
125 786499 : for (j=1; j <= p2; j++)
126 : {
127 766955 : GEN V = FlxV_Fl2_eval_pre(t, mkvecsmall2(0, j), D, p, pi);
128 766955 : for (i=0;; i++)
129 97965182 : {
130 98732137 : GEN r2 = gel(V, n+1);
131 197464274 : c += uel(r2,2) ?
132 97786633 : (uel(r2,1) ? uel(k,1+Fl2_norm_pre(r2, D, p, pi)): e)
133 196518770 : : !!uel(r2,1);
134 98732137 : if (i == p-1) break;
135 97965182 : V = Fl2V_next(V, p);
136 : }
137 766955 : set_avma(av2);
138 : }
139 19544 : set_avma(av);
140 19544 : return mkvecsmalln(6, 0UL, p*p, a*p, (b+2*c+a*a)>>1, a, 1UL);
141 : }
142 :
143 : static long
144 94234 : Flx_genus2trace_naive(GEN H, ulong p)
145 : {
146 94234 : pari_sp av = avma;
147 : ulong i, j;
148 94234 : long a, n = degpol(H);
149 94234 : GEN k = const_vecsmall(p, -1), d;
150 94234 : k[1] = 0;
151 40651520 : for (i=1, j=1; i < p; i += 2, j = Fl_add(j, i, p))
152 40557286 : k[j+1] = 1;
153 94234 : a = n == 5 ? 0: k[1+Flx_lead(H)];
154 94234 : d = Flx_difftable_constant(H, p);
155 81303040 : for (i=0; i < p; i++)
156 : {
157 81208806 : a += k[1+uel(d,n+1)];
158 81208806 : if (n==6)
159 65853242 : uel(d,7) = Fl_add(uel(d,7), uel(d,6), p);
160 81208806 : uel(d,6) = Fl_add(uel(d,6), uel(d,5), p);
161 81208806 : uel(d,5) = Fl_add(uel(d,5), uel(d,4), p);
162 81208806 : uel(d,4) = Fl_add(uel(d,4), uel(d,3), p);
163 81208806 : uel(d,3) = Fl_add(uel(d,3), uel(d,2), p);
164 81208806 : uel(d,2) = Fl_add(uel(d,2), uel(d,1), p);
165 : }
166 94234 : return gc_long(av, a);
167 : }
168 :
169 : static GEN
170 98336 : dirgenus2(GEN Q, GEN p, long n)
171 : {
172 98336 : pari_sp av = avma;
173 : GEN f;
174 98336 : if (n > 2)
175 4102 : f = RgX_recip(hyperellcharpoly(gmul(Q,gmodulo(gen_1, p))));
176 : else
177 : {
178 94234 : ulong pp = itou(p);
179 94234 : GEN Qp = ZX_to_Flx(Q, pp);
180 94234 : long t = Flx_genus2trace_naive(Qp, pp);
181 94234 : f = deg1pol_shallow(stoi(t), gen_1, 0);
182 : }
183 98336 : return gc_upto(av, RgXn_inv_i(f, n));
184 : }
185 :
186 : GEN
187 8470 : dirgenus2_worker(GEN P, ulong X, GEN Q)
188 : {
189 8470 : pari_sp av = avma;
190 8470 : long i, l = lg(P);
191 8470 : GEN V = cgetg(l, t_VEC);
192 106806 : for(i = 1; i < l; i++)
193 : {
194 98336 : ulong p = uel(P,i);
195 98336 : long d = ulogint(X, p) + 1; /* minimal d such that p^d > X */
196 98336 : gel(V,i) = dirgenus2(Q, utoi(uel(P,i)), d);
197 : }
198 8470 : return gc_GEN(av, mkvec2(P,V));
199 : }
200 :
201 : GEN
202 553 : vecan_genus2(GEN an, long L)
203 : {
204 553 : GEN Q = gel(an,1), bad = gel(an, 2);
205 553 : GEN worker = snm_closure(is_entry("_dirgenus2_worker"), mkvec(Q));
206 553 : return pardireuler(worker, gen_2, stoi(L), NULL, bad);
207 : }
208 :
209 : /* Implementation of Kedlaya Algorithm for counting point on hyperelliptic
210 : curves by Bill Allombert based on a GP script by Bernadette Perrin-Riou.
211 :
212 : References:
213 : Pierrick Gaudry and Nicolas G\"urel
214 : Counting Points in Medium Characteristic Using Kedlaya's Algorithm
215 : Experiment. Math. Volume 12, Number 4 (2003), 395-402.
216 : http://projecteuclid.org/euclid.em/1087568016
217 :
218 : Harrison, M. An extension of Kedlaya's algorithm for hyperelliptic
219 : curves. Journal of Symbolic Computation, 47 (1) (2012), 89-101.
220 : http://arxiv.org/pdf/1006.4206v3.pdf
221 : */
222 :
223 : /* We use the basis of differentials (x^i*dx/y^k) (i=1 to 2*g-1),
224 : with k either 1 or 3, depending on p and d, see Harrison paper */
225 :
226 : static long
227 1764 : get_basis(long p, long d)
228 : {
229 1764 : if (odd(d))
230 868 : return p < d-1 ? 3 : 1;
231 : else
232 896 : return 2*p <= d-2 ? 3 : 1;
233 : }
234 :
235 : static GEN
236 20265 : FpXXQ_red(GEN S, GEN T, GEN p)
237 : {
238 20265 : pari_sp av = avma;
239 20265 : long i, dS = degpol(S);
240 : GEN A, C;
241 20265 : if (signe(S)==0) return pol_0(varn(T));
242 20265 : A = cgetg(dS+3, t_POL);
243 20265 : C = pol_0(varn(T));
244 1520393 : for(i=dS; i>0; i--)
245 : {
246 1500128 : GEN Si = FpX_add(C, gel(S,i+2), p);
247 1500128 : GEN R, Q = FpX_divrem(Si, T, p, &R);
248 1500128 : gel(A,i+2) = R;
249 1500128 : C = Q;
250 : }
251 20265 : gel(A,2) = FpX_add(C, gel(S,2), p);
252 20265 : A[1] = S[1];
253 20265 : return gc_GEN(av, FpXX_renormalize(A,dS+3));
254 : }
255 :
256 : static GEN
257 3402 : FpXXQ_sqr(GEN x, GEN T, GEN p)
258 : {
259 3402 : pari_sp av = avma;
260 3402 : long n = degpol(T);
261 3402 : GEN z = FpX_red(ZXX_sqr_Kronecker(x, n), p);
262 3402 : z = Kronecker_to_ZXX(z, n, varn(T));
263 3402 : return gc_upto(av, FpXXQ_red(z, T, p));
264 : }
265 :
266 : static GEN
267 16863 : FpXXQ_mul(GEN x, GEN y, GEN T, GEN p)
268 : {
269 16863 : pari_sp av = avma;
270 16863 : long n = degpol(T);
271 16863 : GEN z = FpX_red(ZXX_mul_Kronecker(x, y, n), p);
272 16863 : z = Kronecker_to_ZXX(z, n, varn(T));
273 16863 : return gc_upto(av, FpXXQ_red(z, T, p));
274 : }
275 :
276 : static GEN
277 1309 : ZpXXQ_invsqrt(GEN S, GEN T, ulong p, long e)
278 : {
279 1309 : pari_sp av = avma, av2;
280 : ulong mask;
281 1309 : long v = varn(S), n=1;
282 1309 : GEN a = pol_1(v);
283 1309 : if (e <= 1) return gc_GEN(av, a);
284 1309 : mask = quadratic_prec_mask(e);
285 1309 : av2 = avma;
286 4676 : for (;mask>1;)
287 : {
288 : GEN q, q2, q22, f, fq, afq;
289 3367 : long n2 = n;
290 3367 : n<<=1; if (mask & 1) n--;
291 3367 : mask >>= 1;
292 3367 : q = powuu(p,n); q2 = powuu(p,n2);
293 3367 : f = RgX_sub(FpXXQ_mul(FpXX_red(S, q), FpXXQ_sqr(a, T, q), T, q), pol_1(v));
294 3367 : fq = ZXX_Z_divexact(f, q2);
295 3367 : q22 = shifti(addiu(q2,1),-1);
296 3367 : afq = FpXX_Fp_mul(FpXXQ_mul(a, fq, T, q2), q22, q2);
297 3367 : a = RgX_sub(a, ZXX_Z_mul(afq, q2));
298 3367 : if (gc_needed(av2,1))
299 : {
300 0 : if(DEBUGMEM>1) pari_warn(warnmem,"ZpXXQ_invsqrt, e = %ld", n);
301 0 : a = gc_upto(av2, a);
302 : }
303 : }
304 1309 : return gc_upto(av, a);
305 : }
306 :
307 : static GEN
308 1029749 : to_ZX(GEN a, long v) { return typ(a)==t_INT? scalarpol(a,v): a; }
309 :
310 : static void
311 14 : is_sing(GEN H, ulong p)
312 : {
313 14 : pari_err_DOMAIN("hyperellpadicfrobenius","H","is singular at",utoi(p),H);
314 0 : }
315 :
316 : static void
317 1309 : get_UV(GEN *U, GEN *V, GEN T, ulong p, long e)
318 : {
319 1309 : GEN q = powuu(p,e), d;
320 1309 : GEN dT = FpX_deriv(T, q);
321 1309 : GEN R = polresultantext(T, dT);
322 1309 : long v = varn(T);
323 1309 : if (dvdiu(gel(R,3),p)) is_sing(T, p);
324 1309 : d = Zp_inv(gel(R,3), utoi(p), e);
325 1309 : *U = FpX_Fp_mul(FpX_red(to_ZX(gel(R,1),v),q),d,q);
326 1309 : *V = FpX_Fp_mul(FpX_red(to_ZX(gel(R,2),v),q),d,q);
327 1309 : }
328 :
329 : static GEN
330 133847 : frac_to_Fp(GEN a, GEN b, GEN p)
331 : {
332 133847 : GEN d = gcdii(a, b);
333 133847 : return Fp_div(diviiexact(a, d), diviiexact(b, d), p);
334 : }
335 :
336 : static GEN
337 10094 : ZpXXQ_frob(GEN S, GEN U, GEN V, long k, GEN T, ulong p, long e)
338 : {
339 10094 : pari_sp av = avma, av2;
340 10094 : long i, pr = degpol(S), dT = degpol(T), vT = varn(T);
341 10094 : GEN q = powuu(p,e);
342 10094 : GEN Tp = FpX_deriv(T, q), Tp1 = RgX_shift_shallow(Tp, 1);
343 10094 : GEN M = to_ZX(gel(S,pr+2),vT), R;
344 10094 : av2 = avma;
345 987868 : for(i = pr-1; i>=k; i--)
346 : {
347 : GEN A, B, H, Bc;
348 : ulong v, r;
349 977774 : H = FpX_divrem(FpX_mul(V,M,q), T, q, &B);
350 977774 : A = FpX_add(FpX_mul(U,M,q), FpX_mul(H, Tp, q),q);
351 977774 : v = u_lvalrem(2*i+1,p,&r);
352 977774 : Bc = ZX_deriv(B);
353 977774 : Bc = FpX_Fp_mul(ZX_divuexact(Bc,upowuu(p,v)),Fp_divu(gen_2, r, q), q);
354 977774 : M = FpX_add(to_ZX(gel(S,i+2),vT), FpX_add(A, Bc, q), q);
355 977774 : if (gc_needed(av2,1))
356 : {
357 0 : if(DEBUGMEM>1) pari_warn(warnmem,"ZpXXQ_frob, step 1, i = %ld", i);
358 0 : M = gc_upto(av2, M);
359 : }
360 : }
361 10094 : if (degpol(M)<dT-1)
362 5488 : return gc_upto(av, M);
363 4606 : R = RgX_shift_shallow(M,dT-degpol(M)-2);
364 4606 : av2 = avma;
365 237629 : for(i = degpol(M)-dT+2; i>=1; i--)
366 : {
367 : GEN B, c;
368 233023 : R = RgX_shift_shallow(R, 1);
369 233023 : gel(R,2) = gel(M, i+1);
370 233023 : if (degpol(R) < dT) continue;
371 130935 : B = FpX_add(FpX_mulu(T, 2*i, q), Tp1, q);
372 130935 : c = frac_to_Fp(leading_coeff(R), leading_coeff(B), q);
373 130935 : R = FpX_sub(R, FpX_Fp_mul(B, c, q), q);
374 130935 : if (gc_needed(av2,1))
375 : {
376 0 : if(DEBUGMEM>1) pari_warn(warnmem,"ZpXXQ_frob, step 2, i = %ld", i);
377 0 : R = gc_upto(av2, R);
378 : }
379 : }
380 4606 : if (degpol(R)==dT-1)
381 : {
382 2912 : GEN c = frac_to_Fp(leading_coeff(R), leading_coeff(Tp), q);
383 2912 : R = FpX_sub(R, FpX_Fp_mul(Tp, c, q), q);
384 2912 : return gc_upto(av, R);
385 : } else
386 1694 : return gc_GEN(av, R);
387 : }
388 :
389 : static GEN
390 12026 : revdigits(GEN v)
391 : {
392 12026 : long i, n = lg(v)-1;
393 12026 : GEN w = cgetg(n+2, t_POL);
394 12026 : w[1] = evalsigne(1)|evalvarn(0);
395 168784 : for (i=0; i<n; i++)
396 156758 : gel(w,i+2) = gel(v,n-i);
397 12026 : return FpXX_renormalize(w, n+2);
398 : }
399 :
400 : static GEN
401 10094 : diff_red(GEN s, GEN A, long m, GEN T, GEN p)
402 : {
403 10094 : long v, n, vT = varn(T);
404 : GEN Q, sQ, qS;
405 : pari_timer ti;
406 10094 : if (DEBUGLEVEL>1) timer_start(&ti);
407 10094 : Q = revdigits(FpX_digits(A,T,p));
408 10094 : n = degpol(Q);
409 10094 : if (DEBUGLEVEL>1) timer_printf(&ti,"reddigits");
410 10094 : sQ = FpXXQ_mul(s,Q,T,p);
411 10094 : if (DEBUGLEVEL>1) timer_printf(&ti,"redmul");
412 10094 : qS = RgX_shift_shallow(sQ,m-n);
413 10094 : v = ZX_val(sQ);
414 10094 : if (n > m + v)
415 : {
416 4564 : long i, l = n-m-v;
417 4564 : GEN rS = cgetg(l+1,t_VEC);
418 29190 : for (i = l-1; i >=0 ; i--)
419 24626 : gel(rS,i+1) = to_ZX(gel(sQ, 1+v+l-i), vT);
420 4564 : rS = FpXV_FpX_fromdigits(rS,T,p);
421 4564 : gel(qS,2) = FpX_add(FpX_mul(rS, T, p), gel(qS, 2), p);
422 4564 : if (DEBUGLEVEL>1) timer_printf(&ti,"redadd");
423 : }
424 10094 : return qS;
425 : }
426 :
427 : static GEN
428 10094 : ZC_to_padic(GEN C, GEN q)
429 : {
430 10094 : long i, l = lg(C);
431 10094 : GEN V = cgetg(l,t_COL);
432 102914 : for(i = 1; i < l; i++)
433 92820 : gel(V, i) = gadd(gel(C, i), q);
434 10094 : return V;
435 : }
436 :
437 : static GEN
438 1309 : ZM_to_padic(GEN M, GEN q)
439 : {
440 1309 : long i, l = lg(M);
441 1309 : GEN V = cgetg(l,t_MAT);
442 11403 : for(i = 1; i < l; i++)
443 10094 : gel(V, i) = ZC_to_padic(gel(M, i), q);
444 1309 : return V;
445 : }
446 :
447 : static GEN
448 1743 : ZX_to_padic(GEN P, GEN q)
449 : {
450 1743 : long i, l = lg(P);
451 1743 : GEN Q = cgetg(l, t_POL);
452 1743 : Q[1] = P[1];
453 5978 : for (i=2; i<l ;i++)
454 4235 : gel(Q,i) = gadd(gel(P,i), q);
455 1743 : return normalizepol(Q);
456 : }
457 :
458 : static GEN
459 469 : ZXC_to_padic(GEN x, GEN q)
460 2212 : { pari_APPLY_type(t_COL, ZX_to_padic(gel(x, i), q)) }
461 :
462 : static GEN
463 147 : ZXM_to_padic(GEN x, GEN q)
464 616 : { pari_APPLY_same(ZXC_to_padic(gel(x, i), q)) }
465 :
466 : static GEN
467 1309 : ZlX_hyperellpadicfrobenius(GEN H, ulong p, long n)
468 : {
469 1309 : pari_sp av = avma;
470 : long k, N, i, d;
471 : GEN F, s, Q, pN1, U, V;
472 : pari_timer ti;
473 1309 : if (typ(H) != t_POL) pari_err_TYPE("hyperellpadicfrobenius",H);
474 1309 : if (p == 2) is_sing(H, 2);
475 1309 : d = degpol(H);
476 1309 : if (d <= 0)
477 0 : pari_err_CONSTPOL("hyperellpadicfrobenius");
478 1309 : if (n < 1)
479 0 : pari_err_DOMAIN("hyperellpadicfrobenius","n","<", gen_1, utoi(n));
480 1309 : k = get_basis(p, d);
481 1309 : N = n + ulogint(2*n, p) + 1;
482 1309 : pN1 = powuu(p,N+1);
483 1309 : Q = RgX_to_FpX(H, pN1);
484 1309 : if (dvdiu(leading_coeff(Q),p)) is_sing(H, p);
485 1309 : setvarn(Q,1);
486 1309 : if (DEBUGLEVEL>1) timer_start(&ti);
487 1309 : s = revdigits(FpX_digits(RgX_inflate(Q, p), Q, pN1));
488 1309 : if (DEBUGLEVEL>1) timer_printf(&ti,"s1");
489 1309 : s = ZpXXQ_invsqrt(s, Q, p, N);
490 1309 : if (k==3)
491 35 : s = FpXXQ_mul(s, FpXXQ_sqr(s, Q, pN1), Q, pN1);
492 1309 : if (DEBUGLEVEL>1) timer_printf(&ti,"invsqrt");
493 1309 : get_UV(&U, &V, Q, p, N+1);
494 1309 : F = cgetg(d, t_MAT);
495 11403 : for (i = 1; i < d; i++)
496 : {
497 10094 : pari_sp av2 = avma;
498 : GEN M, D;
499 10094 : D = diff_red(s, monomial(utoipos(p),p*i-1,1),(k*p-1)>>1, Q, pN1);
500 10094 : if (DEBUGLEVEL>1) timer_printf(&ti,"red");
501 10094 : M = ZpXXQ_frob(D, U, V, (k-1)>>1, Q, p, N + 1);
502 10094 : if (DEBUGLEVEL>1) timer_printf(&ti,"frob");
503 10094 : gel(F, i) = gc_GEN(av2, RgX_to_RgC(M, d-1));
504 : }
505 1309 : return gc_upto(av, F);
506 : }
507 :
508 : GEN
509 1309 : hyperellpadicfrobenius(GEN H, ulong p, long n)
510 : {
511 1309 : pari_sp av = avma;
512 1309 : GEN M = ZlX_hyperellpadicfrobenius(H, p, n);
513 1309 : GEN q = zeropadic_shallow(utoipos(p),n);
514 1309 : return gc_upto(av, ZM_to_padic(M, q));
515 : }
516 :
517 : INLINE GEN
518 2247 : FpXXX_renormalize(GEN x, long lx) { return ZXX_renormalize(x,lx); }
519 :
520 : static GEN
521 1806 : ZpXQXXQ_red(GEN F, GEN S, GEN T, GEN q, GEN p, long e)
522 : {
523 1806 : pari_sp av = avma;
524 1806 : long i, dF = degpol(F);
525 : GEN A, C;
526 1806 : if (signe(F)==0) return pol_0(varn(S));
527 1806 : A = cgetg(dF+3, t_POL);
528 1806 : C = pol_0(varn(S));
529 96404 : for(i=dF; i>0; i--)
530 : {
531 94598 : GEN Fi = FpXX_add(C, gel(F,i+2), q);
532 94598 : GEN R, Q = ZpXQX_divrem(Fi, S, T, q, p, e, &R);
533 94598 : gel(A,i+2) = R;
534 94598 : C = Q;
535 : }
536 1806 : gel(A,2) = FpXX_add(C, gel(F,2), q);
537 1806 : A[1] = F[1];
538 1806 : return gc_GEN(av, FpXXX_renormalize(A,dF+3));
539 : }
540 :
541 : static GEN
542 448 : ZpXQXXQ_sqr(GEN x, GEN S, GEN T, GEN q, GEN p, long e)
543 : {
544 448 : pari_sp av = avma;
545 : GEN z, kx;
546 448 : long n = degpol(S), vx = varn(S);
547 448 : kx = RgXX_to_Kronecker_var(x, n, vx);
548 448 : z = Kronecker_to_ZXX(FpXQX_sqr(kx, T, q), n, vx);
549 448 : setvarn(z, varn(x));
550 448 : return gc_upto(av, ZpXQXXQ_red(z, S, T, q, p, e));
551 : }
552 :
553 : static GEN
554 1358 : ZpXQXXQ_mul(GEN x, GEN y, GEN S, GEN T, GEN q, GEN p, long e)
555 : {
556 1358 : pari_sp av = avma;
557 : GEN z, kx, ky;
558 1358 : long n = degpol(S), vx = varn(S);
559 1358 : kx = RgXX_to_Kronecker_var(x, n, vx);
560 1358 : ky = RgXX_to_Kronecker_var(y, n, vx);
561 1358 : z = Kronecker_to_ZXX(FpXQX_mul(ky, kx, T, q), n, vx);
562 1358 : setvarn(z, varn(x));
563 1358 : return gc_upto(av, ZpXQXXQ_red(z, S, T, q, p, e));
564 : }
565 :
566 : static GEN
567 441 : FpXXX_red(GEN z, GEN p)
568 : {
569 : GEN res;
570 441 : long i, l = lg(z);
571 441 : res = cgetg(l,t_POL); res[1] = z[1];
572 17388 : for (i=2; i<l; i++)
573 : {
574 16947 : GEN zi = gel(z,i);
575 16947 : if (typ(zi)==t_INT)
576 0 : gel(res,i) = modii(zi,p);
577 : else
578 16947 : gel(res,i) = FpXX_red(zi,p);
579 : }
580 441 : return FpXXX_renormalize(res,lg(res));
581 : }
582 :
583 : static GEN
584 441 : FpXXX_Fp_mul(GEN z, GEN a, GEN p)
585 : {
586 441 : return FpXXX_red(RgX_Rg_mul(z, a), p);
587 : }
588 :
589 : static GEN
590 154 : ZpXQXXQ_invsqrt(GEN F, GEN S, GEN T, ulong p, long e)
591 : {
592 154 : pari_sp av = avma, av2, av3;
593 : ulong mask;
594 154 : long v = varn(F), n=1;
595 : pari_timer ti;
596 154 : GEN a = pol_1(v), pp = utoipos(p);
597 154 : if (DEBUGLEVEL>1) timer_start(&ti);
598 154 : if (e <= 1) return gc_GEN(av, a);
599 154 : mask = quadratic_prec_mask(e);
600 154 : av2 = avma;
601 595 : for (;mask>1;)
602 : {
603 : GEN q, q2, q22, f, fq, afq;
604 441 : long n2 = n;
605 441 : n<<=1; if (mask & 1) n--;
606 441 : mask >>= 1;
607 441 : q = powuu(p,n); q2 = powuu(p,n2);
608 441 : av3 = avma;
609 441 : f = RgX_sub(ZpXQXXQ_mul(F, ZpXQXXQ_sqr(a, S, T, q, pp, n), S, T, q, pp, n), pol_1(v));
610 441 : fq = gc_upto(av3, RgX_Rg_divexact(f, q2));
611 441 : q22 = shifti(addiu(q2,1),-1);
612 441 : afq = FpXXX_Fp_mul(ZpXQXXQ_mul(a, fq, S, T, q2, pp, n2), q22, q2);
613 441 : a = RgX_sub(a, RgX_Rg_mul(afq, q2));
614 441 : if (gc_needed(av2,1))
615 : {
616 0 : if(DEBUGMEM>1) pari_warn(warnmem,"ZpXQXXQ_invsqrt, e = %ld", n);
617 0 : a = gc_upto(av2, a);
618 : }
619 : }
620 154 : return gc_upto(av, a);
621 : }
622 :
623 : static GEN
624 6573 : frac_to_Fq(GEN a, GEN b, GEN T, GEN q, GEN p, long e)
625 : {
626 6573 : GEN d = gcdii(ZX_content(a), ZX_content(b));
627 6573 : return ZpXQ_div(ZX_Z_divexact(a, d), ZX_Z_divexact(b, d), T, q, p, e);
628 : }
629 :
630 : static GEN
631 469 : ZpXQXXQ_frob(GEN F, GEN U, GEN V, long k, GEN S, GEN T, ulong p, long e)
632 : {
633 469 : pari_sp av = avma, av2;
634 469 : long i, pr = degpol(F), dS = degpol(S), v = varn(T);
635 469 : GEN q = powuu(p,e), pp = utoipos(p);
636 469 : GEN Sp = RgX_deriv(S), Sp1 = RgX_shift_shallow(Sp, 1);
637 469 : GEN M = gel(F,pr+2), R;
638 469 : av2 = avma;
639 52311 : for(i = pr-1; i>=k; i--)
640 : {
641 : GEN A, B, H, Bc;
642 : ulong v, r;
643 51842 : H = ZpXQX_divrem(FpXQX_mul(V, M, T, q), S, T, q, utoipos(p), e, &B);
644 51842 : A = FpXX_add(FpXQX_mul(U, M, T, q), FpXQX_mul(H, Sp, T, q),q);
645 51842 : v = u_lvalrem(2*i+1,p,&r);
646 51842 : Bc = RgX_deriv(B);
647 51842 : Bc = FpXX_Fp_mul(ZXX_Z_divexact(Bc,powuu(p,v)), Fp_divu(gen_2, r, q), q);
648 51842 : M = FpXX_add(gel(F,i+2), FpXX_add(A, Bc, q), q);
649 51842 : if (gc_needed(av2,1))
650 : {
651 0 : if(DEBUGMEM>1) pari_warn(warnmem,"ZpXQXXQ_frob, step 1, i = %ld", i);
652 0 : M = gc_upto(av2, M);
653 : }
654 : }
655 469 : if (degpol(M)<dS-1)
656 266 : return gc_upto(av, M);
657 203 : R = RgX_shift_shallow(M,dS-degpol(M)-2);
658 203 : av2 = avma;
659 7175 : for(i = degpol(M)-dS+2; i>=1; i--)
660 : {
661 : GEN B, c;
662 6972 : R = RgX_shift_shallow(R, 1);
663 6972 : gel(R,2) = gel(M, i+1);
664 6972 : if (degpol(R) < dS) continue;
665 6412 : B = FpXX_add(FpXX_mulu(S, 2*i, q), Sp1, q);
666 6412 : c = frac_to_Fq(to_ZX(leading_coeff(R),v), to_ZX(leading_coeff(B),v), T, q, pp, e);
667 6412 : R = FpXX_sub(R, FpXQX_FpXQ_mul(B, c, T, q), q);
668 6412 : if (gc_needed(av2,1))
669 : {
670 0 : if(DEBUGMEM>1) pari_warn(warnmem,"ZpXXQ_frob, step 2, i = %ld", i);
671 0 : R = gc_upto(av2, R);
672 : }
673 : }
674 203 : if (degpol(R)==dS-1)
675 : {
676 161 : GEN c = frac_to_Fq(to_ZX(leading_coeff(R),v), to_ZX(leading_coeff(Sp),v), T, q, pp, e);
677 161 : R = FpXX_sub(R, FpXQX_FpXQ_mul(Sp, c, T, q), q);
678 161 : return gc_upto(av, R);
679 : } else
680 42 : return gc_GEN(av, R);
681 : }
682 :
683 : static GEN
684 469 : Fq_diff_red(GEN s, GEN A, long m, GEN S, GEN T, GEN q, GEN p, long e)
685 : {
686 : long v, n;
687 : GEN Q, sQ, qS;
688 : pari_timer ti;
689 469 : if (DEBUGLEVEL>1) timer_start(&ti);
690 469 : Q = revdigits(ZpXQX_digits(A, S, T, q, p, e));
691 469 : n = degpol(Q);
692 469 : if (DEBUGLEVEL>1) timer_printf(&ti,"reddigits");
693 469 : sQ = ZpXQXXQ_mul(s, Q, S, T, q, p, e);
694 469 : if (DEBUGLEVEL>1) timer_printf(&ti,"redmul");
695 469 : qS = RgX_shift_shallow(sQ,m-n);
696 469 : v = ZX_val(sQ);
697 469 : if (n > m + v)
698 : {
699 189 : long i, l = n-m-v;
700 189 : GEN rS = cgetg(l+1,t_VEC);
701 1547 : for (i = l-1; i >=0 ; i--)
702 1358 : gel(rS,i+1) = gel(sQ, 1+v+l-i);
703 189 : rS = FpXQXV_FpXQX_fromdigits(rS, S, T, q);
704 189 : gel(qS,2) = FpXX_add(FpXQX_mul(rS, S, T, q), gel(qS, 2), q);
705 189 : if (DEBUGLEVEL>1) timer_printf(&ti,"redadd");
706 : }
707 469 : return qS;
708 : }
709 :
710 : static void
711 154 : Fq_get_UV(GEN *U, GEN *V, GEN S, GEN T, ulong p, long e)
712 : {
713 154 : GEN q = powuu(p, e), pp = utoipos(p), d;
714 154 : GEN dS = RgX_deriv(S), R = polresultantext(S, dS), C;
715 154 : long v = varn(S);
716 154 : if (signe(FpX_red(to_ZX(gel(R,3),v), pp))==0) is_sing(S, p);
717 147 : C = FpXQ_red(to_ZX(gel(R, 3),v), T, q);
718 147 : d = ZpXQ_inv(C, T, pp, e);
719 147 : *U = FpXQX_FpXQ_mul(FpXQX_red(to_ZX(gel(R,1),v),T,q),d,T,q);
720 147 : *V = FpXQX_FpXQ_mul(FpXQX_red(to_ZX(gel(R,2),v),T,q),d,T,q);
721 147 : }
722 :
723 : static GEN
724 469 : ZXX_to_FpXC(GEN x, long N, GEN p, long v)
725 : {
726 : long i, l;
727 : GEN z;
728 469 : l = lg(x)-1; x++;
729 469 : if (l > N+1) l = N+1; /* truncate higher degree terms */
730 469 : z = cgetg(N+1,t_COL);
731 2170 : for (i=1; i<l ; i++)
732 : {
733 1701 : GEN xi = gel(x, i);
734 1701 : gel(z,i) = typ(xi)==t_INT? scalarpol(Fp_red(xi, p), v): FpX_red(xi, p);
735 : }
736 511 : for ( ; i<=N ; i++)
737 42 : gel(z,i) = pol_0(v);
738 469 : return z;
739 : }
740 :
741 : GEN
742 154 : ZlXQX_hyperellpadicfrobenius(GEN H, GEN T, ulong p, long n)
743 : {
744 154 : pari_sp av = avma;
745 : long k, N, i, d, N1, v0;
746 : GEN xp, F, s, q, Q, pN1, U, V, pp;
747 : pari_timer ti;
748 154 : if (typ(H) != t_POL) pari_err_TYPE("hyperellpadicfrobenius",H);
749 154 : if (p == 2) is_sing(H, 2);
750 154 : d = degpol(H);
751 154 : if (d <= 0) pari_err_CONSTPOL("hyperellpadicfrobenius");
752 154 : if (n < 1) pari_err_DOMAIN("hyperellpadicfrobenius","n","<", gen_1, utoi(n));
753 154 : k = get_basis(p, d); pp = utoipos(p);
754 154 : N = n + ulogint(2*n, p) + 1;
755 154 : q = powuu(p,n); N1 = N+1;
756 154 : pN1 = powuu(p,N1); T = FpX_get_red(T, pN1);
757 154 : Q = RgX_to_FqX(H, T, pN1);
758 154 : if (signe(FpX_red(to_ZX(leading_coeff(Q),varn(Q)),pp))==0) is_sing(H, p);
759 154 : if (DEBUGLEVEL>1) timer_start(&ti);
760 154 : xp = ZpX_Frobenius(T, pp, N1);
761 154 : s = RgX_inflate(FpXY_FpXQ_evalx(Q, xp, T, pN1), p);
762 154 : v0 = fetch_var_higher();
763 154 : s = revdigits(ZpXQX_digits(s, Q, T, pN1, pp, N1));
764 154 : setvarn(s, v0);
765 154 : if (DEBUGLEVEL>1) timer_printf(&ti,"s1");
766 154 : s = ZpXQXXQ_invsqrt(s, Q, T, p, N);
767 154 : if (k==3)
768 7 : s = ZpXQXXQ_mul(s, ZpXQXXQ_sqr(s, Q, T, pN1, pp, N1), Q, T, pN1, pp, N1);
769 154 : if (DEBUGLEVEL>1) timer_printf(&ti,"invsqrt");
770 154 : Fq_get_UV(&U, &V, Q, T, p, N+1);
771 147 : if (DEBUGLEVEL>1) timer_printf(&ti,"get_UV");
772 147 : F = cgetg(d, t_MAT);
773 616 : for (i = 1; i < d; i++)
774 : {
775 469 : pari_sp av2 = avma;
776 : GEN M, D;
777 469 : D = Fq_diff_red(s, monomial(pp,p*i-1,0),(k*p-1)>>1, Q, T, pN1, pp, N1);
778 469 : if (DEBUGLEVEL>1) timer_printf(&ti,"red");
779 469 : M = ZpXQXXQ_frob(D, U, V, (k - 1)>>1, Q, T, p, N1);
780 469 : if (DEBUGLEVEL>1) timer_printf(&ti,"frob");
781 469 : gel(F, i) = gc_upto(av2, ZXX_to_FpXC(M, d-1, q, varn(T)));
782 : }
783 147 : delete_var();
784 147 : return gc_upto(av, F);
785 : }
786 :
787 : GEN
788 154 : nfhyperellpadicfrobenius(GEN H, GEN T, ulong p, long n)
789 : {
790 154 : pari_sp av = avma;
791 154 : GEN pp = utoipos(p), q = zeropadic_shallow(pp, n);
792 154 : GEN M = ZlXQX_hyperellpadicfrobenius(lift_shallow(H),T,p,n);
793 147 : GEN MM = ZpXQM_prodFrobenius(M, T, pp, n);
794 147 : GEN m = gmul(ZXM_to_padic(MM, q), gmodulo(gen_1, T));
795 147 : return gc_upto(av, m);
796 : }
797 :
798 : GEN
799 595 : hyperellpadicfrobenius0(GEN H, GEN Tp, long n)
800 : {
801 : GEN T, p;
802 595 : if (!ff_parse_Tp(Tp, &T,&p,0)) pari_err_TYPE("hyperellpadicfrobenius", Tp);
803 595 : if (lgefint(p) > 3) pari_err_IMPL("large prime in hyperellpadicfrobenius");
804 7 : return T? nfhyperellpadicfrobenius(H, T, itou(p), n)
805 602 : : hyperellpadicfrobenius(H, itou(p), n);
806 : }
807 :
808 : static GEN
809 679 : charpoly_funceq(GEN P, GEN q)
810 : {
811 679 : long i, l, g = degpol(P)>>1;
812 679 : GEN R, Q = gpowers0(q, g-1, q); /* Q[i] = q^i, i <= g */
813 679 : R = cgetg_copy(P, &l); R[1] = P[1];
814 3164 : for (i=0; i<g; i++) gel(R, i+2) = mulii(gel(P, 2*g-i+2), gel(Q, g-i));
815 3843 : for (; i<=2*g; i++) gel(R, i+2) = icopy(gel(P, i+2));
816 679 : return R;
817 : }
818 :
819 : static long
820 686 : hyperell_Weil_bound(GEN q, ulong g, GEN p)
821 : {
822 686 : pari_sp av = avma;
823 686 : GEN w = mulii(binomialuu(2*g,g),sqrtint(shifti(powiu(q, g),2)));
824 686 : return gc_long(av, logint(w,p) + 1);
825 : }
826 :
827 : /* return 4P + Q^2 */
828 : static GEN
829 409364 : check_hyperell(GEN PQ)
830 : {
831 : GEN H;
832 409364 : if (is_vec_t(typ(PQ)) && lg(PQ)==3)
833 292095 : H = gadd(gsqr(gel(PQ, 2)), gmul2n(gel(PQ, 1), 2));
834 : else
835 117269 : H = gmul2n(PQ, 2);
836 409364 : return typ(H) == t_POL? H: NULL;
837 : }
838 :
839 : static long
840 544674 : hyperellgenus(GEN H)
841 544674 : { long d = degpol(H); return ((d+1)>>1)-1; }
842 :
843 : static void
844 155078 : check_hyperell_Rg(const char *fun, GEN *pW, GEN *pF)
845 : {
846 155078 : GEN W = *pW, F = check_hyperell(W);
847 : long v;
848 155078 : if (!F)
849 7 : pari_err_TYPE(fun, W);
850 155071 : if (degpol(F) <= 0) pari_err_CONSTPOL(fun);
851 155064 : v = varn(F);
852 155064 : if (typ(W)==t_POL) W = mkvec2(W, pol_0(v));
853 : else
854 : {
855 154532 : GEN P = gel(W, 1), Q = gel(W, 2);
856 154532 : long g = hyperellgenus(F);
857 154532 : if( typ(P)!=t_POL) P = scalarpol(P, v);
858 154532 : if( typ(Q)!=t_POL) Q = scalarpol(Q, v);
859 154532 : if (degpol(P) > 2*g+2)
860 0 : pari_err_DOMAIN(fun, "poldegree(P)", ">", utoi(2*g+2), P);
861 154532 : if (degpol(Q) > g+1)
862 0 : pari_err_DOMAIN(fun, "poldegree(Q)", ">", utoi(g+1), Q);
863 :
864 154532 : W = mkvec2(P, Q);
865 : }
866 155064 : if (pF) *pF = F;
867 155064 : *pW = W;
868 155064 : }
869 :
870 : GEN
871 20629 : hyperellcharpoly(GEN PQ)
872 : {
873 20629 : pari_sp av = avma;
874 20629 : GEN M, R, T=NULL, pp=NULL, q;
875 20629 : long d, n, eps = 0;
876 : ulong p;
877 20629 : GEN H = check_hyperell(PQ);
878 20629 : if (!H || !RgX_is_FpXQX(H, &T, &pp) || !pp)
879 0 : pari_err_TYPE("hyperellcharpoly", PQ);
880 20629 : p = itou(pp);
881 20629 : if (!T)
882 : {
883 20482 : if (p==2 && is_vec_t(typ(PQ)))
884 : {
885 392 : long dP, dQ, v = varn(H);
886 392 : GEN P = gel(PQ,1), Q = gel(PQ,2);
887 392 : if (typ(P)!=t_POL) P = scalarpol(P, v);
888 392 : if (typ(Q)!=t_POL) Q = scalarpol(Q, v);
889 392 : dP = degpol(P); dQ = degpol(Q);
890 392 : if (dP<=6 && dQ <=3 && (dQ==3 || dP>=5))
891 : {
892 392 : GEN P2 = RgX_to_F2x(P), Q2 = RgX_to_F2x(Q);
893 392 : GEN D = F2x_add(F2x_mul(P2, F2x_sqr(F2x_deriv(Q2))), F2x_sqr(F2x_deriv(P2)));
894 392 : if (F2x_degree(F2x_gcd(D, Q2))) is_sing(PQ, 2);
895 392 : if (dP==6 && dQ<3 && F2x_coeff(P2,5)==F2x_coeff(Q2,2))
896 0 : is_sing(PQ, 2); /* The curve is singular at infinity */
897 392 : R = zx_to_ZX(F2x_genus2charpoly_naive(P2, Q2));
898 392 : return gc_upto(av, R);
899 : }
900 : }
901 20090 : H = RgX_to_FpX(H, pp);
902 20090 : d = degpol(H);
903 20090 : if (d <= 0) is_sing(H, p);
904 20090 : if (p > 2 && ((d == 5 && p < 17500) || (d == 6 && p < 24500)))
905 : {
906 19551 : GEN Hp = ZX_to_Flx(H, p);
907 19551 : if (!Flx_is_squarefree(Hp, p)) is_sing(H, p);
908 19544 : R = zx_to_ZX(Flx_genus2charpoly_naive(Hp, p));
909 19544 : return gc_upto(av, R);
910 : }
911 539 : n = hyperell_Weil_bound(pp, (d-1)>>1, pp);
912 539 : eps = odd(d)? 0: Fp_issquare(leading_coeff(H), pp);
913 539 : M = hyperellpadicfrobenius(H, p, n);
914 539 : R = centerlift(carberkowitz(M, 0));
915 539 : q = pp;
916 : }
917 : else
918 : {
919 : int fixvar;
920 147 : T = typ(T)==t_FFELT? FF_mod(T): RgX_to_FpX(T, pp);
921 147 : q = powuu(p, degpol(T));
922 147 : fixvar = (varncmp(varn(T),varn(H)) <= 0);
923 147 : if (fixvar) setvarn(T, fetch_var());
924 147 : H = RgX_to_FpXQX(H, T, pp);
925 147 : d = degpol(H);
926 147 : if (d <= 0) is_sing(H, p);
927 147 : eps = odd(d)? 0: Fq_issquare(leading_coeff(H), T, pp);
928 147 : n = hyperell_Weil_bound(q, (d-1)>>1, pp);
929 147 : M = nfhyperellpadicfrobenius(H, T, p, n);
930 140 : R = simplify_shallow(centerlift(liftpol_shallow(carberkowitz(M, 0))));
931 140 : if (fixvar) (void)delete_var();
932 : }
933 679 : if (!odd(d))
934 : {
935 301 : GEN b = get_basis(p, d) == 3 ? gen_1 : q;
936 301 : GEN pn = powuu(p, n);
937 301 : R = FpX_div_by_X_x(R, eps? b: negi(b), pn, NULL);
938 301 : R = FpX_center_i(R, pn, shifti(pn,-1));
939 : }
940 679 : return gc_upto(av, charpoly_funceq(R, q));
941 : }
942 :
943 : GEN
944 70 : hyperellordinate(GEN W, GEN x)
945 : {
946 70 : pari_sp av = avma;
947 70 : if (typ(W)==t_POL)
948 : {
949 : GEN d, y;
950 42 : if (typ(x)==t_INFINITY)
951 : {
952 14 : long dW = degpol(W);
953 14 : d = odd(dW) ? gen_0: gel(W,dW+2);
954 : } else
955 28 : d = poleval(W,x);
956 42 : if (gequal0(d)) { return gc_GEN(av, mkvec(d)); }
957 35 : if (!issquareall(d, &y)) retgc_const(av, cgetg(1, t_VEC));
958 14 : return gc_GEN(av, mkvec2(y, gneg(y)));
959 : }
960 : else
961 : {
962 : GEN b, c, d, rd, y, P, Q, F;
963 28 : check_hyperell_Rg("hyperellordinate", &W, &F);
964 28 : P = gel(W,1); Q = gel(W,2);
965 28 : if (typ(x)==t_INFINITY)
966 : {
967 7 : long dP = degpol(P), dQ = degpol(Q), g = hyperellgenus(F);
968 7 : c = dP < 2*g+2 ? gen_0: gel(P,dP+2);
969 7 : b = dQ < g+1 ? gen_0: gel(Q,dQ+2);
970 : } else
971 21 : { b = poleval(Q, x); c = poleval(P, x); }
972 28 : d = gadd(gsqr(b), gmul2n(c, 2));
973 28 : if (gequal0(d)) { return gc_GEN(av, mkvec(gmul2n(gneg(b),-1))); }
974 21 : if (!issquareall(d, &rd)) retgc_const(av, cgetg(1, t_VEC));
975 14 : y = gmul2n(gsub(rd, b), -1);
976 14 : return gc_GEN(av, mkvec2(y, gsub(y,rd)));
977 : }
978 : }
979 :
980 : GEN
981 122534 : hyperelldisc(GEN PQ)
982 : {
983 122534 : pari_sp av = avma;
984 122534 : GEN D, H = check_hyperell(PQ);
985 : long g;
986 122534 : if (!H || signe(H)==0) pari_err_TYPE("hyperelldisc",PQ);
987 122534 : g = hyperellgenus(H);
988 122534 : D = gmul2n(RgX_disc(H),-4*(g+1));
989 122534 : if (odd(degpol(H))) D = gmul(D, gsqr(leading_coeff(H)));
990 122534 : return gc_upto(av, D);
991 : }
992 :
993 : static long
994 136183 : get_ep(GEN W)
995 : {
996 136183 : GEN P = gel(W,1), Q = gel(W,2);
997 136183 : if (signe(Q)==0) return ZX_lval(P,2);
998 91123 : return minss(ZX_lval(P,2), ZX_lval(Q,2));
999 : }
1000 :
1001 : static GEN
1002 55355 : algo51(GEN W, GEN M)
1003 : {
1004 55355 : GEN P = gel(W,1), Q = gel(W,2);
1005 : for(;;)
1006 10654 : {
1007 66009 : long vP = ZX_lval(P,2);
1008 66009 : long vQ = signe(Q) ? ZX_lval(Q,2): vP+1;
1009 : long r;
1010 : /* 1 */
1011 66009 : if (vQ==0) break;
1012 : /* 2 */
1013 39095 : if (vP==0)
1014 : {
1015 : GEN H, H1;
1016 : /* a */
1017 32171 : RgX_even_odd(FpX_red(P,gen_2),&H, &H1);
1018 32171 : if (signe(H1)) break;
1019 : /* b */
1020 15546 : P = ZX_add(P, ZX_mul(H, ZX_sub(Q, H)));
1021 15546 : Q = ZX_sub(Q, ZX_shifti(H, 1));
1022 15546 : vP = ZX_lval(P,2);
1023 15546 : vQ = signe(Q) ? ZX_lval(Q,2): vP+1;
1024 : }
1025 : /* 2c */
1026 22470 : if (vP==1) break;
1027 : /* 2d */
1028 10654 : r = minss(2*vQ, vP)>>1;
1029 10654 : if (M) gel(M,1) = shifti(gel(M,1), r);
1030 10654 : P = ZX_shifti(P, -2*r);
1031 10654 : Q = ZX_shifti(Q, -r);
1032 : }
1033 55355 : return mkvec2(P,Q);
1034 : }
1035 :
1036 : static GEN
1037 111190 : algo52(GEN W, GEN c, long *pt_lambda)
1038 : {
1039 : long lambda;
1040 111190 : GEN P = gel(W,1), Q = gel(W,2);
1041 : for(;;)
1042 120613 : {
1043 : GEN H, H1;
1044 : /* 1 */
1045 231803 : GEN Pc = ZX_affine(P,gen_2,c), Qc = ZX_affine(Q,gen_2,c);
1046 231803 : long mP = ZX_lval(Pc,2), mQ = signe(Qc) ? ZX_lval(Qc,2): mP+1;
1047 : /* 2 */
1048 231803 : if (2*mQ <= mP) { lambda = 2*mQ; break; }
1049 : /* 3 */
1050 198713 : if (odd(mP)) { lambda = mP; break; }
1051 : /* 4 */
1052 132758 : RgX_even_odd(FpX_red(ZX_shifti(Pc, -mP),gen_2),&H, &H1);
1053 132758 : if (signe(H1)) { lambda = mP; break; }
1054 : /* 5 */
1055 120613 : P = ZX_add(P, ZX_mul(H, ZX_sub(Q, H)));
1056 120613 : Q = ZX_sub(Q, ZX_shifti(H, 1));
1057 : }
1058 111190 : *pt_lambda = lambda;
1059 111190 : return mkvec2(P,Q);
1060 : }
1061 :
1062 : static long
1063 152517 : test53(long lambda, long ep, long g)
1064 : {
1065 152517 : return (lambda <= g+1) || (odd(g) && lambda<g+3 && ep==1);
1066 : }
1067 :
1068 : static long
1069 202829 : test55(GEN W, long ep, long g)
1070 : {
1071 202829 : GEN P = gel(W,1), Q = gel(W,2);
1072 202829 : GEN Pe = FpX_red(ep ? ZX_shifti(P,-1): P, gen_2);
1073 202829 : GEN Qe = FpX_red(ep ? ZX_shifti(Q,-1): Q, gen_2);
1074 202829 : if (ep==0)
1075 : {
1076 159867 : if (signe(Qe)!=0) return ZX_val(Qe) >= (g + 3)>>1;
1077 98619 : else return ZX_val(FpX_deriv(Pe, gen_2)) >= g+1;
1078 : }
1079 : else
1080 42962 : return ZX_val(Qe) >= (g+1)>>1 && ZX_val(Pe) >= g + 1;
1081 : }
1082 :
1083 : static GEN
1084 54907 : hyperell_reverse(GEN W, long g)
1085 : {
1086 54907 : return mkvec2(RgXn_recip_shallow(gel(W,1),2*g+3),
1087 54907 : RgXn_recip_shallow(gel(W,2),g+2));
1088 : }
1089 :
1090 : /* [P,Q] -> [P(2x)/4^r, Q(2x)/2^r] */
1091 : static GEN
1092 180246 : ZX2_unscale(GEN W, long r)
1093 : {
1094 180246 : GEN P = ZX_unscale2n(gel(W,1), 1);
1095 180246 : GEN Q = ZX_unscale2n(gel(W,2), 1);
1096 180246 : if (r)
1097 : {
1098 31421 : P = ZX_shifti(P, -2*r);
1099 31421 : Q = ZX_shifti(Q, -r);
1100 : }
1101 180246 : return mkvec2(P,Q);
1102 : }
1103 : /* [P,Q] -> [P(2x+c)/4^r, Q(2x+c)/2^r] */
1104 : static GEN
1105 173325 : ZX2_affine_unscale(GEN W, long c, long r)
1106 : {
1107 253257 : if (c) W = mkvec2(ZX_Z_translate(gel(W,1), gen_1),
1108 79932 : ZX_Z_translate(gel(W,2), gen_1));
1109 173325 : return ZX2_unscale(W, r);
1110 : }
1111 :
1112 : static GEN
1113 54004 : algo56(GEN W, long g)
1114 : {
1115 : long ep;
1116 54004 : GEN M = mkvec2(gen_1, matid(2)), Woo;
1117 54004 : W = algo51(W, M);
1118 54004 : Woo = hyperell_reverse(W, g);
1119 54004 : ep = get_ep(Woo);
1120 54004 : if (test55(Woo,ep,g))
1121 : {
1122 : long lambda;
1123 12052 : Woo = algo52(Woo, gen_0, &lambda);
1124 12052 : if (!test53(lambda,ep,g))
1125 : {
1126 5969 : long r = lambda>>1;
1127 5969 : gel(M,1) = shifti(gel(M,1), r);
1128 5969 : gel(M,2) = ZM2_mul(gel(M,2), mkmat22(gen_0, gen_1, gen_2, gen_0));
1129 5969 : W = ZX2_unscale(Woo, r);
1130 : }
1131 : }
1132 : for(;;)
1133 25004 : {
1134 79008 : long j, ep = get_ep(W);
1135 199224 : for (j = 0; j < 2; j++)
1136 145220 : if (test55(ZX2_affine_unscale(W, j, 0), ep, g))
1137 : {
1138 : long lambda;
1139 96380 : GEN c = utoi(j), Wc = algo52(W, c, &lambda);
1140 96380 : if (!test53(lambda,ep,g))
1141 : {
1142 25004 : long r = lambda>>1;
1143 25004 : gel(M,1) = shifti(gel(M,1), r);
1144 25004 : gel(M,2) = ZM2_mul(gel(M,2), mkmat22(gen_2, c, gen_0, gen_1));
1145 25004 : W = ZX2_affine_unscale(Wc, j, r);
1146 25004 : break;
1147 : }
1148 : }
1149 79008 : if (j==2) break;
1150 : }
1151 54004 : return mkvec2(W, M);
1152 : }
1153 :
1154 : static GEN
1155 1351 : algo56bis(GEN W, long g, long inf, long thr)
1156 : {
1157 1351 : pari_sp av = avma;
1158 1351 : GEN vl = cgetg(3,t_VEC);
1159 1351 : long nl = 1;
1160 1351 : W = algo51(W, NULL);
1161 1351 : if (inf)
1162 : {
1163 903 : GEN Woo = hyperell_reverse(W, g);
1164 903 : long ep = get_ep(Woo);
1165 903 : if (test55(ZX2_unscale(Woo, 0), ep, g))
1166 : {
1167 : long lambda;
1168 658 : Woo = algo52(Woo, gen_0, &lambda);
1169 658 : if (lambda == thr) gel(vl,nl++) = ZX2_unscale(Woo, lambda>>1);
1170 : }
1171 : }
1172 : {
1173 1351 : long j, ep = get_ep(W);
1174 4053 : for (j = 0; j < 2; j++)
1175 2702 : if (test55(ZX2_affine_unscale(W, j, 0), ep, g))
1176 : {
1177 : long lambda;
1178 2100 : GEN Wc = algo52(W, utoi(j), &lambda);
1179 2100 : if (lambda == thr) gel(vl,nl++) = ZX2_affine_unscale(Wc, j, lambda>>1);
1180 : }
1181 : }
1182 1351 : setlg(vl, nl);
1183 1351 : return gc_GEN(av,vl);
1184 : }
1185 :
1186 : /* return the (degree 2) apolar invariant (the nth transvectant of P and P) */
1187 : static GEN
1188 1505 : ZX_apolar(GEN P, long n)
1189 : {
1190 1505 : pari_sp av = avma;
1191 1505 : long d = degpol(P), i;
1192 1505 : GEN s = gen_0, g = cgetg(n+2,t_VEC);
1193 1505 : gel(g,1) = gen_1;
1194 10563 : for (i = 1; i <= n; i++) gel(g,i+1) = muliu(gel(g,i),i); /* g[i+1] = i! */
1195 11466 : for (i = n-d; i <= d; i++)
1196 : {
1197 9961 : GEN a = mulii(mulii(gel(g,i+1),gel(g,n-i+1)),
1198 9961 : mulii(gel(P,i+2),gel(P,n-i+2)));
1199 9961 : s = odd(i)? subii(s, a): addii(s, a);
1200 : }
1201 1505 : return gc_INT(av,s);
1202 : }
1203 :
1204 : static GEN
1205 56251 : algo57(GEN F, long g, GEN pr)
1206 : {
1207 : long i, l;
1208 56251 : GEN D, C = content(F);
1209 56251 : GEN e = gel(core2(shifti(C,-vali(C))),2);
1210 56251 : GEN M = mkvec2(e, matid(2));
1211 56251 : long minvd = (2*g+1)>>(odd(g) ? 4:2);
1212 56251 : F = ZX_Z_divexact(F, sqri(e));
1213 56251 : D = absi(hyperelldisc(F));
1214 56251 : if (!pr)
1215 : {
1216 1505 : GEN A = gcdii(D, ZX_apolar(F, 2*g+2));
1217 1505 : pr = gel(factor(shifti(A, -vali(A))),1);
1218 : }
1219 56251 : l = lg(pr);
1220 318421 : for (i = 1; i < l; i++)
1221 : {
1222 : long ep;
1223 262170 : GEN p = gel(pr, i), ps2 = shifti(p,-1), Fe;
1224 262170 : if (equaliu(p,2) || Z_pval(D,p) < minvd) continue;
1225 198149 : ep = ZX_pvalrem(F,p, &Fe); Fe = FpX_red(Fe, p);
1226 198149 : if (degpol(Fe) < g+1+ep)
1227 : {
1228 6406 : GEN Fi = ZX_unscale(RgXn_recip_shallow(F,2*g+3), p);
1229 6406 : long lambda = ZX_pval(Fi,p);
1230 6406 : if (!test53(lambda,ep,g))
1231 : {
1232 3815 : GEN ppr = powiu(p,lambda>>1);
1233 3815 : F = ZX_Z_divexact(Fi,sqri(ppr));
1234 3815 : gel(M,1) = mulii(gel(M,1), ppr);
1235 3815 : gel(M,2) = ZM2_mul(gel(M,2), mkmat22(gen_0,gen_1,p,gen_0));
1236 : }
1237 : }
1238 : for(;;)
1239 25186 : {
1240 : GEN Fe, R;
1241 223335 : long j, lR, ep = ZX_pvalrem(F,p, &Fe);
1242 223335 : R = FpX_roots_mult(FpX_red(Fe, p), g+2-ep, p); lR = lg(R);
1243 235828 : for (j = 1; j<lR; j++)
1244 : {
1245 37679 : GEN c = Fp_center(gel(R,j), p, ps2);
1246 37679 : GEN Fi = ZX_affine(F,p,c);
1247 37679 : long lambda = ZX_pval(Fi,p);
1248 37679 : if (!test53(lambda,ep,g))
1249 : {
1250 25186 : GEN ppr = powiu(p,lambda>>1);
1251 25186 : F = ZX_Z_divexact(Fi, sqri(ppr));
1252 25186 : gel(M,1) = mulii(gel(M,1), ppr);
1253 25186 : gel(M,2) = ZM2_mul(gel(M,2), mkmat22(p,c,gen_0,gen_1));
1254 25186 : break;
1255 : }
1256 : }
1257 223335 : if (j==lR) break;
1258 : }
1259 : }
1260 56251 : return mkvec2(F, M);
1261 : }
1262 :
1263 : /* if inf=0, ignore point at infinity */
1264 : static GEN
1265 3563 : algo57bis(GEN F, long g, GEN p, long inf, long thr)
1266 : {
1267 3563 : pari_sp av = avma;
1268 3563 : GEN vl = cgetg(3,t_VEC), Fe;
1269 3563 : long nl = 1, ep = ZX_pvalrem(F,p, &Fe);
1270 3563 : Fe = FpX_red(Fe, p);
1271 : {
1272 3563 : GEN R = FpX_roots_mult(Fe, thr-ep, p);
1273 3563 : long j, lR = lg(R);
1274 6496 : for (j = 1; j<lR; j++)
1275 : {
1276 2933 : GEN Fj = ZX_affine(F, p, gel(R,j));
1277 2933 : long lambda = ZX_pvalrem(Fj, p, &Fj);
1278 2933 : if (lambda == thr) gel(vl,nl++) = odd(lambda)? ZX_Z_mul(Fj, p): Fj;
1279 : }
1280 : }
1281 3563 : if (inf==1 && 2*g+2-degpol(Fe) >= thr-ep)
1282 : {
1283 0 : GEN Fj = ZX_unscale(RgXn_recip_shallow(F,2*g+3), p);
1284 0 : long lambda = ZX_pvalrem(Fj, p, &Fj);
1285 0 : if (lambda == thr) gel(vl,nl++) = odd(lambda)? ZX_Z_mul(Fj, p): Fj;
1286 : }
1287 3563 : setlg(vl, nl);
1288 3563 : return gc_GEN(av,vl);
1289 : }
1290 :
1291 : static GEN
1292 4914 : next_model(GEN G, long g, GEN p, long inf, long thr)
1293 : {
1294 6265 : return equaliu(p,2) ? algo56bis(G, g, inf, thr)
1295 6265 : : algo57bis(G, g, p, inf, thr);
1296 : }
1297 :
1298 : static GEN
1299 1855 : get_extremal_even(GEN F, GEN G, long g, GEN p, long *nb)
1300 : {
1301 : while (1)
1302 1274 : {
1303 1855 : GEN Wi = next_model(G, g, p, 0, g+2);
1304 1855 : if (lg(Wi)==1) return F;
1305 1386 : F = gel(Wi,1); ++*nb;
1306 1386 : if (DEBUGLEVEL>1) err_printf("model %ld: %Ps\n", *nb, F);
1307 1386 : Wi = next_model(F, g, p, 0, g+1);
1308 1386 : if (lg(Wi)==1) return F;
1309 1274 : G = gel(Wi,1);
1310 : }
1311 : }
1312 :
1313 : static GEN
1314 0 : get_extremal_odd(GEN F, long g, GEN p, long *nb)
1315 : {
1316 : while (1)
1317 0 : {
1318 0 : GEN Wi = next_model(F, g, p, 0, g+2);
1319 0 : if (lg(Wi)==1) return F;
1320 0 : F = gel(Wi,1); ++*nb;
1321 0 : if (DEBUGLEVEL>1) err_printf("model %ld: %Ps\n", *nb, F);
1322 : }
1323 : }
1324 :
1325 : static GEN
1326 1722 : hyperellextremalmodels_nb(GEN F, long g, GEN p, long *nb)
1327 : {
1328 1722 : pari_sp av = avma;
1329 : GEN W, A, B;
1330 : long l;
1331 :
1332 1722 : *nb = 1;
1333 1722 : if (equaliu(p,2))
1334 : {
1335 917 : if (get_ep(F) > 0) retmkvec(gcopy(F));
1336 : } else
1337 : {
1338 805 : F = check_hyperell(F);
1339 805 : if (ZX_pval(F, p) > 0) return gc_GEN(av, mkvec(F));
1340 : }
1341 1673 : if (DEBUGLEVEL>1) err_printf("model %ld: %Ps\n", *nb, F);
1342 1673 : W = next_model(F, g, p, 1, odd(g)? g+2: g+1);
1343 1673 : l = lg(W); if (l==1) return gc_GEN(av, mkvec(F));
1344 518 : if (odd(g))
1345 : {
1346 0 : *nb = l-1;
1347 0 : A = get_extremal_odd(gel(W,1), g, p, nb);
1348 0 : B = l==3 ? get_extremal_odd(gel(W,2), g, p, nb) : F;
1349 : }
1350 : else
1351 : {
1352 518 : A = get_extremal_even(F, gel(W,1), g, p, nb);
1353 518 : B = l==3 ? get_extremal_even(F, gel(W,2), g, p, nb) : F;
1354 : }
1355 518 : return gc_GEN(av, A == B? mkvec(A): mkvec2(A, B));
1356 : }
1357 :
1358 : static GEN
1359 1715 : hyperellextremalmodels_i(GEN F, long g, GEN p)
1360 : {
1361 : long nb;
1362 1715 : return hyperellextremalmodels_nb(F, g, p, &nb);
1363 : }
1364 :
1365 : GEN
1366 7 : hyperellextremalmodels(GEN PQ, GEN p)
1367 : {
1368 7 : pari_sp av = avma;
1369 7 : GEN H = check_hyperell(PQ), W, v;
1370 : long g, nb;
1371 7 : if (!H || signe(H)==0) pari_err_TYPE("hyperellextremalmodels",PQ);
1372 7 : if (typ(p)!=t_INT || signe(p)<=0) pari_err_TYPE("hyperellextremalmodels",p);
1373 7 : g = hyperellgenus(H);
1374 7 : W = hyperellminimalmodel(H,NULL,mkvec(p));
1375 7 : v = cgetg(3, t_VEC);
1376 7 : gel(v, 2) = hyperellextremalmodels_nb(W, g, p, &nb);
1377 7 : gel(v, 1) = stoi(nb);
1378 7 : return gc_upto(av, v);
1379 : }
1380 :
1381 : static GEN
1382 56265 : minimalmodel_merge(GEN W2, GEN Modd, long g, long v)
1383 : {
1384 56265 : GEN P = gel(W2,1), Q = gel(W2,2);
1385 56265 : GEN e = gel(Modd,1), M = gel(Modd,2);
1386 56265 : GEN A = deg1pol_shallow(gcoeff(M,1,1), gcoeff(M,1,2), v);
1387 56265 : GEN B = deg1pol_shallow(gcoeff(M,2,1), gcoeff(M,2,2), v);
1388 56265 : GEN Bp = gpowers(B, 2*g+2);
1389 56265 : long f = mod4(e)==1 ? 1: -1;
1390 56265 : GEN m = shifti(f > 0 ? subui(1,e): addui(1,e), -2);
1391 56265 : GEN m24 = subii(shifti(m,1), shifti(sqri(m),2));
1392 56265 : P = RgX_homogenous_evalpow(P, A, Bp, 2*g+2);
1393 56265 : Q = RgX_homogenous_evalpow(Q, A, Bp, g+1);
1394 56265 : P = ZX_Z_divexact(ZX_add(P, ZX_Z_mul(ZX_sqr(Q), m24)),sqri(e));
1395 56265 : if (f < 0) Q = ZX_neg(Q);
1396 56265 : return mkvec2(P,Q);
1397 : }
1398 :
1399 : static GEN
1400 112516 : hyperell_redQ(GEN W)
1401 : {
1402 112516 : GEN P = gel(W,1), Q = gel(W,2);
1403 112516 : GEN Pr, Qr = FpX_red(Q, gen_2);
1404 112516 : Pr = ZX_add(P, ZX_shifti(ZX_mul(ZX_sub(Q, Qr),ZX_add(Q, Qr)),-2));
1405 112516 : return mkvec2(Pr, Qr);
1406 : }
1407 :
1408 : static GEN
1409 52219 : hyperellisom_finalize(GEN W1, GEN W2, GEN e, GEN M, long g, long v)
1410 : {
1411 52219 : GEN Q1 = gel(W1,2), Q2 = gel(W2,2);
1412 52219 : GEN A = deg1pol_shallow(gcoeff(M,1,1), gcoeff(M,1,2), v);
1413 52219 : GEN B = deg1pol_shallow(gcoeff(M,2,1), gcoeff(M,2,2), v);
1414 52219 : GEN Hp = RgX_homogenous_eval(Q1, A, B, g+1);
1415 52219 : GEN H = RgX_mul2n(RgX_sub(RgX_Rg_mul(Q2,e), Hp),-1);
1416 52219 : return mkvec3(e, M, H);
1417 : }
1418 :
1419 : static void
1420 56307 : check_hyperell_Q(const char *fun, GEN *pW, GEN *pF)
1421 : {
1422 56307 : GEN W = *pW, F = check_hyperell(W);
1423 : long v, g;
1424 56307 : if (!F || !signe(F) || !RgX_is_ZX(F)) pari_err_TYPE(fun, W);
1425 56300 : if (!signe(ZX_disc(F))) pari_err_DOMAIN(fun,"disc(W)","==",gen_0,W);
1426 56293 : v = varn(F); g = hyperellgenus(F);
1427 56293 : if (g == 0) pari_err_DOMAIN(fun, "genus", "=", gen_0, gen_0);
1428 56279 : if (typ(W)==t_POL) W = mkvec2(W, pol_0(v));
1429 : else
1430 : {
1431 45038 : GEN P = gel(W, 1), Q = gel(W, 2);
1432 45038 : if (typ(P)!=t_POL) P = scalarpol_shallow(P, v);
1433 45038 : if (typ(Q)!=t_POL) Q = scalarpol_shallow(Q, v);
1434 45038 : if (!RgX_is_ZX(P) || !RgX_is_ZX(Q)) pari_err_TYPE(fun,W);
1435 45038 : if (degpol(P) > 2*g+2) pari_err_DOMAIN(fun, "deg(P)", ">", utoi(2*g+2), P);
1436 45038 : if (degpol(Q) > g+1) pari_err_DOMAIN(fun, "deg(Q)", ">", utoi(g+1), Q);
1437 45038 : W = mkvec2(P, Q);
1438 : }
1439 56279 : *pW = W; *pF = F;
1440 56279 : }
1441 :
1442 : GEN
1443 56265 : hyperellminimalmodel(GEN W, GEN *pM, GEN pr)
1444 : {
1445 56265 : pari_sp av = avma;
1446 : GEN Wr, F, WM2, F2, W2, M2, Modd, Wf, ef, Mf;
1447 : long g, v;
1448 56265 : check_hyperell_Q("hyperellminimalmodel",&W, &F);
1449 56265 : if (pr && (!is_vec_t(typ(pr)) || !RgV_is_ZV(pr)))
1450 14 : pari_err_TYPE("hyperellminimalmodel",pr);
1451 56251 : g = hyperellgenus(F); v = varn(F);
1452 56251 : Wr = hyperell_redQ(W);
1453 56251 : if (!pr || RgV_isin(pr, gen_2))
1454 : {
1455 54004 : WM2 = algo56(Wr,g); W2 = gel(WM2, 1); M2 = gel(WM2, 2);
1456 54004 : F2 = check_hyperell(W2);
1457 : }
1458 : else
1459 : {
1460 2247 : W2 = Wr; F2 = F; M2 = mkvec2(gen_1, matid(2));
1461 : }
1462 56251 : Modd = gel(algo57(F2, g, pr), 2);
1463 56251 : Wf = hyperell_redQ(minimalmodel_merge(W2, Modd, g, v));
1464 56251 : if (!pM) return gc_GEN(av, Wf);
1465 50721 : ef = mulii(gel(M2,1), gel(Modd,1));
1466 50721 : Mf = ZM2_mul(gel(M2,2), gel(Modd,2));
1467 50721 : *pM = hyperellisom_finalize(W, Wf, ef, Mf, g, v);
1468 50721 : return gc_all(av, 2, &Wf, pM);
1469 : }
1470 :
1471 : GEN
1472 14 : hyperellminimaldisc(GEN W, GEN pr)
1473 : {
1474 14 : pari_sp av = avma;
1475 14 : GEN C = hyperellminimalmodel(W, NULL, pr);
1476 14 : return gc_INT(av, hyperelldisc(C));
1477 : }
1478 :
1479 : static GEN
1480 35 : redqfbsplit(GEN a, GEN b, GEN c, GEN d)
1481 : {
1482 35 : GEN p = subii(d,b), q = shifti(a,1);
1483 35 : GEN U, Q, u, v, w = bezout(p, q, &u, &v);
1484 :
1485 35 : if (!equali1(w)) { p = diviiexact(p, w); q = diviiexact(q, w); }
1486 35 : U = mkmat22(p, negi(v), q, u);
1487 35 : Q = qfb3_SL2_apply(mkvec3(a,b,c), U);
1488 35 : b = gel(Q, 2); c = gel(Q,3);
1489 35 : if (signe(b) < 0) gel(U,2) = mkcol2(v, negi(u));
1490 35 : gel(U,2) = ZC_lincomb(gen_1, truedivii(negi(c), d), gel(U,2), gel(U,1));
1491 35 : return U;
1492 : }
1493 :
1494 : static GEN
1495 16386 : polreduce(GEN P, GEN M)
1496 : {
1497 16386 : long v = varn(P), dP = degpol(P), d = odd(dP) ? dP+1: dP;
1498 16386 : GEN A = deg1pol_shallow(gcoeff(M,1,1), gcoeff(M,1,2), v);
1499 16386 : GEN B = deg1pol_shallow(gcoeff(M,2,1), gcoeff(M,2,2), v);
1500 16386 : return RgX_homogenous_evalpow(P, A, gpowers(B, d), d);
1501 : }
1502 :
1503 : /* assume deg(P) > 2 */
1504 : static GEN
1505 8193 : red_Cremona_Stoll(GEN P, GEN *pM)
1506 : {
1507 : GEN q1, q2, q3, M, R;
1508 8193 : long i, prec = nbits2prec(2*gexpo(P)) + EXTRAPRECWORD, d = degpol(P);
1509 8193 : GEN dP = ZX_deriv(P);
1510 : for (;;)
1511 0 : {
1512 8193 : GEN r = QX_complex_roots(P, prec);
1513 8193 : q1 = gen_0; q2 = gen_0; q3 = gen_0;
1514 41000 : for (i = 1; i <= d; i++)
1515 : {
1516 32807 : GEN ri = gel(r,i);
1517 32807 : GEN s = ginv(gabs(RgX_cxeval(dP,ri,NULL), prec));
1518 32807 : if (d!=4) s = gpow(s, gdivgs(gen_2,d-2), prec);
1519 32807 : q1 = gadd(q1, s);
1520 32807 : q2 = gsub(q2, gmul(real_i(ri), s));
1521 32807 : q3 = gadd(q3, gmul(gnorm(ri), s));
1522 : }
1523 8193 : M = lllgram(mkmat22(q1,q2,q2,q3));
1524 8193 : if (M && lg(M) == 3) break;
1525 0 : prec = precdbl(prec);
1526 : }
1527 8193 : R = polreduce(P, M);
1528 8193 : *pM = M;
1529 8193 : return R;
1530 : }
1531 :
1532 : /* assume deg(P) > 2 */
1533 : GEN
1534 8193 : ZX_hyperellred(GEN P, GEN *pM)
1535 : {
1536 8193 : pari_sp av = avma;
1537 8193 : long d = degpol(P);
1538 : GEN q1, q2, q3, D, vD;
1539 8193 : GEN a = gel(P,d+2), b = gel(P,d+1), c = gel(P, d);
1540 : GEN M, R, M2;
1541 :
1542 8193 : q1 = muliu(sqri(a), d);
1543 8193 : q2 = shifti(mulii(a,b), 1);
1544 8193 : q3 = subii(sqri(b), shifti(mulii(a,c), 1));
1545 8193 : D = gcdii(gcdii(q1, q2), q3);
1546 8193 : if (!equali1(D))
1547 : {
1548 8172 : q1 = diviiexact(q1, D);
1549 8172 : q2 = diviiexact(q2, D);
1550 8172 : q3 = diviiexact(q3, D);
1551 : }
1552 8193 : D = qfb_disc3(q1, q2, q3);
1553 8193 : if (!signe(D))
1554 49 : M = mkmat22(gen_1, truedivii(negi(q2),shifti(q1,1)), gen_0, gen_1);
1555 8144 : else if (issquareall(D,&vD))
1556 35 : M = redqfbsplit(q1, q2, q3, vD);
1557 : else
1558 8109 : M = gel(qfbredsl2(mkqfb(q1,q2,q3,D), NULL), 2);
1559 8193 : R = red_Cremona_Stoll(polreduce(P, M), &M2);
1560 8193 : if (pM) *pM = gmul(M, M2);
1561 8193 : return gc_all(av, pM ? 2: 1, &R, pM);
1562 : }
1563 :
1564 : GEN
1565 42 : hyperellred(GEN W, GEN *pM)
1566 : {
1567 42 : pari_sp av = avma;
1568 : long g, v;
1569 : GEN F, M, Wf;
1570 42 : check_hyperell_Q("hyperellred", &W, &F);
1571 14 : g = hyperellgenus(F); v = varn(F);
1572 14 : (void) ZX_hyperellred(F, &M);
1573 14 : Wf = hyperell_redQ(minimalmodel_merge(W, mkvec2(gen_1, M), g, v));
1574 14 : if (pM) *pM = hyperellisom_finalize(W, Wf, gen_1, M, g, v);
1575 14 : return gc_all(av, pM ? 2: 1, &Wf, pM);
1576 : }
1577 :
1578 : static void
1579 154511 : check_hyperell_vc(const char *fun, GEN C, long v, GEN *e, GEN *M, GEN *H)
1580 : {
1581 154511 : if (typ(C) != t_VEC || lg(C) != 4) pari_err_TYPE(fun,C);
1582 154504 : *e = gel(C,1); *M = gel(C,2); *H = gel(C,3);
1583 154504 : if (typ(*M) != t_MAT || lg(*M) != 3 || lgcols(*M) != 3) pari_err_TYPE(fun,C);
1584 154497 : if (typ(*H) != t_POL || varncmp(varn(*H),v) > 0) *H = scalarpol_shallow(*H,v);
1585 154497 : if (varncmp(gvar(*M),v) <= 0) pari_err_PRIORITY(fun,*M,"<=",v);
1586 154497 : }
1587 :
1588 : GEN
1589 110509 : hyperellchangecurve(GEN W, GEN C)
1590 : {
1591 110509 : pari_sp av = avma;
1592 : GEN F, P, Q, A, B, Bp, e, M, H;
1593 : long g, v;
1594 :
1595 110509 : check_hyperell_Rg("hyperellchangecurve",&W,&F);
1596 110495 : P = gel(W,1); Q = gel(W,2);
1597 110495 : g = hyperellgenus(F); v = varn(F);
1598 110495 : check_hyperell_vc("hyperellchangecurve", C, v, &e, &M, &H);
1599 110481 : A = deg1pol_shallow(gcoeff(M,1,1), gcoeff(M,1,2), v);
1600 110481 : B = deg1pol_shallow(gcoeff(M,2,1), gcoeff(M,2,2), v);
1601 110481 : Bp = gpowers(B, 2*g+2);
1602 110481 : P = RgX_homogenous_evalpow(P, A, Bp, 2*g+2);
1603 110481 : Q = RgX_homogenous_evalpow(Q, A, Bp, g+1);
1604 110481 : P = RgX_Rg_div(RgX_sub(P, RgX_mul(H,RgX_add(Q,H))), gsqr(e));
1605 110481 : Q = RgX_Rg_div(RgX_add(Q, RgX_mul2n(H,1)), e);
1606 110481 : return gc_GEN(av, mkvec2(P,Q));
1607 : }
1608 :
1609 : static int
1610 413 : checkhyperellpt_i(GEN pt, GEN *x, GEN *y, GEN *z)
1611 : {
1612 413 : if (typ(pt) != t_VEC || lg(pt)<2 || lg(pt)>4)
1613 0 : { *x=NULL; *y=NULL; *z=NULL; return 0; }
1614 413 : if (lg(pt) == 2)
1615 : {
1616 0 : *x = gen_1; *y = gel(pt,1); *z = gen_0;
1617 : } else
1618 : {
1619 413 : *x = gel(pt,1); *y = gel(pt,2);
1620 413 : *z = lg(pt)==3 ? gen_1: gel(pt, 3);
1621 : }
1622 413 : return 1;
1623 : }
1624 :
1625 : static GEN
1626 329 : wprojtoaff(GEN X, GEN Y, GEN Z, GEN pt, long g)
1627 : {
1628 329 : if (lg(pt)==4) return mkvec3(X,Y,Z);
1629 329 : return gequal0(Z) ? mkvec(gequal0(Y) ? gen_0: gdiv(Y,gpowgs(X,g+1)))
1630 329 : : mkvec2(gdiv(X,Z),gequal0(Y) ? gen_0: gdiv(Y, gpowgs(Z,g+1)));
1631 : }
1632 :
1633 : GEN
1634 70 : hyperellchangepointinv(GEN W, GEN pt, GEN C)
1635 : {
1636 70 : pari_sp av = avma;
1637 : GEN F, e, M, H, x, y, z, X, Y,Z;
1638 : long g, v;
1639 :
1640 70 : check_hyperell_Rg("hyperellchangepointinv",&W,&F);
1641 70 : g = hyperellgenus(F); v = varn(F);
1642 70 : check_hyperell_vc("hyperellchangepointinv", C, v, &e, &M, &H);
1643 70 : if (!checkhyperellpt_i(pt,&x,&y,&z))
1644 0 : pari_err_TYPE("hyperellchangepointinv",pt);
1645 70 : X = gadd(gmul(gcoeff(M,1,1), x), gmul(gcoeff(M,1,2),z));
1646 70 : Z = gadd(gmul(gcoeff(M,2,1), x), gmul(gcoeff(M,2,2),z));
1647 70 : Y = gadd(gmul(e, y), RgX_homogenous_eval(H, x, z, g+1));
1648 70 : return gc_GEN(av, wprojtoaff(X,Y,Z,pt,g));
1649 : }
1650 :
1651 : GEN
1652 259 : hyperellchangepoint(GEN W, GEN pt, GEN C)
1653 : {
1654 259 : pari_sp av = avma;
1655 : GEN F, e, M, H, x, y, z, X, Y, Z;
1656 : GEN a, b, c, d, D;
1657 : long g, v;
1658 :
1659 259 : check_hyperell_Rg("hyperellchangepoint",&W,&F);
1660 259 : g = hyperellgenus(F); v = varn(F);
1661 259 : check_hyperell_vc("hyperellchangepoint", C, v, &e, &M, &H);
1662 259 : if (!checkhyperellpt_i(pt,&x,&y,&z))
1663 0 : pari_err_TYPE("hyperellchangepoint",pt);
1664 259 : a = gcoeff(M,1,1); b = gcoeff(M,1,2);
1665 259 : c = gcoeff(M,2,1); d = gcoeff(M,2,2);
1666 259 : Z = gsub(gmul(a, z), gmul(c, x));
1667 259 : X = gsub(gmul(d, x), gmul(b, z));
1668 259 : D = gsub(gmul(a, d), gmul(b, c));
1669 259 : Y = gdiv(gsub(gmul(y, gpowgs(D, g+1)), RgX_homogenous_eval(H,X,Z,g+1)), e);
1670 259 : return gc_GEN(av, wprojtoaff(X,Y,Z,pt,g));
1671 : }
1672 :
1673 : GEN
1674 43561 : hyperellchangeinvert(GEN W, GEN C)
1675 : {
1676 43561 : pari_sp av = avma;
1677 : GEN F, e, M, H, ei, Mi, Hi, X, Z, Zp;
1678 : long g, v;
1679 43561 : check_hyperell_Rg("hyperellchangeinvert",&W,&F);
1680 43561 : g = hyperellgenus(F); v = varn(F);
1681 43561 : check_hyperell_vc("hyperellchangeinvert", C, v, &e, &M, &H);
1682 43561 : ei = ginv(e);
1683 43561 : Mi = RgM_inv(M);
1684 43561 : X = deg1pol_shallow(gcoeff(Mi,1,1), gcoeff(Mi,1,2), v);
1685 43561 : Z = deg1pol_shallow(gcoeff(Mi,2,1), gcoeff(Mi,2,2), v);
1686 43561 : Zp = gpowers(Z, g+1);
1687 43561 : Hi = gmul(ei, gneg(RgX_homogenous_evalpow(H, X, Zp, g+1)));
1688 43561 : return gc_GEN(av, mkvec3(ei, Mi, Hi));
1689 : }
1690 :
1691 : GEN
1692 63 : hyperellchangecompose(GEN W, GEN C1, GEN C2)
1693 : {
1694 63 : pari_sp av = avma;
1695 : GEN F, e1, M1, H1, e2, M2, H2, H, X, Z, Zp;
1696 : long g, v;
1697 63 : check_hyperell_Rg("hyperellchangecompose",&W,&F);
1698 63 : g = hyperellgenus(F); v = varn(F);
1699 63 : check_hyperell_vc("hyperellchangecompose", C1, v, &e1, &M1, &H1);
1700 63 : check_hyperell_vc("hyperellchangecompose", C2, v, &e2, &M2, &H2);
1701 63 : X = deg1pol_shallow(gcoeff(M2,1,1), gcoeff(M2,1,2), v);
1702 63 : Z = deg1pol_shallow(gcoeff(M2,2,1), gcoeff(M2,2,2), v);
1703 63 : Zp = gpowers(Z, g+1);
1704 63 : H = gadd(gmul(e1,H2),RgX_homogenous_evalpow(H1, X, Zp, g+1));
1705 63 : return gc_GEN(av, mkvec3(gmul(e1,e2), gmul(M1, M2), H));
1706 : }
1707 :
1708 : int
1709 84 : hyperellisoncurve(GEN W, GEN P)
1710 : {
1711 84 : pari_sp av = avma;
1712 : GEN x, y, z, F;
1713 : long g;
1714 : int res;
1715 84 : check_hyperell_Rg("hyperellisoncurve",&W,&F);
1716 84 : g = hyperellgenus(F);
1717 84 : if (!checkhyperellpt_i(P,&x,&y,&z)) pari_err_TYPE("hyperellisoncurve",P);
1718 84 : if (typ(W)==t_POL)
1719 0 : res = gequal(gsqr(y), RgX_homogenous_eval(W,x,z,2*g+2));
1720 : else
1721 : {
1722 : GEN zp;
1723 84 : if (typ(W)!=t_VEC || lg(W)!=3) pari_err_TYPE("hyperellisoncurve",W);
1724 84 : zp = gpowers(z, 2*g+2);
1725 84 : res = gequal(gmul(y, gadd(y,RgX_homogenous_evalpow(gel(W,2), x,zp,g+1))),
1726 84 : RgX_homogenous_evalpow(gel(W,1),x,zp,2*(g+1)));
1727 : }
1728 84 : return gc_int(av, res);
1729 : }
1730 :
1731 : /****************************************************************************/
1732 : /*** ***/
1733 : /*** genus2charpoly ***/
1734 : /*** ***/
1735 : /****************************************************************************/
1736 :
1737 : /* Half stable reduction */
1738 :
1739 : static long
1740 588 : Zst_val(GEN P, GEN f, GEN p, long vt, GEN *pR)
1741 : {
1742 588 : pari_sp av = avma;
1743 588 : long v = varn(P);
1744 : while(1)
1745 1260 : {
1746 1848 : long i, j, dm = LONG_MAX;
1747 1848 : GEN Pm = NULL;
1748 1848 : long dP = degpol(P);
1749 7532 : for (i = 0; i <= minss(dP, dm); i++)
1750 : {
1751 5684 : GEN Py = gel(P, i+2);
1752 5684 : if (signe(Py))
1753 : {
1754 4186 : if (typ(Py)==t_POL)
1755 : {
1756 3864 : long dPy = degpol(Py);
1757 12502 : for (j = 0; j <= minss(dPy, dm-i); j++)
1758 : {
1759 8638 : GEN c = gel(Py, j+2);
1760 8638 : if (signe(c))
1761 : {
1762 3556 : if (i+j < dm)
1763 : {
1764 1848 : dm = i+j;
1765 1848 : Pm = monomial(gen_1, dm, v);
1766 1848 : gel(Pm,dm+2) = gen_0;
1767 : }
1768 3556 : gel(Pm,i+2) = c;
1769 : }
1770 : }
1771 : } else
1772 : {
1773 322 : if (i < dm)
1774 : {
1775 77 : dm = i;
1776 77 : Pm = monomial(Py, dm, v);
1777 : }
1778 : else
1779 245 : gel(Pm, i+2) = Py;
1780 : }
1781 : }
1782 : }
1783 1848 : Pm = RgX_renormalize(Pm);
1784 1848 : if (ZX_pval(Pm,p)==0)
1785 : {
1786 588 : *pR = gc_GEN(av, P);
1787 588 : return dm;
1788 : }
1789 1260 : Pm = RgX_homogenize_deg(Pm, dm, vt);
1790 1260 : P = gadd(gsub(P, Pm), gmul(f, ZXX_Z_divexact(Pm, p)));
1791 : }
1792 : }
1793 :
1794 : static long
1795 588 : Zst_normval(GEN P, GEN f, GEN p, long vt, GEN *pR)
1796 : {
1797 588 : long v = Zst_val(P, f, p, vt, pR);
1798 588 : long e = RgX_val(*pR)>>1;
1799 588 : if (e > 0)
1800 : {
1801 0 : v -= 2*e;
1802 0 : *pR = RgX_shift(*pR, -2*e);
1803 : }
1804 588 : return v;
1805 : }
1806 :
1807 : static GEN
1808 1176 : RgXY_swapsafe(GEN P, long v1, long v2)
1809 : {
1810 1176 : if (varn(P)==v2)
1811 : {
1812 77 : P = shallowcopy(P); setvarn(P,v1); return P;
1813 : } else
1814 1099 : return RgXY_swap(P, RgXY_degreex(P), v2);
1815 : }
1816 :
1817 : static GEN
1818 588 : Zst_red1(GEN P, GEN f, GEN p, long vt)
1819 : {
1820 588 : pari_sp av = avma;
1821 : GEN r, f1, f2, P1, P2;
1822 588 : long vs = varn(P);
1823 588 : long w = Zst_normval(P, f, p, vt, &r), ww = w-odd(w);
1824 588 : GEN st = monomial(pol_x(vt), 1, vs);
1825 588 : f1 = gsubst(f, vt, st);
1826 588 : P1 = gsubst(gdiv(r, monomial(gen_1,ww,vs)),vt,st);
1827 588 : f2 = gsubst(f, vs, st);
1828 588 : P2 = gsubst(gdiv(r, monomial(gen_1,ww,vt)),vs,st);
1829 588 : f2 = RgXY_swapsafe(f2, vs, vt);
1830 588 : P2 = RgXY_swapsafe(P2, vs, vt);
1831 588 : return gc_GEN(av, mkvec4(P1, f1, P2, f2));
1832 : }
1833 :
1834 : static GEN
1835 1176 : Zst_reduce(GEN P, GEN p, long vt, long *pv)
1836 : {
1837 : GEN C;
1838 1176 : long v = RgX_val(P);
1839 1176 : *pv = v + ZXX_pvalrem(RgX_shift(P, -v), p, &P);
1840 1176 : C = constant_coeff(P);
1841 1176 : C = typ(C) == t_POL ? C: scalarpol_shallow(C, vt);
1842 1176 : return FpX_red(C, p);
1843 : }
1844 :
1845 : static GEN
1846 588 : Zst_red3(GEN C, GEN p, long vt)
1847 : {
1848 : while(1)
1849 511 : {
1850 588 : GEN P1 = gel(C,1), f1 = gel(C,2), Poo = gel(C,3), foo= gel(C,4);
1851 : long e;
1852 588 : GEN Qoop = Zst_reduce(Poo, p, vt, &e), Qp, R;
1853 588 : if (RgX_val(Qoop) >= 3-e)
1854 : {
1855 0 : C = Zst_red1(Poo, foo, p, vt);
1856 511 : continue;
1857 : }
1858 588 : Qp = Zst_reduce(P1, p, vt, &e);
1859 588 : R = FpX_roots_mult(Qp, 3-e, p);
1860 588 : if (lg(R) > 1)
1861 511 : {
1862 511 : GEN xz = deg1pol_shallow(gen_1, gel(R,1), vt);
1863 511 : C = Zst_red1(gsubst(P1, vt, xz), gsubst(f1, vt, xz), p, vt);
1864 511 : continue;
1865 : }
1866 77 : return Qp;
1867 : }
1868 : }
1869 :
1870 : static GEN
1871 77 : genus2_halfstablemodel_i(GEN P, GEN p, long vt)
1872 : {
1873 : GEN Qp, R, Poo, Qoop;
1874 77 : long e = ZX_pvalrem(P, p, &Qp);
1875 77 : R = FpX_roots_mult(FpX_red(Qp,p), 4-e, p);
1876 77 : if (lg(R) > 1)
1877 : {
1878 77 : GEN C = Zst_red1(ZX_Z_translate(P, gel(R,1)), pol_x(vt), p, vt);
1879 77 : return Zst_red3(C, p, vt);
1880 : }
1881 0 : Poo = RgXn_recip_shallow(P, 7);
1882 0 : e = ZX_pvalrem(Poo, p, &Qoop);
1883 0 : Qoop = FpX_red(Qoop,p);
1884 0 : if (RgX_val(Qoop)>=4-e)
1885 : {
1886 0 : GEN C = Zst_red1(Poo, pol_x(vt), p, vt);
1887 0 : return Zst_red3(C, p, vt);
1888 : }
1889 0 : return gcopy(P);
1890 : }
1891 :
1892 : static GEN
1893 77 : genus2_halfstablemodel(GEN P, GEN p)
1894 : {
1895 77 : pari_sp av = avma;
1896 77 : long vt = fetch_var(), vs = varn(P);
1897 77 : GEN S = genus2_halfstablemodel_i(P, p, vt);
1898 77 : setvarn(S, vs); delete_var();
1899 77 : return gc_GEN(av, S);
1900 : }
1901 :
1902 : /* semi-stable reduction */
1903 :
1904 : static GEN
1905 1015 : genus2_redmodel(GEN P, GEN p)
1906 : {
1907 : GEN LP, U, F;
1908 : long i, k, r;
1909 1015 : if (degpol(P) < 0) return mkvec2(cgetg(1, t_COL), P);
1910 980 : F = FpX_factor_squarefree(P, p);
1911 980 : r = lg(F); U = NULL;
1912 3416 : for (i = k = 1; i < r; i++)
1913 : {
1914 2436 : GEN f = gel(F,i);
1915 2436 : long df = degpol(f);
1916 2436 : if (!df) continue;
1917 1687 : if (odd(i)) U = U? FpX_mul(U, f, p): f;
1918 1687 : if (i > 1) gel(F,k++) = df == 1? mkcol(f): gel(FpX_factor(f, p), 1);
1919 : }
1920 980 : LP = leading_coeff(P);
1921 980 : if (!U)
1922 154 : U = scalarpol_shallow(LP, varn(P));
1923 : else
1924 : {
1925 826 : GEN LU = leading_coeff(U);
1926 826 : if (!equalii(LU, LP)) U = FpX_Fp_mul(U, Fp_div(LP, LU, p), p);
1927 : }
1928 980 : setlg(F,k); if (k > 1) F = shallowconcat1(F);
1929 980 : return mkvec2(F, U);
1930 : }
1931 :
1932 : static GEN
1933 8834 : xdminusone(long d)
1934 : {
1935 8834 : return gsub(pol_xn(d, 0),gen_1);
1936 : }
1937 :
1938 : static GEN
1939 637 : ellfromeqncharpoly(GEN P, GEN Q, GEN p)
1940 : {
1941 : long v;
1942 : GEN E, F, t, y;
1943 637 : v = fetch_var();
1944 637 : y = pol_x(v);
1945 637 : F = gsub(gadd(ZX_sqr(y), gmul(y, Q)), P);
1946 637 : E = ellinit(ellfromeqn(F), p, DEFAULTPREC);
1947 637 : delete_var();
1948 637 : t = ellcharpoly(E, p);
1949 637 : obj_free(E);
1950 637 : return t;
1951 : }
1952 :
1953 : static GEN
1954 1722 : RgX_remswap(GEN P, GEN f, long vy)
1955 : {
1956 1722 : GEN R = RgX_rem(RgXY_swap(P, 3, vy), gsub(f, pol_x(vy)));
1957 1722 : return RgXY_swap(R, 3, vy);
1958 : }
1959 :
1960 : static GEN
1961 721 : ftrans(GEN f, GEN r, GEN p)
1962 : {
1963 721 : r = shallowcopy(r);
1964 721 : setvarn(r, varn(f));
1965 721 : return RgX_Rg_div(gsub(f, r), p);
1966 : }
1967 :
1968 : static GEN
1969 651 : algo52_F4(GEN W, GEN T, GEN c, GEN f, long *pt_lambda)
1970 : {
1971 651 : GEN P = gel(W,1), Q = gel(W,2);
1972 651 : long lambda, vy = varn(T);
1973 651 : GEN fc = ftrans(f,c,gen_2);
1974 : for(;;)
1975 105 : {
1976 : GEN H, H1;
1977 : /* 1 */
1978 756 : GEN Pc = RgX_remswap(RgX_affine(P,gen_2,c), fc, vy);
1979 756 : GEN Qc = RgX_remswap(RgX_affine(Q,gen_2,c), fc, vy);
1980 756 : long mP = ZXX_pval(Pc,gen_2), mQ = signe(Qc) ? ZXX_pval(Qc,gen_2): mP+1;
1981 : /* 2 */
1982 756 : if (2*mQ <= mP) { lambda = 2*mQ; break; }
1983 : /* 3 */
1984 693 : if (odd(mP)) { lambda = mP; break; }
1985 : /* 4 */
1986 280 : RgX_even_odd(FpXX_red(ZXX_shifti(Pc, -mP),gen_2),&H, &H1);
1987 280 : if (signe(H1)) { lambda = mP; break; }
1988 : /* 5 */
1989 105 : H = RgX_deflate(FpXQX_sqr(H, T, gen_2), 2);
1990 105 : P = RgX_add(P, RgX_mul(H, RgX_sub(Q, H)));
1991 105 : P = RgX_remswap(P, f, vy);
1992 105 : Q = RgX_sub(Q, RgX_mul2n(H, 1));
1993 : }
1994 651 : *pt_lambda = lambda;
1995 651 : return mkvec2(P,Q);
1996 : }
1997 :
1998 : static GEN
1999 651 : genus2_tr2(GEN W, GEN r, GEN *f, GEN T)
2000 : {
2001 651 : long lambda, v, vy = varn(T);
2002 : GEN P, Q;
2003 651 : W = algo52_F4(W, T, r, *f, &lambda);
2004 651 : if (lambda < 2) return NULL;
2005 273 : v = lambda>>1;
2006 273 : if (signe(r))
2007 : {
2008 35 : *f = ftrans(*f, r, gen_2);
2009 35 : P = RgX_remswap(RgX_affine(gel(W,1), gen_2, r), *f, vy);
2010 35 : Q = RgX_remswap(RgX_affine(gel(W,2), gen_2, r), *f, vy);
2011 : } else
2012 : {
2013 238 : *f = RgX_mul2n(*f, -1);
2014 238 : P = RgX_unscale(gel(W,1), gen_2);
2015 238 : Q = RgX_unscale(gel(W,2), gen_2);
2016 : }
2017 273 : return mkvec2(ZXX_shifti(P, -2*v), ZXX_shifti(Q, -v));
2018 : }
2019 :
2020 : static GEN
2021 273 : genus2_red2(GEN W, GEN T, GEN f, GEN p)
2022 : {
2023 : while(1)
2024 56 : {
2025 : long i, l;
2026 273 : GEN P = gel(W,1), Q = gel(W,2), Pr, R;
2027 273 : (void) ZXX_pvalrem(P, p, &Pr);
2028 273 : R = FpXQX_roots(FpXQX_gcd(Pr,Q,T,p), T, p);
2029 273 : l = lg(R);
2030 273 : if (l < 2) break;
2031 581 : for (i = 1; i < l; i++)
2032 : {
2033 399 : GEN W2 = genus2_tr2(W, gel(R,i), &f, T);
2034 399 : if (!W2) continue;
2035 56 : W = W2;
2036 56 : break;
2037 : }
2038 238 : if (i == l) break;
2039 : }
2040 217 : return W;
2041 : }
2042 :
2043 : static GEN
2044 651 : cf(GEN P, long i, long v)
2045 651 : { return i <= degpol(P) ? to_ZX(gel(P,i+2), v) : pol_0(v); }
2046 :
2047 : static GEN
2048 651 : cfu(GEN P, long i, GEN u, long v)
2049 651 : { return i <= degpol(P) ? ZX_mul(to_ZX(gel(P,i+2), v), u) : pol_0(v); }
2050 :
2051 : static GEN
2052 1127 : genus2_type5ns_2(GEN P, GEN Q, GEN p)
2053 : {
2054 1127 : pari_sp av = avma;
2055 : GEN FP, FQ, F, T, P1, Q1, E, W;
2056 : GEN a1, a2, a3, a4, a6, u, f;
2057 1127 : long v, vy = varn(P);
2058 1127 : (void) ZXX_pvalrem(P, p, &FP);
2059 1127 : if (signe(Q)) (void) ZXX_pvalrem(Q, p, &FQ);
2060 1127 : FP = FpX_red(FP, p);
2061 1127 : FQ = signe(Q) ? FpX_red(FQ, p): Q;
2062 1127 : F = FpX_gcd(FP, FpX_sqr(FQ, p), p);
2063 1127 : T = deg2pol_shallow(gen_1, gen_1, gen_1, vy);
2064 1127 : if (signe(FpX_rem(F,FpX_sqr(T, p), p))) return NULL;
2065 252 : v = fetch_var_higher();
2066 252 : P1 = RgV_to_RgX(ZX_digits(P, T), v);
2067 252 : Q1 = RgV_to_RgX(ZX_digits(Q, T), v);
2068 252 : f = shallowcopy(T); setvarn(f, v);
2069 252 : W = genus2_tr2(mkvec2(P1,Q1), pol_0(vy), &f, T);
2070 252 : if (!W) { delete_var(); return NULL; }
2071 217 : E = genus2_red2(W, T, f, p);
2072 217 : P = FpXX_red(gel(E,1), gen_2); Q = FpXX_red(gel(E,2), gen_2);
2073 217 : u = cf(P, 3, vy);
2074 217 : a1 = cf(Q, 1, vy);
2075 217 : a2 = cf(P, 2, vy);
2076 217 : a3 = cfu(Q, 0, u, vy);
2077 217 : a4 = cfu(P, 1, u, vy);
2078 217 : a6 = ZX_mul(cfu(P, 0, u, vy), u);
2079 217 : delete_var();
2080 217 : E = mkvec5(a1, a2, a3, a4, a6);
2081 217 : return gc_GEN(av, RgX_inflate(FpXQV_ellcharpoly(E, T, p), 2));
2082 : }
2083 :
2084 : static GEN
2085 14 : genus2_red5(GEN P, GEN T, GEN p)
2086 : {
2087 14 : long vx = varn(P), vy = varn(T);
2088 14 : GEN f = shallowcopy(T), pi = shifti(p,-1);
2089 14 : setvarn(f, vx);
2090 : while(1)
2091 21 : {
2092 : GEN Pr, R, r;
2093 35 : long v = ZXX_pvalrem(P, p, &Pr);
2094 35 : R = FpXQX_roots_mult(Pr, 2-v, T, p);
2095 49 : if (lg(R)==1) return P;
2096 35 : r = FpX_center(gel(R,1), p, pi);
2097 35 : Pr = RgX_affine(P, p, r);
2098 35 : f = ftrans(f, r, p);
2099 35 : Pr = RgX_remswap(Pr, f, vy);
2100 35 : if (ZXX_pvalrem(Pr, sqri(p), &Pr)==0) return P;
2101 21 : P = Pr;
2102 : }
2103 : }
2104 :
2105 : static GEN
2106 819 : genus2_type5ns(GEN P, GEN p)
2107 : {
2108 819 : pari_sp av = avma;
2109 : GEN E, F, T, Q, u, a2, a4, a6;
2110 819 : long v, vy = varn(P);
2111 819 : if (equaliu(p, 2))
2112 0 : (void) ZXX_pvalrem(P, sqri(p), &P);
2113 819 : (void) ZX_pvalrem(P, p, &F);
2114 819 : F = FpX_red(F, p);
2115 819 : if (degpol(F) < 1) return NULL;
2116 819 : F = FpX_factor(F, p);
2117 819 : if (mael(F,2,1) != 3 || degpol(gmael(F,1,1)) != 2) return NULL;
2118 14 : T = gmael(F, 1, 1);
2119 14 : v = fetch_var_higher();
2120 14 : Q = RgV_to_RgX(ZX_digits(P, T), v);
2121 14 : Q = genus2_red5(Q, T, p);
2122 14 : u = to_ZX(gel(Q,5), vy);
2123 14 : a2 = to_ZX(gel(Q,4), vy);
2124 14 : a4 = ZX_mul(to_ZX(gel(Q,3),vy), u);
2125 14 : a6 = ZX_mul(to_ZX(gel(Q,2),vy), ZX_sqr(u));
2126 14 : E = mkvec5(pol_0(vy), a2, pol_0(vy), a4, a6);
2127 14 : delete_var();
2128 14 : return gc_GEN(av, RgX_inflate(FpXQV_ellcharpoly(E, T, p), 2));
2129 : }
2130 :
2131 : /* Assume P has semistable reduction at p */
2132 : static GEN
2133 1015 : genus2_eulerfact_semistable(GEN P, GEN p)
2134 : {
2135 1015 : GEN Pp = FpX_red(P, p);
2136 1015 : GEN GU = genus2_redmodel(Pp, p);
2137 1015 : long d = 6-degpol(Pp), v = d/2, w = odd(d);
2138 : GEN abe, tor;
2139 1015 : GEN ki, kp = pol_1(0), kq = pol_1(0);
2140 1015 : GEN F = gel(GU,1), Q = gel(GU,2);
2141 1015 : long dQ = degpol(Q), lF = lg(F)-1;
2142 :
2143 7 : abe = dQ >= 5 ? hyperellcharpoly(gmul(Q,gmodulo(gen_1,p)))
2144 2023 : : dQ >= 3 ? ellfromeqncharpoly(Q,gen_0,p)
2145 1008 : : pol_1(0);
2146 861 : ki = dQ != 0 ? xdminusone(1)
2147 1169 : : Fp_issquare(gel(Q,2),p) ? ZX_sqr(xdminusone(1))
2148 154 : : xdminusone(2);
2149 1015 : if (lF)
2150 : {
2151 : long i;
2152 2100 : for(i=1; i <= lF; i++)
2153 : {
2154 1183 : GEN Fi = gel(F, i);
2155 1183 : long d = degpol(Fi);
2156 1183 : GEN e = FpX_rem(Q, Fi, p);
2157 2100 : GEN kqf = lgpol(e)==0 ? xdminusone(d):
2158 1561 : FpXQ_issquare(e, Fi, p) ? ZX_sqr(xdminusone(d))
2159 917 : : xdminusone(2*d);
2160 1183 : kp = gmul(kp, xdminusone(d));
2161 1183 : kq = gmul(kq, kqf);
2162 : }
2163 : }
2164 1015 : if (v)
2165 : {
2166 273 : GEN kqoo = w==1 ? xdminusone(1):
2167 21 : Fp_issquare(leading_coeff(Q), p)? ZX_sqr(xdminusone(1))
2168 14 : : xdminusone(2);
2169 259 : kp = gmul(kp, xdminusone(1));
2170 259 : kq = gmul(kq, kqoo);
2171 : }
2172 1015 : tor = RgX_div(ZX_mul(xdminusone(1), kq), ZX_mul(ki, kp));
2173 1015 : return ZX_mul(abe, tor);
2174 : }
2175 :
2176 : GEN
2177 1442 : genus2_eulerfact(GEN P, GEN p, long ra, long rt)
2178 : {
2179 1442 : pari_sp av = avma;
2180 : GEN W, R, E;
2181 1442 : long d = 2*ra+rt;
2182 1442 : if (d == 0) return pol_1(0);
2183 819 : R = genus2_type5ns(P, p);
2184 819 : if (R) return R;
2185 805 : W = hyperellextremalmodels_i(P, 2, p);
2186 805 : if (lg(W) < 3)
2187 : {
2188 672 : GEN F = genus2_eulerfact_semistable(P,p);
2189 672 : if (degpol(F)!=d)
2190 : {
2191 77 : GEN S = genus2_halfstablemodel(P, p);
2192 77 : F = genus2_eulerfact_semistable(S, p);
2193 77 : if (degpol(F)!=d) pari_err_BUG("genus2charpoly");
2194 : }
2195 672 : return F;
2196 : }
2197 133 : E = gmul(genus2_eulerfact_semistable(gel(W,1),p),
2198 133 : genus2_eulerfact_semistable(gel(W,2),p));
2199 133 : return gc_upto(av, E);
2200 : }
2201 :
2202 : /* p = 2 */
2203 :
2204 : static GEN
2205 1106 : F2x_genus2_find_trans(GEN P, GEN Q, GEN F)
2206 : {
2207 1106 : pari_sp av = avma;
2208 1106 : long i, d = F2x_degree(F), v = P[1];
2209 : GEN M, C, V;
2210 1106 : M = cgetg(d+1, t_MAT);
2211 3416 : for (i=1; i<=d; i++)
2212 : {
2213 2310 : GEN Mi = F2x_rem(F2x_add(F2x_shift(Q,i-1), monomial_F2x(2*i-2,v)), F);
2214 2310 : gel(M,i) = F2x_to_F2v(Mi, d);
2215 : }
2216 1106 : C = F2x_to_F2v(F2x_rem(P, F), d);
2217 1106 : V = F2m_F2c_invimage(M, C);
2218 1106 : return gc_leaf(av, F2v_to_F2x(V, v));
2219 : }
2220 :
2221 : static GEN
2222 1547 : F2x_genus2_trans(GEN P, GEN Q, GEN H)
2223 : {
2224 1547 : return F2x_add(P,F2x_add(F2x_mul(H,Q), F2x_sqr(H)));
2225 : }
2226 :
2227 : static GEN
2228 2814 : F2x_genus_redoo(GEN P, GEN Q, long k)
2229 : {
2230 2814 : if (F2x_degree(P)==2*k)
2231 : {
2232 700 : long c = F2x_coeff(P,2*k-1), dQ = F2x_degree(Q);
2233 700 : if ((dQ==k-1 && c==1) || (dQ<k-1 && c==0))
2234 441 : return F2x_genus2_trans(P, Q, monomial_F2x(k, P[1]));
2235 : }
2236 2373 : return P;
2237 : }
2238 :
2239 : static GEN
2240 1834 : F2x_pseudodisc(GEN P, GEN Q)
2241 : {
2242 1834 : GEN dP = F2x_deriv(P), dQ = F2x_deriv(Q);
2243 1834 : return F2x_gcd(Q, F2x_add(F2x_mul(P, F2x_sqr(dQ)), F2x_sqr(dP)));
2244 : }
2245 :
2246 : static GEN
2247 938 : F2x_genus_red(GEN P, GEN Q)
2248 : {
2249 : long dP, dQ;
2250 : GEN F, FF;
2251 938 : P = F2x_genus_redoo(P, Q, 3);
2252 938 : P = F2x_genus_redoo(P, Q, 2);
2253 938 : P = F2x_genus_redoo(P, Q, 1);
2254 938 : dP = F2x_degree(P);
2255 938 : dQ = F2x_degree(Q);
2256 938 : FF = F = F2x_pseudodisc(P,Q);
2257 1834 : while(F2x_degree(F)>0)
2258 : {
2259 896 : GEN M = gel(F2x_factor(F),1);
2260 896 : long i, l = lg(M);
2261 2002 : for(i=1; i<l; i++)
2262 : {
2263 1106 : GEN R = F2x_sqr(gel(M,i));
2264 1106 : GEN H = F2x_genus2_find_trans(P, Q, R);
2265 1106 : P = F2x_div(F2x_genus2_trans(P, Q, H), R);
2266 1106 : Q = F2x_div(Q, gel(M,i));
2267 : }
2268 896 : F = F2x_pseudodisc(P, Q);
2269 : }
2270 938 : return mkvec4(P,Q,FF,mkvecsmall2(dP,dQ));
2271 : }
2272 :
2273 : /* Number of solutions of x^2+b*x+c */
2274 : static long
2275 896 : F2xqX_quad_nbroots(GEN b, GEN c, GEN T)
2276 : {
2277 896 : if (lgpol(b) > 0)
2278 : {
2279 154 : GEN d = F2xq_div(c, F2xq_sqr(b, T), T);
2280 154 : return F2xq_trace(d, T)? 0: 2;
2281 : }
2282 : else
2283 742 : return 1;
2284 : }
2285 :
2286 : static GEN
2287 938 : genus2_eulerfact2_semistable(GEN PQ)
2288 : {
2289 938 : GEN V = F2x_genus_red(ZX_to_F2x(gel(PQ, 1)), ZX_to_F2x(gel(PQ, 2)));
2290 938 : GEN P = gel(V, 1), Q = gel(V, 2);
2291 938 : GEN F = gel(V, 3), v = gel(V, 4);
2292 : GEN abe, tor;
2293 938 : GEN ki, kp = pol_1(0), kq = pol_1(0);
2294 938 : long dP = F2x_degree(P), dQ = F2x_degree(Q), d = maxss(dP, 2*dQ);
2295 938 : if (!lgpol(F)) return pol_1(0);
2296 840 : ki = dQ!=0 || dP>0 ? xdminusone(1):
2297 42 : dP==-1 ? ZX_sqr(xdminusone(1)): xdminusone(2);
2298 1596 : abe = d>=5? hyperellcharpoly(gmul(PQ,gmodulss(1,2))):
2299 798 : d>=3? ellfromeqncharpoly(F2x_to_ZX(P), F2x_to_ZX(Q), gen_2):
2300 644 : pol_1(0);
2301 798 : if (lgpol(F))
2302 : {
2303 798 : GEN M = gel(F2x_factor(F), 1);
2304 798 : long i, lF = lg(M)-1;
2305 1694 : for(i=1; i <= lF; i++)
2306 : {
2307 896 : GEN Fi = gel(M, i);
2308 896 : long d = F2x_degree(Fi);
2309 896 : long nb = F2xqX_quad_nbroots(F2x_rem(Q, Fi), F2x_rem(P, Fi), Fi);
2310 1050 : GEN kqf = nb==1 ? xdminusone(d):
2311 56 : nb==2 ? ZX_sqr(xdminusone(d))
2312 154 : : xdminusone(2*d);
2313 896 : kp = gmul(kp, xdminusone(d));
2314 896 : kq = gmul(kq, kqf);
2315 : }
2316 : }
2317 798 : if (maxss(v[1],2*v[2])<5)
2318 : {
2319 287 : GEN kqoo = v[1]>2*v[2] ? xdminusone(1):
2320 7 : v[1]<2*v[2] ? ZX_sqr(xdminusone(1))
2321 21 : : xdminusone(2);
2322 266 : kp = gmul(kp, xdminusone(1));
2323 266 : kq = gmul(kq, kqoo);
2324 : }
2325 798 : tor = RgX_div(ZX_mul(xdminusone(1),kq), ZX_mul(ki, kp));
2326 798 : return ZX_mul(abe, tor);
2327 : }
2328 :
2329 : GEN
2330 1127 : genus2_eulerfact2(GEN PQ)
2331 : {
2332 1127 : pari_sp av = avma;
2333 1127 : GEN W, R = genus2_type5ns_2(gel(PQ,1), gel(PQ,2), gen_2), E;
2334 1127 : if (R) return R;
2335 910 : W = hyperellextremalmodels_i(PQ, 2, gen_2);
2336 910 : if (lg(W) < 3) return genus2_eulerfact2_semistable(PQ);
2337 28 : E = gmul(genus2_eulerfact2_semistable(gel(W,1)),
2338 28 : genus2_eulerfact2_semistable(gel(W,2)));
2339 28 : return gc_upto(av, E);
2340 : }
2341 :
2342 : GEN
2343 1806 : genus2charpoly(GEN G, GEN p)
2344 : {
2345 1806 : pari_sp av = avma;
2346 1806 : GEN gr = genus2red(G, p), F;
2347 1806 : GEN PQ = gel(gr, 3), L = gel(gr, 4), r = gel(L, 4);
2348 1806 : if (equaliu(p,2))
2349 910 : F = genus2_eulerfact2(PQ);
2350 : else
2351 : {
2352 896 : GEN P = gadd(gsqr(gel(PQ, 2)), gmul2n(gel(PQ, 1), 2));
2353 896 : F = genus2_eulerfact(P,p, r[1],r[2]);
2354 : }
2355 1806 : return gc_upto(av, F);
2356 : }
2357 :
2358 : /****************************************************************************/
2359 : /** **/
2360 : /** hyperellisisom **/
2361 : /** **/
2362 : /****************************************************************************/
2363 :
2364 : /* Based on a magma script isgl2equiv.m from
2365 : R. Lercier, C. Ritzenthaler & J. Sijsling
2366 : https://github.com/JRSijsling/hyperelliptic/blob/main/magma/toolbox/isgl2equiv.m
2367 : based on the paper
2368 : R. Lercier, C. Ritzenthaler & J. Sijsling
2369 : Fast computation of isomorphisms of hyperelliptic curves and explicit Galois descent.
2370 : ANTS X pages 463-486. Mathematical Sciences Publishers, 2013.
2371 : https://msp.org/obs/2013/1-1/obs-v1-n1-p23-s.pdf
2372 : https://arxiv.org/pdf/1203.5440v1
2373 : */
2374 :
2375 : static long
2376 3759 : hyperelldegree(GEN f)
2377 3759 : { long d = degpol(f); return d + odd(d); }
2378 :
2379 : static GEN
2380 2681 : hyperellchangevar(GEN P, GEN M, long v)
2381 : {
2382 2681 : long d = hyperelldegree(P);
2383 2681 : GEN A = deg1pol_shallow(gcoeff(M,1,1), gcoeff(M,1,2), v);
2384 2681 : GEN B = deg1pol_shallow(gcoeff(M,2,1), gcoeff(M,2,2), v);
2385 2681 : return RgX_homogenous_eval(P, A, B, d);
2386 : }
2387 :
2388 : static GEN
2389 2107 : checkisom(GEN nf, GEN f, GEN g, GEN s)
2390 : {
2391 2107 : GEN fs = gdiv(hyperellchangevar(f, s, varn(f)), g), z;
2392 2107 : if (typ(fs)!=t_POL || degpol(fs)!=0) return NULL;
2393 2107 : if (!nf) return issquareall(gel(fs,2), &z) ? z: NULL;
2394 945 : return nfissquare(nf, gel(fs,2), &z) ? nf_to_scalar_or_polmod(nf, z):NULL;
2395 : }
2396 :
2397 : static GEN
2398 539 : hyperellisisom_M0(GEN nf, GEN f, GEN g)
2399 : {
2400 539 : pari_sp av = avma;
2401 : GEN EQ2, EQ3, PG;
2402 539 : long d = hyperelldegree(f), d2 = d*d, dg = degpol(g);
2403 539 : GEN a0 = gel(f, 2), a1 = gel(f, 3), a2 = gel(f, 4), a3 = gel(f, 5);
2404 539 : GEN bm0 = dg<d ? gen_0: gel(g, d+2), bm1 = gel(g, d+1), bm2 = gel(g, d), bm3 = gel(g, d-1);
2405 : GEN L, R;
2406 539 : long l, i, k = 1;
2407 : GEN a1_2, a0_2, a0_3, bm0_2, bm0_3;
2408 539 : if (gequal0(a0))
2409 147 : return NULL;
2410 392 : a1_2 = gsqr(a1); a0_2 = gsqr(a0); a0_3 = gmul(a0,a0_2);
2411 392 : bm0_2 = gsqr(bm0); bm0_3 = gmul(bm0, bm0_2);
2412 392 : EQ2 = mkpoln(3,
2413 : gmul(bm0_2, gadd(gmulsg(1-d, a1_2), gmul(gmulgs(gmulsg(2, a2), d), a0))),
2414 : gen_0,
2415 : gmul(gneg(a0_2), gadd(gmul(gmulgs(gmulsg(2, bm2), d), bm0),
2416 : gmulsg(1-d, gsqr(bm1)))));
2417 1568 : EQ3 = mkpoln(4,
2418 : gmul(gmulsg(2, bm0_3), gadd(gmul(gadd(gmul(gmulsg(3*d2, a3), a0),
2419 392 : gmul(gmulsg(6*d-3*d2, a2), a1)), a0), gmulsg(d2-3*d+2, gpowgs(a1, 3)))),
2420 392 : gmul(gmul(gmul(gadd(gmul(gmulsg(6*d2-12*d, a2), a0),
2421 392 : gmulsg(-3*d2+9*d-6, a1_2)), a0), bm1), gsqr(bm0)),
2422 : gen_0,
2423 392 : gmul(gadd(gmul(gmulsg(-6*d2, bm3), gsqr(bm0)), gmulsg(d2-3*d+2, gpowgs(bm1, 3))), a0_3));
2424 392 : PG = RgX_gcd(EQ2, EQ3);
2425 392 : if (gequal0(PG)) return NULL;
2426 322 : if (degpol(PG)==0) retgc_const(av, cgetg(1, t_VEC));
2427 70 : R = nfroots(nf, PG); l = lg(R);
2428 70 : L = cgetg(l, t_VEC);
2429 140 : for (i = 1; i < l; i++)
2430 : {
2431 70 : GEN B = gel(R, i), D, M;
2432 70 : if (gequal0(B))
2433 35 : continue;
2434 35 : D = gdiv(gsub(gmul(bm1, a0), gmul(gmul(a1, B), bm0)), gmulgs(gmul(bm0, a0), d));
2435 35 : M = mkmat22(gen_0, B, gen_1, D);
2436 35 : if (checkisom(nf, f, g, M))
2437 28 : gel(L, k++) = M;
2438 : }
2439 70 : setlg(L, k);
2440 70 : return gc_GEN(av, L);
2441 : }
2442 :
2443 : static GEN
2444 2331 : RgX_hom_evaly(GEN F, GEN y, long d)
2445 2331 : { return poleval(RgXn_recip_shallow(F, d+1), y); }
2446 :
2447 : #define Dy(a,b,c) RgX_homogenous_derivn(a,b,c)
2448 :
2449 : static GEN
2450 539 : hyperellisisom_gen(GEN nf, GEN f1, GEN f2)
2451 : {
2452 539 : GEN F2 = f2;
2453 539 : long x = varn(f2);
2454 : GEN Nm2, Nm3, Nm4, EQ1, EQ2, PG;
2455 : GEN M12, bm0, bm2, bm3, F1, dF1, d2F1, d3F1, d4F1, L, R;
2456 : GEN dF1_2, dF1_3, dF1_4;
2457 : GEN dyF1, dyF1_2, dyF1_3, dyF1_4;
2458 539 : long lR, i, Li = 1;
2459 539 : long d2 = degpol(f2), d = hyperelldegree(f1);
2460 539 : M12 = gel(f2, d2+1);
2461 539 : if (!gequal0(M12))
2462 : {
2463 322 : M12 = gdiv(gneg(M12), gmulsg(d2, gel(f2, 2+d2)));
2464 322 : f2 = RgX_Rg_translate(f2, M12);
2465 : }
2466 539 : bm0 = d2 < d ? gen_0: gel(f2, d + 2);
2467 539 : bm2 = gel(f2, d);
2468 539 : bm3 = gel(f2, d - 1);
2469 539 : F1 = f1;
2470 539 : dF1 = RgX_deriv(F1); d2F1 = RgX_deriv(dF1); d3F1 = RgX_deriv(d2F1); d4F1 = RgX_deriv(d3F1);
2471 539 : dF1_2 = gsqr(dF1); dF1_3 = gmul(dF1_2, dF1); dF1_4 = gsqr(dF1_2);
2472 539 : dyF1 = Dy(F1, 1, d); dyF1_2 = gsqr(dyF1); dyF1_3 = gmul(dyF1_2, dyF1); dyF1_4 = gsqr(dyF1_2);
2473 539 : Nm2 = gadd(gsub(gmul(d2F1, dyF1_2), gmul(gmul(gmulsg(2, dF1), dyF1), Dy(dF1, 1, d-1))), gmul(Dy(F1, 2, d), dF1_2));
2474 539 : Nm2 = RgX_div(RgXn_recip_shallow(Nm2, 3*d-4+1), RgXn_recip_shallow(F1, d+1));
2475 539 : if (d > 3)
2476 : {
2477 539 : GEN bm4 = gel(f2, d - 2);
2478 539 : Nm4 = gadd(gsub(gadd(gsub(gmul(d4F1, dyF1_4),
2479 : gmul(gmul(gmulsg(4, Dy(d3F1, 1, d-3)), dyF1_3), dF1)),
2480 : gmul(gmul(gmulsg(6, Dy(d2F1, 2, d-2)), dyF1_2), dF1_2)),
2481 : gmul(gmul(gmulsg(4, Dy(dF1, 3, d-1)), dyF1), dF1_3)),
2482 : gmul(Dy(F1, 4, d), dF1_4));
2483 539 : Nm4 = RgX_div(RgXn_recip_shallow(Nm4, 5*d - 8+1), RgXn_recip_shallow(F1, d+1));
2484 539 : EQ1 = gsub(gmul(gmul(gmulsg(6, bm0), bm4), gsqr(Nm2)), gmul(gsqr(bm2), Nm4));
2485 : }
2486 539 : Nm3 = gadd(gsub(gadd(gmul(gneg(d3F1), dyF1_3),
2487 : gmul(gmul(gmulsg(3, Dy(d2F1, 1, d-2)), dyF1_2), dF1)),
2488 : gmul(gmul(gmulsg(3, Dy(dF1, 2, d-1)), dyF1), dF1_2)),
2489 : gmul(Dy(F1, 3, d), dF1_3));
2490 539 : Nm3 = RgX_div(RgXn_recip_shallow(Nm3, 4*d - 6+1), RgXn_recip_shallow(F1, d+1));
2491 539 : EQ2 = gsub(gmul(gmul(gmulsg(9, bm0), gsqr(bm3)), gpowgs(Nm2, 3)),
2492 : gmul(gmulsg(2, gpowgs(bm2, 3)), gsqr(Nm3)));
2493 :
2494 539 : PG = d > 3 ? RgX_gcd(EQ1, EQ2): EQ2;
2495 539 : if (gequal0(PG)) return NULL;
2496 336 : if (lgpol(PG)==0) return cgetg(1, t_VEC);
2497 336 : R = nfroots(nf, PG); lR = lg(R);
2498 336 : L = cgetg(lR, t_VEC);
2499 1498 : for (i = 1; i < lR; i++)
2500 : {
2501 1169 : long j, LRi = 1, k=1, lR, g;
2502 : GEN gU, U, m, Rd;
2503 : GEN degs, LR, g1, m12;
2504 1169 : GEN m21 = gel(R, i), d12 = RgX_hom_evaly(dF1, m21, d - 1), N, D;
2505 1169 : if (gequal0(d12))
2506 7 : return NULL;
2507 1162 : m12 = gdiv(RgX_hom_evaly(Dy(F1, 1, d), m21, d - 1), gneg(d12));
2508 1162 : g1 = RgX_homogenous_eval(F1, deg1pol_shallow(gen_1, m12, x), deg1pol_shallow(m21, gen_1, x), d);
2509 7665 : for (j = 1; j < d; j++)
2510 : {
2511 6524 : GEN am1 = gel(g1, j+1), am0 = gel(g1, j+2), ap1 = gel(g1, j+3);
2512 6524 : GEN bm1 = gel(f2, j+1), bm0 = gel(f2, j+2), bp1 = gel(f2, j+3);
2513 6524 : if (!gequal(gmul(gmul(am1, ap1), gsqr(bm0)), gmul(gmul(bm1, bp1), gsqr(am0))))
2514 21 : break;
2515 : }
2516 1162 : if (j < d) continue;
2517 1141 : degs = cgetg(d, t_VEC);
2518 7644 : for (j = 2; j <= d; ++j)
2519 6503 : if (!gequal0(gel(f2, d-j+2)))
2520 5621 : gel(degs, k++) = stoi(j);
2521 1141 : setlg(degs,k);
2522 1141 : gU = mathnf0(degs, 1);
2523 1141 : g = itos(gmael3(gU,1, 1, 1));
2524 1141 : U = gel(gU, 2); U = vec_to_vecsmall(gel(U, lg(U)-1));
2525 1141 : degs = vec_to_vecsmall(degs);
2526 9716 : for (j = 2; j <= d+2; j++)
2527 8610 : if (gequal0(gel(g1, j)) != gequal0(gel(f2, j)))
2528 35 : break;
2529 1141 : if (j < d) continue;
2530 1106 : m = gdiv(gel(g1, d+2), gel(f2, d+2));
2531 1106 : N = gen_1; D = gen_1;
2532 6552 : for (j = 1; j < k; ++j)
2533 : {
2534 5446 : long c = 2+d-degs[j];
2535 5446 : N = gmul(N, gpowgs(gmul(m, gel(f2, c)), U[j]));
2536 5446 : D = gmul(D, gpowgs(gel(g1, c), U[j]));
2537 : }
2538 1106 : Rd = nfroots(nf, gsub(gmul(D, pol_xn(g, x)), N)); lR = lg(Rd);
2539 1106 : LR = cgetg(lR, t_VEC);
2540 2436 : for (j = 1; j < lR; j++)
2541 : {
2542 1330 : GEN RL = gel(Rd, j);
2543 1330 : GEN M = mkmat22(gen_1, gsub(gmul(m12,RL),M12), m21, gsub(RL,gmul(M12,m21)));
2544 1330 : if (checkisom(nf, f1, F2, M))
2545 903 : gel(LR, LRi++) = M;
2546 : }
2547 1106 : setlg(LR, LRi);
2548 1106 : gel(L, Li++) = LR;
2549 : }
2550 329 : setlg(L,Li);
2551 329 : return Li>1 ? shallowconcat1(L): cgetg(1,t_VEC);
2552 : }
2553 :
2554 : #undef Dy
2555 :
2556 : static GEN
2557 581 : random_SL2(GEN B)
2558 : {
2559 : GEN a, b, d, u, v;
2560 581 : do { a = randomi(B); b = randomi(B); } while (signe(a)==0 && signe(b)==0);
2561 574 : d = bezout(a,b,&u,&v);
2562 574 : if(!is_pm1(d)) { a = diviiexact(a,d); b = diviiexact(b,d); }
2563 574 : return mkmat22(a, b, v, negi(u));
2564 : }
2565 :
2566 : static GEN
2567 322 : nfM_primpart(GEN nf, GEN M)
2568 : {
2569 322 : pari_sp av = avma;
2570 322 : GEN d, id = nfV_idealhnf(nf, RgM_flatten_RgC(M), &d);
2571 322 : GEN c = gel(idealred(nf, mkvec2(id, gen_1)), 2);
2572 322 : GEN A = gdiv(M, nf_to_scalar_or_polmod(nf, c));
2573 322 : return gc_upto(av, d ? gmul(A,d): A);
2574 : }
2575 :
2576 : static GEN
2577 252 : nfhyperellisisom(GEN nf, GEN W1, GEN W2)
2578 : {
2579 252 : pari_sp av = avma, av2;
2580 : long i, v, g;
2581 252 : GEN f1, f2, F1, F2, M1 = NULL, M2 = NULL;
2582 252 : if (nf)
2583 : {
2584 : GEN u;
2585 84 : checknf(nf);
2586 84 : u = gmodulo(gen_1, nf_get_pol(nf));
2587 84 : W1 = gmul(W1, u);
2588 84 : W2 = gmul(W2, u);
2589 : }
2590 252 : check_hyperell_Rg("hyperellisisom", &W1, &F1); f1 = F1;
2591 252 : check_hyperell_Rg("hyperellisisom", &W2, &F2); f2 = F2;
2592 252 : g = hyperellgenus(F1); v = varn(F1);
2593 252 : if (g < 1) pari_err_DOMAIN("hyperellisisom","genus(C1)","<",gen_1,W1);
2594 252 : if (hyperellgenus(F2) != g) return cgetg(1,t_VEC);
2595 252 : av2 = avma;
2596 252 : for (i = 10; ; i++)
2597 287 : {
2598 539 : GEN Q0 = hyperellisisom_M0(nf, f1, f2);
2599 539 : GEN Qp = hyperellisisom_gen(nf, f1, f2);
2600 539 : if (Q0 && Qp)
2601 : {
2602 252 : GEN Q = shallowconcat(Q0,Qp);
2603 252 : long i, l = lg(Q);
2604 252 : GEN V = cgetg(2*l-1, t_VEC);
2605 994 : for (i = 1; i < l; i++)
2606 : {
2607 742 : GEN Mr = !M1 ? gel(Q,i): RgM_mul(RgM_mul(M1, gel(Q, i)), M2);
2608 742 : GEN M = nf ? nfM_primpart(nf, Mr): Q_primpart(Mr);
2609 742 : GEN e = checkisom(nf, F1, F2, M);
2610 742 : gel(V,2*i-1) = hyperellisom_finalize(W1, W2, e, M, g, v);
2611 742 : gel(V,2*i) = hyperellisom_finalize(W1, W2, gneg(e), M, g, v);
2612 : }
2613 252 : return gc_GEN(av, V);
2614 : }
2615 287 : set_avma(av2);
2616 287 : M1 = random_SL2(stoi(-i));
2617 287 : f1 = hyperellchangevar(F1, M1, v);
2618 287 : M2 = random_SL2(stoi(-i));
2619 287 : f2 = hyperellchangevar(F2, ginv(M2), v);
2620 : }
2621 : }
2622 :
2623 : GEN
2624 84 : hyperellisisom(GEN W1, GEN W2, GEN nf)
2625 84 : { return nfhyperellisisom(nf, W1, W2); }
2626 :
2627 : GEN
2628 168 : hyperellauto(GEN W, GEN nf)
2629 168 : { return nfhyperellisisom(nf, W, W); }
|
