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 : /* Implementation of Kedlaya Algorithm for counting point on hyperelliptic
26 : curves by Bill Allombert based on a GP script by Bernadette Perrin-Riou.
27 :
28 : References:
29 : Pierrick Gaudry and Nicolas G\"urel
30 : Counting Points in Medium Characteristic Using Kedlaya's Algorithm
31 : Experiment. Math. Volume 12, Number 4 (2003), 395-402.
32 : http://projecteuclid.org/euclid.em/1087568016
33 :
34 : Harrison, M. An extension of Kedlaya's algorithm for hyperelliptic
35 : curves. Journal of Symbolic Computation, 47 (1) (2012), 89-101.
36 : http://arxiv.org/pdf/1006.4206v3.pdf
37 : */
38 :
39 : /* We use the basis of differentials (x^i*dx/y^k) (i=1 to 2*g-1),
40 : with k either 1 or 3, depending on p and d, see Harrison paper */
41 :
42 : static long
43 1764 : get_basis(long p, long d)
44 : {
45 1764 : if (odd(d))
46 868 : return p < d-1 ? 3 : 1;
47 : else
48 896 : return 2*p <= d-2 ? 3 : 1;
49 : }
50 :
51 : static GEN
52 20265 : FpXXQ_red(GEN S, GEN T, GEN p)
53 : {
54 20265 : pari_sp av = avma;
55 20265 : long i, dS = degpol(S);
56 : GEN A, C;
57 20265 : if (signe(S)==0) return pol_0(varn(T));
58 20265 : A = cgetg(dS+3, t_POL);
59 20265 : C = pol_0(varn(T));
60 1520393 : for(i=dS; i>0; i--)
61 : {
62 1500128 : GEN Si = FpX_add(C, gel(S,i+2), p);
63 1500128 : GEN R, Q = FpX_divrem(Si, T, p, &R);
64 1500128 : gel(A,i+2) = R;
65 1500128 : C = Q;
66 : }
67 20265 : gel(A,2) = FpX_add(C, gel(S,2), p);
68 20265 : A[1] = S[1];
69 20265 : return gc_GEN(av, FpXX_renormalize(A,dS+3));
70 : }
71 :
72 : static GEN
73 3402 : FpXXQ_sqr(GEN x, GEN T, GEN p)
74 : {
75 3402 : pari_sp av = avma;
76 3402 : long n = degpol(T);
77 3402 : GEN z = FpX_red(ZXX_sqr_Kronecker(x, n), p);
78 3402 : z = Kronecker_to_ZXX(z, n, varn(T));
79 3402 : return gc_upto(av, FpXXQ_red(z, T, p));
80 : }
81 :
82 : static GEN
83 16863 : FpXXQ_mul(GEN x, GEN y, GEN T, GEN p)
84 : {
85 16863 : pari_sp av = avma;
86 16863 : long n = degpol(T);
87 16863 : GEN z = FpX_red(ZXX_mul_Kronecker(x, y, n), p);
88 16863 : z = Kronecker_to_ZXX(z, n, varn(T));
89 16863 : return gc_upto(av, FpXXQ_red(z, T, p));
90 : }
91 :
92 : static GEN
93 1309 : ZpXXQ_invsqrt(GEN S, GEN T, ulong p, long e)
94 : {
95 1309 : pari_sp av = avma, av2;
96 : ulong mask;
97 1309 : long v = varn(S), n=1;
98 1309 : GEN a = pol_1(v);
99 1309 : if (e <= 1) return gc_GEN(av, a);
100 1309 : mask = quadratic_prec_mask(e);
101 1309 : av2 = avma;
102 4676 : for (;mask>1;)
103 : {
104 : GEN q, q2, q22, f, fq, afq;
105 3367 : long n2 = n;
106 3367 : n<<=1; if (mask & 1) n--;
107 3367 : mask >>= 1;
108 3367 : q = powuu(p,n); q2 = powuu(p,n2);
109 3367 : f = RgX_sub(FpXXQ_mul(FpXX_red(S, q), FpXXQ_sqr(a, T, q), T, q), pol_1(v));
110 3367 : fq = ZXX_Z_divexact(f, q2);
111 3367 : q22 = shifti(addiu(q2,1),-1);
112 3367 : afq = FpXX_Fp_mul(FpXXQ_mul(a, fq, T, q2), q22, q2);
113 3367 : a = RgX_sub(a, ZXX_Z_mul(afq, q2));
114 3367 : if (gc_needed(av2,1))
115 : {
116 0 : if(DEBUGMEM>1) pari_warn(warnmem,"ZpXXQ_invsqrt, e = %ld", n);
117 0 : a = gc_upto(av2, a);
118 : }
119 : }
120 1309 : return gc_upto(av, a);
121 : }
122 :
123 : static GEN
124 1029007 : to_ZX(GEN a, long v) { return typ(a)==t_INT? scalarpol(a,v): a; }
125 :
126 : static void
127 14 : is_sing(GEN H, ulong p)
128 : {
129 14 : pari_err_DOMAIN("hyperellpadicfrobenius","H","is singular at",utoi(p),H);
130 0 : }
131 :
132 : static void
133 1309 : get_UV(GEN *U, GEN *V, GEN T, ulong p, long e)
134 : {
135 1309 : GEN q = powuu(p,e), d;
136 1309 : GEN dT = FpX_deriv(T, q);
137 1309 : GEN R = polresultantext(T, dT);
138 1309 : long v = varn(T);
139 1309 : if (dvdiu(gel(R,3),p)) is_sing(T, p);
140 1309 : d = Zp_inv(gel(R,3), utoi(p), e);
141 1309 : *U = FpX_Fp_mul(FpX_red(to_ZX(gel(R,1),v),q),d,q);
142 1309 : *V = FpX_Fp_mul(FpX_red(to_ZX(gel(R,2),v),q),d,q);
143 1309 : }
144 :
145 : static GEN
146 133847 : frac_to_Fp(GEN a, GEN b, GEN p)
147 : {
148 133847 : GEN d = gcdii(a, b);
149 133847 : return Fp_div(diviiexact(a, d), diviiexact(b, d), p);
150 : }
151 :
152 : static GEN
153 10094 : ZpXXQ_frob(GEN S, GEN U, GEN V, long k, GEN T, ulong p, long e)
154 : {
155 10094 : pari_sp av = avma, av2;
156 10094 : long i, pr = degpol(S), dT = degpol(T), vT = varn(T);
157 10094 : GEN q = powuu(p,e);
158 10094 : GEN Tp = FpX_deriv(T, q), Tp1 = RgX_shift_shallow(Tp, 1);
159 10094 : GEN M = to_ZX(gel(S,pr+2),vT), R;
160 10094 : av2 = avma;
161 987868 : for(i = pr-1; i>=k; i--)
162 : {
163 : GEN A, B, H, Bc;
164 : ulong v, r;
165 977774 : H = FpX_divrem(FpX_mul(V,M,q), T, q, &B);
166 977774 : A = FpX_add(FpX_mul(U,M,q), FpX_mul(H, Tp, q),q);
167 977774 : v = u_lvalrem(2*i+1,p,&r);
168 977774 : Bc = ZX_deriv(B);
169 977774 : Bc = FpX_Fp_mul(ZX_divuexact(Bc,upowuu(p,v)),Fp_divu(gen_2, r, q), q);
170 977774 : M = FpX_add(to_ZX(gel(S,i+2),vT), FpX_add(A, Bc, q), q);
171 977774 : if (gc_needed(av2,1))
172 : {
173 0 : if(DEBUGMEM>1) pari_warn(warnmem,"ZpXXQ_frob, step 1, i = %ld", i);
174 0 : M = gc_upto(av2, M);
175 : }
176 : }
177 10094 : if (degpol(M)<dT-1)
178 5488 : return gc_upto(av, M);
179 4606 : R = RgX_shift_shallow(M,dT-degpol(M)-2);
180 4606 : av2 = avma;
181 237629 : for(i = degpol(M)-dT+2; i>=1; i--)
182 : {
183 : GEN B, c;
184 233023 : R = RgX_shift_shallow(R, 1);
185 233023 : gel(R,2) = gel(M, i+1);
186 233023 : if (degpol(R) < dT) continue;
187 130935 : B = FpX_add(FpX_mulu(T, 2*i, q), Tp1, q);
188 130935 : c = frac_to_Fp(leading_coeff(R), leading_coeff(B), q);
189 130935 : R = FpX_sub(R, FpX_Fp_mul(B, c, q), q);
190 130935 : if (gc_needed(av2,1))
191 : {
192 0 : if(DEBUGMEM>1) pari_warn(warnmem,"ZpXXQ_frob, step 2, i = %ld", i);
193 0 : R = gc_upto(av2, R);
194 : }
195 : }
196 4606 : if (degpol(R)==dT-1)
197 : {
198 2912 : GEN c = frac_to_Fp(leading_coeff(R), leading_coeff(Tp), q);
199 2912 : R = FpX_sub(R, FpX_Fp_mul(Tp, c, q), q);
200 2912 : return gc_upto(av, R);
201 : } else
202 1694 : return gc_GEN(av, R);
203 : }
204 :
205 : static GEN
206 12026 : revdigits(GEN v)
207 : {
208 12026 : long i, n = lg(v)-1;
209 12026 : GEN w = cgetg(n+2, t_POL);
210 12026 : w[1] = evalsigne(1)|evalvarn(0);
211 168784 : for (i=0; i<n; i++)
212 156758 : gel(w,i+2) = gel(v,n-i);
213 12026 : return FpXX_renormalize(w, n+2);
214 : }
215 :
216 : static GEN
217 10094 : diff_red(GEN s, GEN A, long m, GEN T, GEN p)
218 : {
219 10094 : long v, n, vT = varn(T);
220 : GEN Q, sQ, qS;
221 : pari_timer ti;
222 10094 : if (DEBUGLEVEL>1) timer_start(&ti);
223 10094 : Q = revdigits(FpX_digits(A,T,p));
224 10094 : n = degpol(Q);
225 10094 : if (DEBUGLEVEL>1) timer_printf(&ti,"reddigits");
226 10094 : sQ = FpXXQ_mul(s,Q,T,p);
227 10094 : if (DEBUGLEVEL>1) timer_printf(&ti,"redmul");
228 10094 : qS = RgX_shift_shallow(sQ,m-n);
229 10094 : v = ZX_val(sQ);
230 10094 : if (n > m + v)
231 : {
232 4564 : long i, l = n-m-v;
233 4564 : GEN rS = cgetg(l+1,t_VEC);
234 29190 : for (i = l-1; i >=0 ; i--)
235 24626 : gel(rS,i+1) = to_ZX(gel(sQ, 1+v+l-i), vT);
236 4564 : rS = FpXV_FpX_fromdigits(rS,T,p);
237 4564 : gel(qS,2) = FpX_add(FpX_mul(rS, T, p), gel(qS, 2), p);
238 4564 : if (DEBUGLEVEL>1) timer_printf(&ti,"redadd");
239 : }
240 10094 : return qS;
241 : }
242 :
243 : static GEN
244 10094 : ZC_to_padic(GEN C, GEN q)
245 : {
246 10094 : long i, l = lg(C);
247 10094 : GEN V = cgetg(l,t_COL);
248 102914 : for(i = 1; i < l; i++)
249 92820 : gel(V, i) = gadd(gel(C, i), q);
250 10094 : return V;
251 : }
252 :
253 : static GEN
254 1309 : ZM_to_padic(GEN M, GEN q)
255 : {
256 1309 : long i, l = lg(M);
257 1309 : GEN V = cgetg(l,t_MAT);
258 11403 : for(i = 1; i < l; i++)
259 10094 : gel(V, i) = ZC_to_padic(gel(M, i), q);
260 1309 : return V;
261 : }
262 :
263 : static GEN
264 1743 : ZX_to_padic(GEN P, GEN q)
265 : {
266 1743 : long i, l = lg(P);
267 1743 : GEN Q = cgetg(l, t_POL);
268 1743 : Q[1] = P[1];
269 5978 : for (i=2; i<l ;i++)
270 4235 : gel(Q,i) = gadd(gel(P,i), q);
271 1743 : return normalizepol(Q);
272 : }
273 :
274 : static GEN
275 469 : ZXC_to_padic(GEN x, GEN q)
276 2212 : { pari_APPLY_type(t_COL, ZX_to_padic(gel(x, i), q)) }
277 :
278 : static GEN
279 147 : ZXM_to_padic(GEN x, GEN q)
280 616 : { pari_APPLY_same(ZXC_to_padic(gel(x, i), q)) }
281 :
282 : static GEN
283 1309 : ZlX_hyperellpadicfrobenius(GEN H, ulong p, long n)
284 : {
285 1309 : pari_sp av = avma;
286 : long k, N, i, d;
287 : GEN F, s, Q, pN1, U, V;
288 : pari_timer ti;
289 1309 : if (typ(H) != t_POL) pari_err_TYPE("hyperellpadicfrobenius",H);
290 1309 : if (p == 2) is_sing(H, 2);
291 1309 : d = degpol(H);
292 1309 : if (d <= 0)
293 0 : pari_err_CONSTPOL("hyperellpadicfrobenius");
294 1309 : if (n < 1)
295 0 : pari_err_DOMAIN("hyperellpadicfrobenius","n","<", gen_1, utoi(n));
296 1309 : k = get_basis(p, d);
297 1309 : N = n + ulogint(2*n, p) + 1;
298 1309 : pN1 = powuu(p,N+1);
299 1309 : Q = RgX_to_FpX(H, pN1);
300 1309 : if (dvdiu(leading_coeff(Q),p)) is_sing(H, p);
301 1309 : setvarn(Q,1);
302 1309 : if (DEBUGLEVEL>1) timer_start(&ti);
303 1309 : s = revdigits(FpX_digits(RgX_inflate(Q, p), Q, pN1));
304 1309 : if (DEBUGLEVEL>1) timer_printf(&ti,"s1");
305 1309 : s = ZpXXQ_invsqrt(s, Q, p, N);
306 1309 : if (k==3)
307 35 : s = FpXXQ_mul(s, FpXXQ_sqr(s, Q, pN1), Q, pN1);
308 1309 : if (DEBUGLEVEL>1) timer_printf(&ti,"invsqrt");
309 1309 : get_UV(&U, &V, Q, p, N+1);
310 1309 : F = cgetg(d, t_MAT);
311 11403 : for (i = 1; i < d; i++)
312 : {
313 10094 : pari_sp av2 = avma;
314 : GEN M, D;
315 10094 : D = diff_red(s, monomial(utoipos(p),p*i-1,1),(k*p-1)>>1, Q, pN1);
316 10094 : if (DEBUGLEVEL>1) timer_printf(&ti,"red");
317 10094 : M = ZpXXQ_frob(D, U, V, (k-1)>>1, Q, p, N + 1);
318 10094 : if (DEBUGLEVEL>1) timer_printf(&ti,"frob");
319 10094 : gel(F, i) = gc_GEN(av2, RgX_to_RgC(M, d-1));
320 : }
321 1309 : return gc_upto(av, F);
322 : }
323 :
324 : GEN
325 1309 : hyperellpadicfrobenius(GEN H, ulong p, long n)
326 : {
327 1309 : pari_sp av = avma;
328 1309 : GEN M = ZlX_hyperellpadicfrobenius(H, p, n);
329 1309 : GEN q = zeropadic_shallow(utoipos(p),n);
330 1309 : return gc_upto(av, ZM_to_padic(M, q));
331 : }
332 :
333 : INLINE GEN
334 2247 : FpXXX_renormalize(GEN x, long lx) { return ZXX_renormalize(x,lx); }
335 :
336 : static GEN
337 1806 : ZpXQXXQ_red(GEN F, GEN S, GEN T, GEN q, GEN p, long e)
338 : {
339 1806 : pari_sp av = avma;
340 1806 : long i, dF = degpol(F);
341 : GEN A, C;
342 1806 : if (signe(F)==0) return pol_0(varn(S));
343 1806 : A = cgetg(dF+3, t_POL);
344 1806 : C = pol_0(varn(S));
345 96404 : for(i=dF; i>0; i--)
346 : {
347 94598 : GEN Fi = FpXX_add(C, gel(F,i+2), q);
348 94598 : GEN R, Q = ZpXQX_divrem(Fi, S, T, q, p, e, &R);
349 94598 : gel(A,i+2) = R;
350 94598 : C = Q;
351 : }
352 1806 : gel(A,2) = FpXX_add(C, gel(F,2), q);
353 1806 : A[1] = F[1];
354 1806 : return gc_GEN(av, FpXXX_renormalize(A,dF+3));
355 : }
356 :
357 : static GEN
358 448 : ZpXQXXQ_sqr(GEN x, GEN S, GEN T, GEN q, GEN p, long e)
359 : {
360 448 : pari_sp av = avma;
361 : GEN z, kx;
362 448 : long n = degpol(S), vx = varn(S);
363 448 : kx = RgXX_to_Kronecker_var(x, n, vx);
364 448 : z = Kronecker_to_ZXX(FpXQX_sqr(kx, T, q), n, vx);
365 448 : setvarn(z, varn(x));
366 448 : return gc_upto(av, ZpXQXXQ_red(z, S, T, q, p, e));
367 : }
368 :
369 : static GEN
370 1358 : ZpXQXXQ_mul(GEN x, GEN y, GEN S, GEN T, GEN q, GEN p, long e)
371 : {
372 1358 : pari_sp av = avma;
373 : GEN z, kx, ky;
374 1358 : long n = degpol(S), vx = varn(S);
375 1358 : kx = RgXX_to_Kronecker_var(x, n, vx);
376 1358 : ky = RgXX_to_Kronecker_var(y, n, vx);
377 1358 : z = Kronecker_to_ZXX(FpXQX_mul(ky, kx, T, q), n, vx);
378 1358 : setvarn(z, varn(x));
379 1358 : return gc_upto(av, ZpXQXXQ_red(z, S, T, q, p, e));
380 : }
381 :
382 : static GEN
383 441 : FpXXX_red(GEN z, GEN p)
384 : {
385 : GEN res;
386 441 : long i, l = lg(z);
387 441 : res = cgetg(l,t_POL); res[1] = z[1];
388 17388 : for (i=2; i<l; i++)
389 : {
390 16947 : GEN zi = gel(z,i);
391 16947 : if (typ(zi)==t_INT)
392 287 : gel(res,i) = modii(zi,p);
393 : else
394 16660 : gel(res,i) = FpXX_red(zi,p);
395 : }
396 441 : return FpXXX_renormalize(res,lg(res));
397 : }
398 :
399 : static GEN
400 441 : FpXXX_Fp_mul(GEN z, GEN a, GEN p)
401 : {
402 441 : return FpXXX_red(RgX_Rg_mul(z, a), p);
403 : }
404 :
405 : static GEN
406 154 : ZpXQXXQ_invsqrt(GEN F, GEN S, GEN T, ulong p, long e)
407 : {
408 154 : pari_sp av = avma, av2, av3;
409 : ulong mask;
410 154 : long v = varn(F), n=1;
411 : pari_timer ti;
412 154 : GEN a = pol_1(v), pp = utoipos(p);
413 154 : if (DEBUGLEVEL>1) timer_start(&ti);
414 154 : if (e <= 1) return gc_GEN(av, a);
415 154 : mask = quadratic_prec_mask(e);
416 154 : av2 = avma;
417 595 : for (;mask>1;)
418 : {
419 : GEN q, q2, q22, f, fq, afq;
420 441 : long n2 = n;
421 441 : n<<=1; if (mask & 1) n--;
422 441 : mask >>= 1;
423 441 : q = powuu(p,n); q2 = powuu(p,n2);
424 441 : av3 = avma;
425 441 : f = RgX_sub(ZpXQXXQ_mul(F, ZpXQXXQ_sqr(a, S, T, q, pp, n), S, T, q, pp, n), pol_1(v));
426 441 : fq = gc_upto(av3, RgX_Rg_divexact(f, q2));
427 441 : q22 = shifti(addiu(q2,1),-1);
428 441 : afq = FpXXX_Fp_mul(ZpXQXXQ_mul(a, fq, S, T, q2, pp, n2), q22, q2);
429 441 : a = RgX_sub(a, RgX_Rg_mul(afq, q2));
430 441 : if (gc_needed(av2,1))
431 : {
432 0 : if(DEBUGMEM>1) pari_warn(warnmem,"ZpXQXXQ_invsqrt, e = %ld", n);
433 0 : a = gc_upto(av2, a);
434 : }
435 : }
436 154 : return gc_upto(av, a);
437 : }
438 :
439 : static GEN
440 6573 : frac_to_Fq(GEN a, GEN b, GEN T, GEN q, GEN p, long e)
441 : {
442 6573 : GEN d = gcdii(ZX_content(a), ZX_content(b));
443 6573 : return ZpXQ_div(ZX_Z_divexact(a, d), ZX_Z_divexact(b, d), T, q, p, e);
444 : }
445 :
446 : static GEN
447 469 : ZpXQXXQ_frob(GEN F, GEN U, GEN V, long k, GEN S, GEN T, ulong p, long e)
448 : {
449 469 : pari_sp av = avma, av2;
450 469 : long i, pr = degpol(F), dS = degpol(S), v = varn(T);
451 469 : GEN q = powuu(p,e), pp = utoipos(p);
452 469 : GEN Sp = RgX_deriv(S), Sp1 = RgX_shift_shallow(Sp, 1);
453 469 : GEN M = gel(F,pr+2), R;
454 469 : av2 = avma;
455 52311 : for(i = pr-1; i>=k; i--)
456 : {
457 : GEN A, B, H, Bc;
458 : ulong v, r;
459 51842 : H = ZpXQX_divrem(FpXQX_mul(V, M, T, q), S, T, q, utoipos(p), e, &B);
460 51842 : A = FpXX_add(FpXQX_mul(U, M, T, q), FpXQX_mul(H, Sp, T, q),q);
461 51842 : v = u_lvalrem(2*i+1,p,&r);
462 51842 : Bc = RgX_deriv(B);
463 51842 : Bc = FpXX_Fp_mul(ZXX_Z_divexact(Bc,powuu(p,v)), Fp_divu(gen_2, r, q), q);
464 51842 : M = FpXX_add(gel(F,i+2), FpXX_add(A, Bc, q), q);
465 51842 : if (gc_needed(av2,1))
466 : {
467 0 : if(DEBUGMEM>1) pari_warn(warnmem,"ZpXQXXQ_frob, step 1, i = %ld", i);
468 0 : M = gc_upto(av2, M);
469 : }
470 : }
471 469 : if (degpol(M)<dS-1)
472 266 : return gc_upto(av, M);
473 203 : R = RgX_shift_shallow(M,dS-degpol(M)-2);
474 203 : av2 = avma;
475 7175 : for(i = degpol(M)-dS+2; i>=1; i--)
476 : {
477 : GEN B, c;
478 6972 : R = RgX_shift_shallow(R, 1);
479 6972 : gel(R,2) = gel(M, i+1);
480 6972 : if (degpol(R) < dS) continue;
481 6412 : B = FpXX_add(FpXX_mulu(S, 2*i, q), Sp1, q);
482 6412 : c = frac_to_Fq(to_ZX(leading_coeff(R),v), to_ZX(leading_coeff(B),v), T, q, pp, e);
483 6412 : R = FpXX_sub(R, FpXQX_FpXQ_mul(B, c, T, q), q);
484 6412 : if (gc_needed(av2,1))
485 : {
486 0 : if(DEBUGMEM>1) pari_warn(warnmem,"ZpXXQ_frob, step 2, i = %ld", i);
487 0 : R = gc_upto(av2, R);
488 : }
489 : }
490 203 : if (degpol(R)==dS-1)
491 : {
492 161 : GEN c = frac_to_Fq(to_ZX(leading_coeff(R),v), to_ZX(leading_coeff(Sp),v), T, q, pp, e);
493 161 : R = FpXX_sub(R, FpXQX_FpXQ_mul(Sp, c, T, q), q);
494 161 : return gc_upto(av, R);
495 : } else
496 42 : return gc_GEN(av, R);
497 : }
498 :
499 : static GEN
500 469 : Fq_diff_red(GEN s, GEN A, long m, GEN S, GEN T, GEN q, GEN p, long e)
501 : {
502 : long v, n;
503 : GEN Q, sQ, qS;
504 : pari_timer ti;
505 469 : if (DEBUGLEVEL>1) timer_start(&ti);
506 469 : Q = revdigits(ZpXQX_digits(A, S, T, q, p, e));
507 469 : n = degpol(Q);
508 469 : if (DEBUGLEVEL>1) timer_printf(&ti,"reddigits");
509 469 : sQ = ZpXQXXQ_mul(s, Q, S, T, q, p, e);
510 469 : if (DEBUGLEVEL>1) timer_printf(&ti,"redmul");
511 469 : qS = RgX_shift_shallow(sQ,m-n);
512 469 : v = ZX_val(sQ);
513 469 : if (n > m + v)
514 : {
515 189 : long i, l = n-m-v;
516 189 : GEN rS = cgetg(l+1,t_VEC);
517 1547 : for (i = l-1; i >=0 ; i--)
518 1358 : gel(rS,i+1) = gel(sQ, 1+v+l-i);
519 189 : rS = FpXQXV_FpXQX_fromdigits(rS, S, T, q);
520 189 : gel(qS,2) = FpXX_add(FpXQX_mul(rS, S, T, q), gel(qS, 2), q);
521 189 : if (DEBUGLEVEL>1) timer_printf(&ti,"redadd");
522 : }
523 469 : return qS;
524 : }
525 :
526 : static void
527 154 : Fq_get_UV(GEN *U, GEN *V, GEN S, GEN T, ulong p, long e)
528 : {
529 154 : GEN q = powuu(p, e), pp = utoipos(p), d;
530 154 : GEN dS = RgX_deriv(S), R = polresultantext(S, dS), C;
531 154 : long v = varn(S);
532 154 : if (signe(FpX_red(to_ZX(gel(R,3),v), pp))==0) is_sing(S, p);
533 147 : C = FpXQ_red(to_ZX(gel(R, 3),v), T, q);
534 147 : d = ZpXQ_inv(C, T, pp, e);
535 147 : *U = FpXQX_FpXQ_mul(FpXQX_red(to_ZX(gel(R,1),v),T,q),d,T,q);
536 147 : *V = FpXQX_FpXQ_mul(FpXQX_red(to_ZX(gel(R,2),v),T,q),d,T,q);
537 147 : }
538 :
539 : static GEN
540 469 : ZXX_to_FpXC(GEN x, long N, GEN p, long v)
541 : {
542 : long i, l;
543 : GEN z;
544 469 : l = lg(x)-1; x++;
545 469 : if (l > N+1) l = N+1; /* truncate higher degree terms */
546 469 : z = cgetg(N+1,t_COL);
547 2170 : for (i=1; i<l ; i++)
548 : {
549 1701 : GEN xi = gel(x, i);
550 1701 : gel(z,i) = typ(xi)==t_INT? scalarpol(Fp_red(xi, p), v): FpX_red(xi, p);
551 : }
552 511 : for ( ; i<=N ; i++)
553 42 : gel(z,i) = pol_0(v);
554 469 : return z;
555 : }
556 :
557 : GEN
558 154 : ZlXQX_hyperellpadicfrobenius(GEN H, GEN T, ulong p, long n)
559 : {
560 154 : pari_sp av = avma;
561 : long k, N, i, d, N1, v0;
562 : GEN xp, F, s, q, Q, pN1, U, V, pp;
563 : pari_timer ti;
564 154 : if (typ(H) != t_POL) pari_err_TYPE("hyperellpadicfrobenius",H);
565 154 : if (p == 2) is_sing(H, 2);
566 154 : d = degpol(H);
567 154 : if (d <= 0) pari_err_CONSTPOL("hyperellpadicfrobenius");
568 154 : if (n < 1) pari_err_DOMAIN("hyperellpadicfrobenius","n","<", gen_1, utoi(n));
569 154 : k = get_basis(p, d); pp = utoipos(p);
570 154 : N = n + ulogint(2*n, p) + 1;
571 154 : q = powuu(p,n); N1 = N+1;
572 154 : pN1 = powuu(p,N1); T = FpX_get_red(T, pN1);
573 154 : Q = RgX_to_FqX(H, T, pN1);
574 154 : if (signe(FpX_red(to_ZX(leading_coeff(Q),varn(Q)),pp))==0) is_sing(H, p);
575 154 : if (DEBUGLEVEL>1) timer_start(&ti);
576 154 : xp = ZpX_Frobenius(T, pp, N1);
577 154 : s = RgX_inflate(FpXY_FpXQ_evalx(Q, xp, T, pN1), p);
578 154 : v0 = fetch_var_higher();
579 154 : s = revdigits(ZpXQX_digits(s, Q, T, pN1, pp, N1));
580 154 : setvarn(s, v0);
581 154 : if (DEBUGLEVEL>1) timer_printf(&ti,"s1");
582 154 : s = ZpXQXXQ_invsqrt(s, Q, T, p, N);
583 154 : if (k==3)
584 7 : s = ZpXQXXQ_mul(s, ZpXQXXQ_sqr(s, Q, T, pN1, pp, N1), Q, T, pN1, pp, N1);
585 154 : if (DEBUGLEVEL>1) timer_printf(&ti,"invsqrt");
586 154 : Fq_get_UV(&U, &V, Q, T, p, N+1);
587 147 : if (DEBUGLEVEL>1) timer_printf(&ti,"get_UV");
588 147 : F = cgetg(d, t_MAT);
589 616 : for (i = 1; i < d; i++)
590 : {
591 469 : pari_sp av2 = avma;
592 : GEN M, D;
593 469 : D = Fq_diff_red(s, monomial(pp,p*i-1,0),(k*p-1)>>1, Q, T, pN1, pp, N1);
594 469 : if (DEBUGLEVEL>1) timer_printf(&ti,"red");
595 469 : M = ZpXQXXQ_frob(D, U, V, (k - 1)>>1, Q, T, p, N1);
596 469 : if (DEBUGLEVEL>1) timer_printf(&ti,"frob");
597 469 : gel(F, i) = gc_upto(av2, ZXX_to_FpXC(M, d-1, q, varn(T)));
598 : }
599 147 : delete_var();
600 147 : return gc_upto(av, F);
601 : }
602 :
603 : GEN
604 154 : nfhyperellpadicfrobenius(GEN H, GEN T, ulong p, long n)
605 : {
606 154 : pari_sp av = avma;
607 154 : GEN pp = utoipos(p), q = zeropadic_shallow(pp, n);
608 154 : GEN M = ZlXQX_hyperellpadicfrobenius(lift_shallow(H),T,p,n);
609 147 : GEN MM = ZpXQM_prodFrobenius(M, T, pp, n);
610 147 : GEN m = gmul(ZXM_to_padic(MM, q), gmodulo(gen_1, T));
611 147 : return gc_upto(av, m);
612 : }
613 :
614 : GEN
615 595 : hyperellpadicfrobenius0(GEN H, GEN Tp, long n)
616 : {
617 : GEN T, p;
618 595 : if (!ff_parse_Tp(Tp, &T,&p,0)) pari_err_TYPE("hyperellpadicfrobenius", Tp);
619 595 : if (lgefint(p) > 3) pari_err_IMPL("large prime in hyperellpadicfrobenius");
620 7 : return T? nfhyperellpadicfrobenius(H, T, itou(p), n)
621 602 : : hyperellpadicfrobenius(H, itou(p), n);
622 : }
623 :
624 : static GEN
625 392 : F2x_genus2charpoly_naive(GEN P, GEN Q)
626 : {
627 392 : long a, b = 1, c = 0;
628 392 : GEN T = mkvecsmall2(P[1], 7);
629 392 : GEN PT = F2x_rem(P, T), QT = F2x_rem(Q, T);
630 392 : long q0 = F2x_eval(Q, 0), q1 = F2x_eval(Q, 1);
631 392 : long dP = F2x_degree(P), dQ = F2x_degree(Q);
632 392 : a= dQ<3 ? 0: dP<=5 ? 1: -1;
633 392 : a += (q0? F2x_eval(P, 0)? -1: 1: 0) + (q1? F2x_eval(P, 1)? -1: 1: 0);
634 392 : b += q0 + q1;
635 392 : if (lgpol(QT))
636 308 : c = (F2xq_trace(F2xq_div(PT, F2xq_sqr(QT, T), T), T)==0 ? 1: -1);
637 392 : return mkvecsmalln(6, 0UL, 4UL, 2*a, (b+2*c+a*a)>>1, a, 1UL);
638 : }
639 :
640 : static GEN
641 19544 : Flx_difftable(GEN P, ulong p)
642 : {
643 19544 : long i, n = degpol(P);
644 19544 : GEN V = cgetg(n+2, t_VEC);
645 19544 : gel(V, n+1) = P;
646 132832 : for(i = n; i >= 1; i--)
647 113288 : gel(V, i) = Flx_diff1(gel(V, i+1), p);
648 19544 : return V;
649 : }
650 :
651 : static GEN
652 766955 : FlxV_Fl2_eval_pre(GEN V, GEN x, ulong D, ulong p, ulong pi)
653 : {
654 766955 : long i, n = lg(V)-1;
655 766955 : GEN r = cgetg(n+1, t_VEC);
656 5981661 : for (i = 1; i <= n; i++)
657 5214706 : gel(r, i) = Flx_Fl2_eval_pre(gel(V, i), x, D, p, pi);
658 766955 : return r;
659 : }
660 :
661 : static GEN
662 97965182 : Fl2V_next(GEN V, ulong p)
663 : {
664 97965182 : long i, n = lg(V)-1;
665 97965182 : GEN r = cgetg(n+1, t_VEC);
666 97965182 : gel(r, 1) = gel(V, 1);
667 666067388 : for (i = 2; i <= n; i++)
668 568102206 : gel(r, i) = Flv_add(gel(V, i), gel(V, i-1), p);
669 97965182 : return r;
670 : }
671 :
672 : static GEN
673 19544 : FlxV_constant(GEN x)
674 152376 : { pari_APPLY_long(Flx_constant(gel(x,i))) }
675 :
676 : static GEN
677 19544 : Flx_genus2charpoly_naive(GEN H, ulong p)
678 : {
679 19544 : pari_sp av = avma, av2;
680 19544 : ulong pi = get_Fl_red(p);
681 19544 : ulong i, j, p2 = p>>1, D = 2, e = ((p&2UL) == 0) ? -1 : 1;
682 19544 : long a, b, c = 0, n = degpol(H);
683 19544 : GEN t, d, k = const_vecsmall(p, -1);
684 19544 : k[1] = 0;
685 786499 : for (i=1, j=1; i < p; i += 2, j = Fl_add(j, i, p)) k[j+1] = 1;
686 35679 : while (k[1+D] >= 0) D++;
687 19544 : b = n == 5 ? 0 : 1;
688 19544 : a = b ? k[1+Flx_lead(H)]: 0;
689 19544 : t = Flx_difftable(H, p);
690 19544 : d = FlxV_constant(t);
691 19544 : av2 = avma;
692 1572998 : for (i=0; i < p; i++)
693 : {
694 1553454 : ulong v = uel(d,n+1);
695 1553454 : a += k[1+v];
696 1553454 : b += !!v;
697 1553454 : if (n==6)
698 1241520 : uel(d,7) = Fl_add(uel(d,7), uel(d,6), p);
699 1553454 : uel(d,6) = Fl_add(uel(d,6), uel(d,5), p);
700 1553454 : uel(d,5) = Fl_add(uel(d,5), uel(d,4), p);
701 1553454 : uel(d,4) = Fl_add(uel(d,4), uel(d,3), p);
702 1553454 : uel(d,3) = Fl_add(uel(d,3), uel(d,2), p);
703 1553454 : uel(d,2) = Fl_add(uel(d,2), uel(d,1), p);
704 : }
705 786499 : for (j=1; j <= p2; j++)
706 : {
707 766955 : GEN V = FlxV_Fl2_eval_pre(t, mkvecsmall2(0, j), D, p, pi);
708 766955 : for (i=0;; i++)
709 97965182 : {
710 98732137 : GEN r2 = gel(V, n+1);
711 197464274 : c += uel(r2,2) ?
712 97786633 : (uel(r2,1) ? uel(k,1+Fl2_norm_pre(r2, D, p, pi)): e)
713 196518770 : : !!uel(r2,1);
714 98732137 : if (i == p-1) break;
715 97965182 : V = Fl2V_next(V, p);
716 : }
717 766955 : set_avma(av2);
718 : }
719 19544 : set_avma(av);
720 19544 : return mkvecsmalln(6, 0UL, p*p, a*p, (b+2*c+a*a)>>1, a, 1UL);
721 : }
722 :
723 : static GEN
724 679 : charpoly_funceq(GEN P, GEN q)
725 : {
726 679 : long i, l, g = degpol(P)>>1;
727 679 : GEN R, Q = gpowers0(q, g-1, q); /* Q[i] = q^i, i <= g */
728 679 : R = cgetg_copy(P, &l); R[1] = P[1];
729 3164 : for (i=0; i<g; i++) gel(R, i+2) = mulii(gel(P, 2*g-i+2), gel(Q, g-i));
730 3843 : for (; i<=2*g; i++) gel(R, i+2) = icopy(gel(P, i+2));
731 679 : return R;
732 : }
733 :
734 : static long
735 686 : hyperell_Weil_bound(GEN q, ulong g, GEN p)
736 : {
737 686 : pari_sp av = avma;
738 686 : GEN w = mulii(binomialuu(2*g,g),sqrtint(shifti(powiu(q, g),2)));
739 686 : return gc_long(av, logint(w,p) + 1);
740 : }
741 :
742 : /* return 4P + Q^2 */
743 : static GEN
744 401783 : check_hyperell(GEN PQ)
745 : {
746 : GEN H;
747 401783 : if (is_vec_t(typ(PQ)) && lg(PQ)==3)
748 332114 : H = gadd(gsqr(gel(PQ, 2)), gmul2n(gel(PQ, 1), 2));
749 : else
750 69669 : H = gmul2n(PQ, 2);
751 401783 : return typ(H) == t_POL? H: NULL;
752 : }
753 :
754 : static long
755 535427 : hyperellgenus(GEN H)
756 535427 : { long d = degpol(H); return ((d+1)>>1)-1; }
757 :
758 : static void
759 152922 : check_hyperell_Rg(const char *fun, GEN *pW, GEN *pF)
760 : {
761 152922 : GEN W = *pW, F = check_hyperell(W);
762 : long v;
763 152922 : if (!F)
764 7 : pari_err_TYPE(fun, W);
765 152915 : if (degpol(F) <= 0) pari_err_CONSTPOL(fun);
766 152908 : v = varn(F);
767 152908 : if (typ(W)==t_POL) W = mkvec2(W, pol_0(v));
768 : else
769 : {
770 152880 : GEN P = gel(W, 1), Q = gel(W, 2);
771 152880 : long g = hyperellgenus(F);
772 152880 : if( typ(P)!=t_POL) P = scalarpol(P, v);
773 152880 : if( typ(Q)!=t_POL) Q = scalarpol(Q, v);
774 152880 : if (degpol(P) > 2*g+2)
775 0 : pari_err_DOMAIN(fun, "poldegree(P)", ">", utoi(2*g+2), P);
776 152880 : if (degpol(Q) > g+1)
777 0 : pari_err_DOMAIN(fun, "poldegree(Q)", ">", utoi(g+1), Q);
778 :
779 152880 : W = mkvec2(P, Q);
780 : }
781 152908 : if (pF) *pF = F;
782 152908 : *pW = W;
783 152908 : }
784 :
785 : GEN
786 20629 : hyperellcharpoly(GEN PQ)
787 : {
788 20629 : pari_sp av = avma;
789 20629 : GEN M, R, T=NULL, pp=NULL, q;
790 20629 : long d, n, eps = 0;
791 : ulong p;
792 20629 : GEN H = check_hyperell(PQ);
793 20629 : if (!H || !RgX_is_FpXQX(H, &T, &pp) || !pp)
794 0 : pari_err_TYPE("hyperellcharpoly", PQ);
795 20629 : p = itou(pp);
796 20629 : if (!T)
797 : {
798 20482 : if (p==2 && is_vec_t(typ(PQ)))
799 : {
800 392 : long dP, dQ, v = varn(H);
801 392 : GEN P = gel(PQ,1), Q = gel(PQ,2);
802 392 : if (typ(P)!=t_POL) P = scalarpol(P, v);
803 392 : if (typ(Q)!=t_POL) Q = scalarpol(Q, v);
804 392 : dP = degpol(P); dQ = degpol(Q);
805 392 : if (dP<=6 && dQ <=3 && (dQ==3 || dP>=5))
806 : {
807 392 : GEN P2 = RgX_to_F2x(P), Q2 = RgX_to_F2x(Q);
808 392 : GEN D = F2x_add(F2x_mul(P2, F2x_sqr(F2x_deriv(Q2))), F2x_sqr(F2x_deriv(P2)));
809 392 : if (F2x_degree(F2x_gcd(D, Q2))) is_sing(PQ, 2);
810 392 : if (dP==6 && dQ<3 && F2x_coeff(P2,5)==F2x_coeff(Q2,2))
811 0 : is_sing(PQ, 2); /* The curve is singular at infinity */
812 392 : R = zx_to_ZX(F2x_genus2charpoly_naive(P2, Q2));
813 392 : return gc_upto(av, R);
814 : }
815 : }
816 20090 : H = RgX_to_FpX(H, pp);
817 20090 : d = degpol(H);
818 20090 : if (d <= 0) is_sing(H, p);
819 20090 : if (p > 2 && ((d == 5 && p < 17500) || (d == 6 && p < 24500)))
820 : {
821 19551 : GEN Hp = ZX_to_Flx(H, p);
822 19551 : if (!Flx_is_squarefree(Hp, p)) is_sing(H, p);
823 19544 : R = zx_to_ZX(Flx_genus2charpoly_naive(Hp, p));
824 19544 : return gc_upto(av, R);
825 : }
826 539 : n = hyperell_Weil_bound(pp, (d-1)>>1, pp);
827 539 : eps = odd(d)? 0: Fp_issquare(leading_coeff(H), pp);
828 539 : M = hyperellpadicfrobenius(H, p, n);
829 539 : R = centerlift(carberkowitz(M, 0));
830 539 : q = pp;
831 : }
832 : else
833 : {
834 : int fixvar;
835 147 : T = typ(T)==t_FFELT? FF_mod(T): RgX_to_FpX(T, pp);
836 147 : q = powuu(p, degpol(T));
837 147 : fixvar = (varncmp(varn(T),varn(H)) <= 0);
838 147 : if (fixvar) setvarn(T, fetch_var());
839 147 : H = RgX_to_FpXQX(H, T, pp);
840 147 : d = degpol(H);
841 147 : if (d <= 0) is_sing(H, p);
842 147 : eps = odd(d)? 0: Fq_issquare(leading_coeff(H), T, pp);
843 147 : n = hyperell_Weil_bound(q, (d-1)>>1, pp);
844 147 : M = nfhyperellpadicfrobenius(H, T, p, n);
845 140 : R = simplify_shallow(centerlift(liftpol_shallow(carberkowitz(M, 0))));
846 140 : if (fixvar) (void)delete_var();
847 : }
848 679 : if (!odd(d))
849 : {
850 301 : GEN b = get_basis(p, d) == 3 ? gen_1 : q;
851 301 : GEN pn = powuu(p, n);
852 301 : R = FpX_div_by_X_x(R, eps? b: negi(b), pn, NULL);
853 301 : R = FpX_center_i(R, pn, shifti(pn,-1));
854 : }
855 679 : return gc_upto(av, charpoly_funceq(R, q));
856 : }
857 :
858 : GEN
859 70 : hyperellordinate(GEN W, GEN x)
860 : {
861 70 : pari_sp av = avma;
862 70 : if (typ(W)==t_POL)
863 : {
864 : GEN d, y;
865 42 : if (typ(x)==t_INFINITY)
866 : {
867 14 : long dW = degpol(W);
868 14 : d = odd(dW) ? gen_0: gel(W,dW+2);
869 : } else
870 28 : d = poleval(W,x);
871 42 : if (gequal0(d)) { return gc_GEN(av, mkvec(d)); }
872 35 : if (!issquareall(d, &y)) retgc_const(av, cgetg(1, t_VEC));
873 14 : return gc_GEN(av, mkvec2(y, gneg(y)));
874 : }
875 : else
876 : {
877 : GEN b, c, d, rd, y, P, Q, F;
878 28 : check_hyperell_Rg("hyperellordinate", &W, &F);
879 28 : P = gel(W,1); Q = gel(W,2);
880 28 : if (typ(x)==t_INFINITY)
881 : {
882 7 : long dP = degpol(P), dQ = degpol(Q), g = hyperellgenus(F);
883 7 : c = dP < 2*g+2 ? gen_0: gel(P,dP+2);
884 7 : b = dQ < g+1 ? gen_0: gel(Q,dQ+2);
885 : } else
886 21 : { b = poleval(Q, x); c = poleval(P, x); }
887 28 : d = gadd(gsqr(b), gmul2n(c, 2));
888 28 : if (gequal0(d)) { return gc_GEN(av, mkvec(gmul2n(gneg(b),-1))); }
889 21 : if (!issquareall(d, &rd)) retgc_const(av, cgetg(1, t_VEC));
890 14 : y = gmul2n(gsub(rd, b), -1);
891 14 : return gc_GEN(av, mkvec2(y, gsub(y,rd)));
892 : }
893 : }
894 :
895 : GEN
896 120721 : hyperelldisc(GEN PQ)
897 : {
898 120721 : pari_sp av = avma;
899 120721 : GEN D, H = check_hyperell(PQ);
900 : long g;
901 120721 : if (!H || signe(H)==0) pari_err_TYPE("hyperelldisc",PQ);
902 120721 : g = hyperellgenus(H);
903 120721 : D = gmul2n(RgX_disc(H),-4*(g+1));
904 120721 : if (odd(degpol(H))) D = gmul(D, gsqr(leading_coeff(H)));
905 120721 : return gc_upto(av, D);
906 : }
907 :
908 : static long
909 130093 : get_ep(GEN W)
910 : {
911 130093 : GEN P = gel(W,1), Q = gel(W,2);
912 130093 : if (signe(Q)==0) return ZX_lval(P,2);
913 89814 : return minss(ZX_lval(P,2), ZX_lval(Q,2));
914 : }
915 :
916 : static GEN
917 52534 : algo51(GEN W, GEN M)
918 : {
919 52534 : GEN P = gel(W,1), Q = gel(W,2);
920 : for(;;)
921 10647 : {
922 63181 : long vP = ZX_lval(P,2);
923 63181 : long vQ = signe(Q) ? ZX_lval(Q,2): vP+1;
924 : long r;
925 : /* 1 */
926 63181 : if (vQ==0) break;
927 : /* 2 */
928 36302 : if (vP==0)
929 : {
930 : GEN H, H1;
931 : /* a */
932 29714 : RgX_even_odd(FpX_red(P,gen_2),&H, &H1);
933 29714 : if (signe(H1)) break;
934 : /* b */
935 15035 : P = ZX_add(P, ZX_mul(H, ZX_sub(Q, H)));
936 15035 : Q = ZX_sub(Q, ZX_shifti(H, 1));
937 15035 : vP = ZX_lval(P,2);
938 15035 : vQ = signe(Q) ? ZX_lval(Q,2): vP+1;
939 : }
940 : /* 2c */
941 21623 : if (vP==1) break;
942 : /* 2d */
943 10647 : r = minss(2*vQ, vP)>>1;
944 10647 : if (M) gel(M,1) = shifti(gel(M,1), r);
945 10647 : P = ZX_shifti(P, -2*r);
946 10647 : Q = ZX_shifti(Q, -r);
947 : }
948 52534 : return mkvec2(P,Q);
949 : }
950 :
951 : static GEN
952 105506 : algo52(GEN W, GEN c, long *pt_lambda)
953 : {
954 : long lambda;
955 105506 : GEN P = gel(W,1), Q = gel(W,2);
956 : for(;;)
957 118338 : {
958 : GEN H, H1;
959 : /* 1 */
960 223844 : GEN Pc = ZX_affine(P,gen_2,c), Qc = ZX_affine(Q,gen_2,c);
961 223844 : long mP = ZX_lval(Pc,2), mQ = signe(Qc) ? ZX_lval(Qc,2): mP+1;
962 : /* 2 */
963 223844 : if (2*mQ <= mP) { lambda = 2*mQ; break; }
964 : /* 3 */
965 191223 : if (odd(mP)) { lambda = mP; break; }
966 : /* 4 */
967 129265 : RgX_even_odd(FpX_red(ZX_shifti(Pc, -mP),gen_2),&H, &H1);
968 129265 : if (signe(H1)) { lambda = mP; break; }
969 : /* 5 */
970 118338 : P = ZX_add(P, ZX_mul(H, ZX_sub(Q, H)));
971 118338 : Q = ZX_sub(Q, ZX_shifti(H, 1));
972 : }
973 105506 : *pt_lambda = lambda;
974 105506 : return mkvec2(P,Q);
975 : }
976 :
977 : static long
978 149115 : test53(long lambda, long ep, long g)
979 : {
980 149115 : return (lambda <= g+1) || (odd(g) && lambda<g+3 && ep==1);
981 : }
982 :
983 : static long
984 194492 : test55(GEN W, long ep, long g)
985 : {
986 194492 : GEN P = gel(W,1), Q = gel(W,2);
987 194492 : GEN Pe = FpX_red(ep ? ZX_shifti(P,-1): P, gen_2);
988 194492 : GEN Qe = FpX_red(ep ? ZX_shifti(Q,-1): Q, gen_2);
989 194492 : if (ep==0)
990 : {
991 153833 : if (signe(Qe)!=0) return ZX_val(Qe) >= (g + 3)>>1;
992 92620 : else return ZX_val(FpX_deriv(Pe, gen_2)) >= g+1;
993 : }
994 : else
995 40659 : return ZX_val(Qe) >= (g+1)>>1 && ZX_val(Pe) >= g + 1;
996 : }
997 :
998 : static GEN
999 52436 : hyperell_reverse(GEN W, long g)
1000 : {
1001 52436 : return mkvec2(RgXn_recip_shallow(gel(W,1),2*g+3),
1002 52436 : RgXn_recip_shallow(gel(W,2),g+2));
1003 : }
1004 :
1005 : /* [P,Q] -> [P(2x)/4^r, Q(2x)/2^r] */
1006 : static GEN
1007 173246 : ZX2_unscale(GEN W, long r)
1008 : {
1009 173246 : GEN P = ZX_unscale2n(gel(W,1), 1);
1010 173246 : GEN Q = ZX_unscale2n(gel(W,2), 1);
1011 173246 : if (r)
1012 : {
1013 30959 : P = ZX_shifti(P, -2*r);
1014 30959 : Q = ZX_shifti(Q, -r);
1015 : }
1016 173246 : return mkvec2(P,Q);
1017 : }
1018 : /* [P,Q] -> [P(2x+c)/4^r, Q(2x+c)/2^r] */
1019 : static GEN
1020 167032 : ZX2_affine_unscale(GEN W, long c, long r)
1021 : {
1022 243807 : if (c) W = mkvec2(ZX_Z_translate(gel(W,1), gen_1),
1023 76775 : ZX_Z_translate(gel(W,2), gen_1));
1024 167032 : return ZX2_unscale(W, r);
1025 : }
1026 :
1027 : static GEN
1028 52205 : algo56(GEN W, long g)
1029 : {
1030 : long ep;
1031 52205 : GEN M = mkvec2(gen_1, matid(2)), Woo;
1032 52205 : W = algo51(W, M);
1033 52205 : Woo = hyperell_reverse(W, g);
1034 52205 : ep = get_ep(Woo);
1035 52205 : if (test55(Woo,ep,g))
1036 : {
1037 : long lambda;
1038 11940 : Woo = algo52(Woo, gen_0, &lambda);
1039 11940 : if (!test53(lambda,ep,g))
1040 : {
1041 5969 : long r = lambda>>1;
1042 5969 : gel(M,1) = shifti(gel(M,1), r);
1043 5969 : gel(M,2) = ZM2_mul(gel(M,2), mkmat22(gen_0, gen_1, gen_2, gen_0));
1044 5969 : W = ZX2_unscale(Woo, r);
1045 : }
1046 : }
1047 : for(;;)
1048 24892 : {
1049 77097 : long j, ep = get_ep(W);
1050 193603 : for (j = 0; j < 2; j++)
1051 141398 : if (test55(ZX2_affine_unscale(W, j, 0), ep, g))
1052 : {
1053 : long lambda;
1054 93090 : GEN c = utoi(j), Wc = algo52(W, c, &lambda);
1055 93090 : if (!test53(lambda,ep,g))
1056 : {
1057 24892 : long r = lambda>>1;
1058 24892 : gel(M,1) = shifti(gel(M,1), r);
1059 24892 : gel(M,2) = ZM2_mul(gel(M,2), mkmat22(gen_2, c, gen_0, gen_1));
1060 24892 : W = ZX2_affine_unscale(Wc, j, r);
1061 24892 : break;
1062 : }
1063 : }
1064 77097 : if (j==2) break;
1065 : }
1066 52205 : return mkvec2(W, M);
1067 : }
1068 :
1069 : static GEN
1070 329 : algo56bis(GEN W, long g, long inf, long thr)
1071 : {
1072 329 : pari_sp av = avma;
1073 329 : GEN vl = cgetg(3,t_VEC);
1074 329 : long nl = 1;
1075 329 : W = algo51(W, NULL);
1076 329 : if (inf)
1077 : {
1078 231 : GEN Woo = hyperell_reverse(W, g);
1079 231 : long ep = get_ep(Woo);
1080 231 : if (test55(ZX2_unscale(Woo, 0), ep, g))
1081 : {
1082 : long lambda;
1083 70 : Woo = algo52(Woo, gen_0, &lambda);
1084 70 : if (lambda == thr) gel(vl,nl++) = ZX2_unscale(Woo, lambda>>1);
1085 : }
1086 : }
1087 : {
1088 329 : long j, ep = get_ep(W);
1089 987 : for (j = 0; j < 2; j++)
1090 658 : if (test55(ZX2_affine_unscale(W, j, 0), ep, g))
1091 : {
1092 : long lambda;
1093 406 : GEN Wc = algo52(W, utoi(j), &lambda);
1094 406 : if (lambda == thr) gel(vl,nl++) = ZX2_affine_unscale(Wc, j, lambda>>1);
1095 : }
1096 : }
1097 329 : setlg(vl, nl);
1098 329 : return gc_GEN(av,vl);
1099 : }
1100 :
1101 : /* return the (degree 2) apolar invariant (the nth transvectant of P and P) */
1102 : static GEN
1103 1491 : ZX_apolar(GEN P, long n)
1104 : {
1105 1491 : pari_sp av = avma;
1106 1491 : long d = degpol(P), i;
1107 1491 : GEN s = gen_0, g = cgetg(n+2,t_VEC);
1108 1491 : gel(g,1) = gen_1;
1109 10437 : for (i = 1; i <= n; i++) gel(g,i+1) = muliu(gel(g,i),i); /* g[i+1] = i! */
1110 11354 : for (i = n-d; i <= d; i++)
1111 : {
1112 9863 : GEN a = mulii(mulii(gel(g,i+1),gel(g,n-i+1)),
1113 9863 : mulii(gel(P,i+2),gel(P,n-i+2)));
1114 9863 : s = odd(i)? subii(s, a): addii(s, a);
1115 : }
1116 1491 : return gc_INT(av,s);
1117 : }
1118 :
1119 : static GEN
1120 54438 : algo57(GEN F, long g, GEN pr)
1121 : {
1122 : long i, l;
1123 54438 : GEN D, C = content(F);
1124 54438 : GEN e = gel(core2(shifti(C,-vali(C))),2);
1125 54438 : GEN M = mkvec2(e, matid(2));
1126 54438 : long minvd = (2*g+1)>>(odd(g) ? 4:2);
1127 54438 : F = ZX_Z_divexact(F, sqri(e));
1128 54438 : D = absi(hyperelldisc(F));
1129 54438 : if (!pr)
1130 : {
1131 1491 : GEN A = gcdii(D, ZX_apolar(F, 2*g+2));
1132 1491 : pr = gel(factor(shifti(A, -vali(A))),1);
1133 : }
1134 54438 : l = lg(pr);
1135 314809 : for (i = 1; i < l; i++)
1136 : {
1137 : long ep;
1138 260371 : GEN p = gel(pr, i), ps2 = shifti(p,-1), Fe;
1139 260371 : if (equaliu(p,2) || Z_pval(D,p) < minvd) continue;
1140 198135 : ep = ZX_pvalrem(F,p, &Fe); Fe = FpX_red(Fe, p);
1141 198135 : if (degpol(Fe) < g+1+ep)
1142 : {
1143 6406 : GEN Fi = ZX_unscale(RgXn_recip_shallow(F,2*g+3), p);
1144 6406 : long lambda = ZX_pval(Fi,p);
1145 6406 : if (!test53(lambda,ep,g))
1146 : {
1147 3815 : GEN ppr = powiu(p,lambda>>1);
1148 3815 : F = ZX_Z_divexact(Fi,sqri(ppr));
1149 3815 : gel(M,1) = mulii(gel(M,1), ppr);
1150 3815 : gel(M,2) = ZM2_mul(gel(M,2), mkmat22(gen_0,gen_1,p,gen_0));
1151 : }
1152 : }
1153 : for(;;)
1154 25186 : {
1155 : GEN Fe, R;
1156 223321 : long j, lR, ep = ZX_pvalrem(F,p, &Fe);
1157 223321 : R = FpX_roots_mult(FpX_red(Fe, p), g+2-ep, p); lR = lg(R);
1158 235814 : for (j = 1; j<lR; j++)
1159 : {
1160 37679 : GEN c = Fp_center(gel(R,j), p, ps2);
1161 37679 : GEN Fi = ZX_affine(F,p,c);
1162 37679 : long lambda = ZX_pval(Fi,p);
1163 37679 : if (!test53(lambda,ep,g))
1164 : {
1165 25186 : GEN ppr = powiu(p,lambda>>1);
1166 25186 : F = ZX_Z_divexact(Fi, sqri(ppr));
1167 25186 : gel(M,1) = mulii(gel(M,1), ppr);
1168 25186 : gel(M,2) = ZM2_mul(gel(M,2), mkmat22(p,c,gen_0,gen_1));
1169 25186 : break;
1170 : }
1171 : }
1172 223321 : if (j==lR) break;
1173 : }
1174 : }
1175 54438 : return mkvec2(F, M);
1176 : }
1177 :
1178 : /* if inf=0, ignore point at infinity */
1179 : static GEN
1180 3563 : algo57bis(GEN F, long g, GEN p, long inf, long thr)
1181 : {
1182 3563 : pari_sp av = avma;
1183 3563 : GEN vl = cgetg(3,t_VEC), Fe;
1184 3563 : long nl = 1, ep = ZX_pvalrem(F,p, &Fe);
1185 3563 : Fe = FpX_red(Fe, p);
1186 : {
1187 3563 : GEN R = FpX_roots_mult(Fe, thr-ep, p);
1188 3563 : long j, lR = lg(R);
1189 6496 : for (j = 1; j<lR; j++)
1190 : {
1191 2933 : GEN Fj = ZX_affine(F, p, gel(R,j));
1192 2933 : long lambda = ZX_pvalrem(Fj, p, &Fj);
1193 2933 : if (lambda == thr) gel(vl,nl++) = odd(lambda)? ZX_Z_mul(Fj, p): Fj;
1194 : }
1195 : }
1196 3563 : if (inf==1 && 2*g+2-degpol(Fe) >= thr-ep)
1197 : {
1198 0 : GEN Fj = ZX_unscale(RgXn_recip_shallow(F,2*g+3), p);
1199 0 : long lambda = ZX_pvalrem(Fj, p, &Fj);
1200 0 : if (lambda == thr) gel(vl,nl++) = odd(lambda)? ZX_Z_mul(Fj, p): Fj;
1201 : }
1202 3563 : setlg(vl, nl);
1203 3563 : return gc_GEN(av,vl);
1204 : }
1205 :
1206 : static GEN
1207 3892 : next_model(GEN G, long g, GEN p, long inf, long thr)
1208 : {
1209 4221 : return equaliu(p,2) ? algo56bis(G, g, inf, thr)
1210 4221 : : algo57bis(G, g, p, inf, thr);
1211 : }
1212 :
1213 : static GEN
1214 1526 : get_extremal_even(GEN F, GEN G, long g, GEN p, long *nb)
1215 : {
1216 : while (1)
1217 1274 : {
1218 1526 : GEN Wi = next_model(G, g, p, 0, g+2);
1219 1526 : if (lg(Wi)==1) return F;
1220 1365 : F = gel(Wi,1); ++*nb;
1221 1365 : if (DEBUGLEVEL>1) err_printf("model %ld: %Ps\n", *nb, F);
1222 1365 : Wi = next_model(F, g, p, 0, g+1);
1223 1365 : if (lg(Wi)==1) return F;
1224 1274 : G = gel(Wi,1);
1225 : }
1226 : }
1227 :
1228 : static GEN
1229 0 : get_extremal_odd(GEN F, long g, GEN p, long *nb)
1230 : {
1231 : while (1)
1232 0 : {
1233 0 : GEN Wi = next_model(F, g, p, 0, g+2);
1234 0 : if (lg(Wi)==1) return F;
1235 0 : F = gel(Wi,1); ++*nb;
1236 0 : if (DEBUGLEVEL>1) err_printf("model %ld: %Ps\n", *nb, F);
1237 : }
1238 : }
1239 :
1240 : static GEN
1241 1036 : hyperellextremalmodels_nb(GEN F, long g, GEN p, long *nb)
1242 : {
1243 1036 : pari_sp av = avma;
1244 : GEN W, A, B;
1245 : long l;
1246 :
1247 1036 : *nb = 1;
1248 1036 : if (equaliu(p,2))
1249 : {
1250 231 : if (get_ep(F) > 0) retmkvec(gcopy(F));
1251 : } else
1252 : {
1253 805 : F = check_hyperell(F);
1254 805 : if (ZX_pval(F, p) > 0) return gc_GEN(av, mkvec(F));
1255 : }
1256 1001 : if (DEBUGLEVEL>1) err_printf("model %ld: %Ps\n", *nb, F);
1257 1001 : W = next_model(F, g, p, 1, odd(g)? g+2: g+1);
1258 1001 : l = lg(W); if (l==1) return gc_GEN(av, mkvec(F));
1259 238 : if (odd(g))
1260 : {
1261 0 : *nb = l-1;
1262 0 : A = get_extremal_odd(gel(W,1), g, p, nb);
1263 0 : B = l==3 ? get_extremal_odd(gel(W,2), g, p, nb) : F;
1264 : }
1265 : else
1266 : {
1267 238 : A = get_extremal_even(F, gel(W,1), g, p, nb);
1268 238 : B = l==3 ? get_extremal_even(F, gel(W,2), g, p, nb) : F;
1269 : }
1270 238 : return gc_GEN(av, A == B? mkvec(A): mkvec2(A, B));
1271 : }
1272 :
1273 : static GEN
1274 1029 : hyperellextremalmodels_i(GEN F, long g, GEN p)
1275 : {
1276 : long nb;
1277 1029 : return hyperellextremalmodels_nb(F, g, p, &nb);
1278 : }
1279 :
1280 : GEN
1281 7 : hyperellextremalmodels(GEN PQ, GEN p)
1282 : {
1283 7 : pari_sp av = avma;
1284 7 : GEN H = check_hyperell(PQ), W, v;
1285 : long g, nb;
1286 7 : if (!H || signe(H)==0) pari_err_TYPE("hyperellextremalmodels",PQ);
1287 7 : if (typ(p)!=t_INT || signe(p)<=0) pari_err_TYPE("hyperellextremalmodels",p);
1288 7 : g = hyperellgenus(H);
1289 7 : W = hyperellminimalmodel(H,NULL,mkvec(p));
1290 7 : v = cgetg(3, t_VEC);
1291 7 : gel(v, 2) = hyperellextremalmodels_nb(W, g, p, &nb);
1292 7 : gel(v, 1) = stoi(nb);
1293 7 : return gc_upto(av, v);
1294 : }
1295 :
1296 : static GEN
1297 54452 : minimalmodel_merge(GEN W2, GEN Modd, long g, long v)
1298 : {
1299 54452 : GEN P = gel(W2,1), Q = gel(W2,2);
1300 54452 : GEN e = gel(Modd,1), M = gel(Modd,2);
1301 54452 : GEN A = deg1pol_shallow(gcoeff(M,1,1), gcoeff(M,1,2), v);
1302 54452 : GEN B = deg1pol_shallow(gcoeff(M,2,1), gcoeff(M,2,2), v);
1303 54452 : GEN Bp = gpowers(B, 2*g+2);
1304 54452 : long f = mod4(e)==1 ? 1: -1;
1305 54452 : GEN m = shifti(f > 0 ? subui(1,e): addui(1,e), -2);
1306 54452 : GEN m24 = subii(shifti(m,1), shifti(sqri(m),2));
1307 54452 : P = RgX_homogenous_evalpow(P, A, Bp, 2*g+2);
1308 54452 : Q = RgX_homogenous_evalpow(Q, A, Bp, g+1);
1309 54452 : P = ZX_Z_divexact(ZX_add(P, ZX_Z_mul(ZX_sqr(Q), m24)),sqri(e));
1310 54452 : if (f < 0) Q = ZX_neg(Q);
1311 54452 : return mkvec2(P,Q);
1312 : }
1313 :
1314 : static GEN
1315 108890 : hyperell_redQ(GEN W)
1316 : {
1317 108890 : GEN P = gel(W,1), Q = gel(W,2);
1318 108890 : GEN Pr, Qr = FpX_red(Q, gen_2);
1319 108890 : Pr = ZX_add(P, ZX_shifti(ZX_mul(ZX_sub(Q, Qr),ZX_add(Q, Qr)),-2));
1320 108890 : return mkvec2(Pr, Qr);
1321 : }
1322 :
1323 : static GEN
1324 50735 : hyperellisom_finalize(GEN W1, GEN W2, GEN e, GEN M, long g, long v)
1325 : {
1326 50735 : GEN Q1 = gel(W1,2), Q2 = gel(W2,2);
1327 50735 : GEN A = deg1pol_shallow(gcoeff(M,1,1), gcoeff(M,1,2), v);
1328 50735 : GEN B = deg1pol_shallow(gcoeff(M,2,1), gcoeff(M,2,2), v);
1329 50735 : GEN Bp = gpowers(B, g+1);
1330 50735 : GEN Hp = RgX_homogenous_evalpow(Q1, A, Bp, g+1);
1331 50735 : GEN H = RgX_mul2n(RgX_sub(RgX_Rg_mul(Q2,e), Hp),-1);
1332 50735 : return mkvec3(e, M, H);
1333 : }
1334 :
1335 : static void
1336 54494 : check_hyperell_Q(const char *fun, GEN *pW, GEN *pF)
1337 : {
1338 54494 : GEN W = *pW, F = check_hyperell(W);
1339 : long v, g;
1340 54494 : if (!F || !signe(F) || !RgX_is_ZX(F)) pari_err_TYPE(fun, W);
1341 54487 : if (!signe(ZX_disc(F))) pari_err_DOMAIN(fun,"disc(W)","==",gen_0,W);
1342 54480 : v = varn(F); g = hyperellgenus(F);
1343 54480 : if (g == 0) pari_err_DOMAIN(fun, "genus", "=", gen_0, gen_0);
1344 54466 : if (typ(W)==t_POL) W = mkvec2(W, pol_0(v));
1345 : else
1346 : {
1347 45010 : GEN P = gel(W, 1), Q = gel(W, 2);
1348 45010 : if (typ(P)!=t_POL) P = scalarpol_shallow(P, v);
1349 45010 : if (typ(Q)!=t_POL) Q = scalarpol_shallow(Q, v);
1350 45010 : if (!RgX_is_ZX(P) || !RgX_is_ZX(Q)) pari_err_TYPE(fun,W);
1351 45010 : if (degpol(P) > 2*g+2) pari_err_DOMAIN(fun, "deg(P)", ">", utoi(2*g+2), P);
1352 45010 : if (degpol(Q) > g+1) pari_err_DOMAIN(fun, "deg(Q)", ">", utoi(g+1), Q);
1353 45010 : W = mkvec2(P, Q);
1354 : }
1355 54466 : *pW = W; *pF = F;
1356 54466 : }
1357 :
1358 : GEN
1359 54452 : hyperellminimalmodel(GEN W, GEN *pM, GEN pr)
1360 : {
1361 54452 : pari_sp av = avma;
1362 : GEN Wr, F, WM2, F2, W2, M2, Modd, Wf, ef, Mf;
1363 : long g, v;
1364 54452 : check_hyperell_Q("hyperellminimalmodel",&W, &F);
1365 54452 : if (pr && (!is_vec_t(typ(pr)) || !RgV_is_ZV(pr)))
1366 14 : pari_err_TYPE("hyperellminimalmodel",pr);
1367 54438 : g = hyperellgenus(F); v = varn(F);
1368 54438 : Wr = hyperell_redQ(W);
1369 54438 : if (!pr || RgV_isin(pr, gen_2))
1370 : {
1371 52205 : WM2 = algo56(Wr,g); W2 = gel(WM2, 1); M2 = gel(WM2, 2);
1372 52205 : F2 = check_hyperell(W2);
1373 : }
1374 : else
1375 : {
1376 2233 : W2 = Wr; F2 = F; M2 = mkvec2(gen_1, matid(2));
1377 : }
1378 54438 : Modd = gel(algo57(F2, g, pr), 2);
1379 54438 : Wf = hyperell_redQ(minimalmodel_merge(W2, Modd, g, v));
1380 54438 : if (!pM) return gc_GEN(av, Wf);
1381 50721 : ef = mulii(gel(M2,1), gel(Modd,1));
1382 50721 : Mf = ZM2_mul(gel(M2,2), gel(Modd,2));
1383 50721 : *pM = hyperellisom_finalize(W, Wf, ef, Mf, g, v);
1384 50721 : return gc_all(av, 2, &Wf, pM);
1385 : }
1386 :
1387 : GEN
1388 14 : hyperellminimaldisc(GEN W, GEN pr)
1389 : {
1390 14 : pari_sp av = avma;
1391 14 : GEN C = hyperellminimalmodel(W, NULL, pr);
1392 14 : return gc_INT(av, hyperelldisc(C));
1393 : }
1394 :
1395 : static GEN
1396 35 : redqfbsplit(GEN a, GEN b, GEN c, GEN d)
1397 : {
1398 35 : GEN p = subii(d,b), q = shifti(a,1);
1399 35 : GEN U, Q, u, v, w = bezout(p, q, &u, &v);
1400 :
1401 35 : if (!equali1(w)) { p = diviiexact(p, w); q = diviiexact(q, w); }
1402 35 : U = mkmat22(p, negi(v), q, u);
1403 35 : Q = qfb3_SL2_apply(mkvec3(a,b,c), U);
1404 35 : b = gel(Q, 2); c = gel(Q,3);
1405 35 : if (signe(b) < 0) gel(U,2) = mkcol2(v, negi(u));
1406 35 : gel(U,2) = ZC_lincomb(gen_1, truedivii(negi(c), d), gel(U,2), gel(U,1));
1407 35 : return U;
1408 : }
1409 :
1410 : static GEN
1411 16386 : polreduce(GEN P, GEN M)
1412 : {
1413 16386 : long v = varn(P), dP = degpol(P), d = odd(dP) ? dP+1: dP;
1414 16386 : GEN A = deg1pol_shallow(gcoeff(M,1,1), gcoeff(M,1,2), v);
1415 16386 : GEN B = deg1pol_shallow(gcoeff(M,2,1), gcoeff(M,2,2), v);
1416 16386 : return RgX_homogenous_evalpow(P, A, gpowers(B, d), d);
1417 : }
1418 :
1419 : /* assume deg(P) > 2 */
1420 : static GEN
1421 8193 : red_Cremona_Stoll(GEN P, GEN *pM)
1422 : {
1423 : GEN q1, q2, q3, M, R;
1424 8193 : long i, prec = nbits2prec(2*gexpo(P)) + EXTRAPRECWORD, d = degpol(P);
1425 8193 : GEN dP = ZX_deriv(P);
1426 : for (;;)
1427 0 : {
1428 8193 : GEN r = QX_complex_roots(P, prec);
1429 8193 : q1 = gen_0; q2 = gen_0; q3 = gen_0;
1430 41000 : for (i = 1; i <= d; i++)
1431 : {
1432 32807 : GEN ri = gel(r,i);
1433 32807 : GEN s = ginv(gabs(RgX_cxeval(dP,ri,NULL), prec));
1434 32807 : if (d!=4) s = gpow(s, gdivgs(gen_2,d-2), prec);
1435 32807 : q1 = gadd(q1, s);
1436 32807 : q2 = gsub(q2, gmul(real_i(ri), s));
1437 32807 : q3 = gadd(q3, gmul(gnorm(ri), s));
1438 : }
1439 8193 : M = lllgram(mkmat22(q1,q2,q2,q3));
1440 8193 : if (M && lg(M) == 3) break;
1441 0 : prec = precdbl(prec);
1442 : }
1443 8193 : R = polreduce(P, M);
1444 8193 : *pM = M;
1445 8193 : return R;
1446 : }
1447 :
1448 : /* assume deg(P) > 2 */
1449 : GEN
1450 8193 : ZX_hyperellred(GEN P, GEN *pM)
1451 : {
1452 8193 : pari_sp av = avma;
1453 8193 : long d = degpol(P);
1454 : GEN q1, q2, q3, D, vD;
1455 8193 : GEN a = gel(P,d+2), b = gel(P,d+1), c = gel(P, d);
1456 : GEN M, R, M2;
1457 :
1458 8193 : q1 = muliu(sqri(a), d);
1459 8193 : q2 = shifti(mulii(a,b), 1);
1460 8193 : q3 = subii(sqri(b), shifti(mulii(a,c), 1));
1461 8193 : D = gcdii(gcdii(q1, q2), q3);
1462 8193 : if (!equali1(D))
1463 : {
1464 8172 : q1 = diviiexact(q1, D);
1465 8172 : q2 = diviiexact(q2, D);
1466 8172 : q3 = diviiexact(q3, D);
1467 : }
1468 8193 : D = qfb_disc3(q1, q2, q3);
1469 8193 : if (!signe(D))
1470 49 : M = mkmat22(gen_1, truedivii(negi(q2),shifti(q1,1)), gen_0, gen_1);
1471 8144 : else if (issquareall(D,&vD))
1472 35 : M = redqfbsplit(q1, q2, q3, vD);
1473 : else
1474 8109 : M = gel(qfbredsl2(mkqfb(q1,q2,q3,D), NULL), 2);
1475 8193 : R = red_Cremona_Stoll(polreduce(P, M), &M2);
1476 8193 : if (pM) *pM = gmul(M, M2);
1477 8193 : return gc_all(av, pM ? 2: 1, &R, pM);
1478 : }
1479 :
1480 : GEN
1481 42 : hyperellred(GEN W, GEN *pM)
1482 : {
1483 42 : pari_sp av = avma;
1484 : long g, v;
1485 : GEN F, M, Wf;
1486 42 : check_hyperell_Q("hyperellred", &W, &F);
1487 14 : g = hyperellgenus(F); v = varn(F);
1488 14 : (void) ZX_hyperellred(F, &M);
1489 14 : Wf = hyperell_redQ(minimalmodel_merge(W, mkvec2(gen_1, M), g, v));
1490 14 : if (pM) *pM = hyperellisom_finalize(W, Wf, gen_1, M, g, v);
1491 14 : return gc_all(av, pM ? 2: 1, &Wf, pM);
1492 : }
1493 :
1494 : static void
1495 152859 : check_hyperell_vc(const char *fun, GEN C, long v, GEN *e, GEN *M, GEN *H)
1496 : {
1497 152859 : if (typ(C) != t_VEC || lg(C) != 4) pari_err_TYPE(fun,C);
1498 152852 : *e = gel(C,1); *M = gel(C,2); *H = gel(C,3);
1499 152852 : if (typ(*M) != t_MAT || lg(*M) != 3 || lgcols(*M) != 3) pari_err_TYPE(fun,C);
1500 152845 : if (typ(*H) != t_POL || varncmp(varn(*H),v) > 0) *H = scalarpol_shallow(*H,v);
1501 152845 : if (varncmp(gvar(*M),v) <= 0) pari_err_PRIORITY(fun,*M,"<=",v);
1502 152845 : }
1503 :
1504 : GEN
1505 108857 : hyperellchangecurve(GEN W, GEN C)
1506 : {
1507 108857 : pari_sp av = avma;
1508 : GEN F, P, Q, A, B, Bp, e, M, H;
1509 : long g, v;
1510 :
1511 108857 : check_hyperell_Rg("hyperellchangecurve",&W,&F);
1512 108843 : P = gel(W,1); Q = gel(W,2);
1513 108843 : g = hyperellgenus(F); v = varn(F);
1514 108843 : check_hyperell_vc("hyperellchangecurve", C, v, &e, &M, &H);
1515 108829 : A = deg1pol_shallow(gcoeff(M,1,1), gcoeff(M,1,2), v);
1516 108829 : B = deg1pol_shallow(gcoeff(M,2,1), gcoeff(M,2,2), v);
1517 108829 : Bp = gpowers(B, 2*g+2);
1518 108829 : P = RgX_homogenous_evalpow(P, A, Bp, 2*g+2);
1519 108829 : Q = RgX_homogenous_evalpow(Q, A, Bp, g+1);
1520 108829 : P = RgX_Rg_div(RgX_sub(P, RgX_mul(H,RgX_add(Q,H))), gsqr(e));
1521 108829 : Q = RgX_Rg_div(RgX_add(Q, RgX_mul2n(H,1)), e);
1522 108829 : return gc_GEN(av, mkvec2(P,Q));
1523 : }
1524 :
1525 : static int
1526 413 : checkhyperellpt_i(GEN pt, GEN *x, GEN *y, GEN *z)
1527 : {
1528 413 : if (typ(pt) != t_VEC || lg(pt)<2 || lg(pt)>4)
1529 0 : { *x=NULL; *y=NULL; *z=NULL; return 0; }
1530 413 : if (lg(pt) == 2)
1531 : {
1532 0 : *x = gen_1; *y = gel(pt,1); *z = gen_0;
1533 : } else
1534 : {
1535 413 : *x = gel(pt,1); *y = gel(pt,2);
1536 413 : *z = lg(pt)==3 ? gen_1: gel(pt, 3);
1537 : }
1538 413 : return 1;
1539 : }
1540 :
1541 : static GEN
1542 329 : wprojtoaff(GEN X, GEN Y, GEN Z, GEN pt, long g)
1543 : {
1544 329 : if (lg(pt)==4) return mkvec3(X,Y,Z);
1545 329 : return gequal0(Z) ? mkvec(gequal0(Y) ? gen_0: gdiv(Y,gpowgs(X,g+1)))
1546 329 : : mkvec2(gdiv(X,Z),gequal0(Y) ? gen_0: gdiv(Y, gpowgs(Z,g+1)));
1547 : }
1548 :
1549 : GEN
1550 70 : hyperellchangepointinv(GEN W, GEN pt, GEN C)
1551 : {
1552 70 : pari_sp av = avma;
1553 : GEN F, e, M, H, x, y, z, X, Y,Z;
1554 : long g, v;
1555 :
1556 70 : check_hyperell_Rg("hyperellchangepointinv",&W,&F);
1557 70 : g = hyperellgenus(F); v = varn(F);
1558 70 : check_hyperell_vc("hyperellchangepointinv", C, v, &e, &M, &H);
1559 70 : if (!checkhyperellpt_i(pt,&x,&y,&z))
1560 0 : pari_err_TYPE("hyperellchangepointinv",pt);
1561 70 : X = gadd(gmul(gcoeff(M,1,1), x), gmul(gcoeff(M,1,2),z));
1562 70 : Z = gadd(gmul(gcoeff(M,2,1), x), gmul(gcoeff(M,2,2),z));
1563 70 : Y = gadd(gmul(e, y), RgX_homogenous_eval(H, x, z, g+1));
1564 70 : return gc_GEN(av, wprojtoaff(X,Y,Z,pt,g));
1565 : }
1566 :
1567 : GEN
1568 259 : hyperellchangepoint(GEN W, GEN pt, GEN C)
1569 : {
1570 259 : pari_sp av = avma;
1571 : GEN F, e, M, H, x, y, z, X, Y, Z;
1572 : GEN a, b, c, d, D;
1573 : long g, v;
1574 :
1575 259 : check_hyperell_Rg("hyperellchangepoint",&W,&F);
1576 259 : g = hyperellgenus(F); v = varn(F);
1577 259 : check_hyperell_vc("hyperellchangepoint", C, v, &e, &M, &H);
1578 259 : if (!checkhyperellpt_i(pt,&x,&y,&z))
1579 0 : pari_err_TYPE("hyperellchangepoint",pt);
1580 259 : a = gcoeff(M,1,1); b = gcoeff(M,1,2);
1581 259 : c = gcoeff(M,2,1); d = gcoeff(M,2,2);
1582 259 : Z = gsub(gmul(a, z), gmul(c, x));
1583 259 : X = gsub(gmul(d, x), gmul(b, z));
1584 259 : D = gsub(gmul(a, d), gmul(b, c));
1585 259 : Y = gdiv(gsub(gmul(y, gpowgs(D, g+1)), RgX_homogenous_eval(H,X,Z,g+1)), e);
1586 259 : return gc_GEN(av, wprojtoaff(X,Y,Z,pt,g));
1587 : }
1588 :
1589 : GEN
1590 43561 : hyperellchangeinvert(GEN W, GEN C)
1591 : {
1592 43561 : pari_sp av = avma;
1593 : GEN F, e, M, H, ei, Mi, Hi, X, Z, Zp;
1594 : long g, v;
1595 43561 : check_hyperell_Rg("hyperellchangeinvert",&W,&F);
1596 43561 : g = hyperellgenus(F); v = varn(F);
1597 43561 : check_hyperell_vc("hyperellchangeinvert", C, v, &e, &M, &H);
1598 43561 : ei = ginv(e);
1599 43561 : Mi = RgM_inv(M);
1600 43561 : X = deg1pol_shallow(gcoeff(Mi,1,1), gcoeff(Mi,1,2), v);
1601 43561 : Z = deg1pol_shallow(gcoeff(Mi,2,1), gcoeff(Mi,2,2), v);
1602 43561 : Zp = gpowers(Z, g+1);
1603 43561 : Hi = gmul(ei, gneg(RgX_homogenous_evalpow(H, X, Zp, g+1)));
1604 43561 : return gc_GEN(av, mkvec3(ei, Mi, Hi));
1605 : }
1606 :
1607 : GEN
1608 63 : hyperellchangecompose(GEN W, GEN C1, GEN C2)
1609 : {
1610 63 : pari_sp av = avma;
1611 : GEN F, e1, M1, H1, e2, M2, H2, H, X, Z, Zp;
1612 : long g, v;
1613 63 : check_hyperell_Rg("hyperellchangecompose",&W,&F);
1614 63 : g = hyperellgenus(F); v = varn(F);
1615 63 : check_hyperell_vc("hyperellchangecompose", C1, v, &e1, &M1, &H1);
1616 63 : check_hyperell_vc("hyperellchangecompose", C2, v, &e2, &M2, &H2);
1617 63 : X = deg1pol_shallow(gcoeff(M2,1,1), gcoeff(M2,1,2), v);
1618 63 : Z = deg1pol_shallow(gcoeff(M2,2,1), gcoeff(M2,2,2), v);
1619 63 : Zp = gpowers(Z, g+1);
1620 63 : H = gadd(gmul(e1,H2),RgX_homogenous_evalpow(H1, X, Zp, g+1));
1621 63 : return gc_GEN(av, mkvec3(gmul(e1,e2), gmul(M1, M2), H));
1622 : }
1623 :
1624 : int
1625 84 : hyperellisoncurve(GEN W, GEN P)
1626 : {
1627 84 : pari_sp av = avma;
1628 : GEN x, y, z, F;
1629 : long g;
1630 : int res;
1631 84 : check_hyperell_Rg("hyperellisoncurve",&W,&F);
1632 84 : g = hyperellgenus(F);
1633 84 : if (!checkhyperellpt_i(P,&x,&y,&z)) pari_err_TYPE("hyperellisoncurve",P);
1634 84 : if (typ(W)==t_POL)
1635 0 : res = gequal(gsqr(y), RgX_homogenous_eval(W,x,z,2*g+2));
1636 : else
1637 : {
1638 : GEN zp;
1639 84 : if (typ(W)!=t_VEC || lg(W)!=3) pari_err_TYPE("hyperellisoncurve",W);
1640 84 : zp = gpowers(z, 2*g+2);
1641 84 : res = gequal(gmul(y, gadd(y,RgX_homogenous_evalpow(gel(W,2), x,zp,g+1))),
1642 84 : RgX_homogenous_evalpow(gel(W,1),x,zp,2*(g+1)));
1643 : }
1644 84 : return gc_int(av, res);
1645 : }
1646 :
1647 : /****************************************************************************/
1648 : /*** ***/
1649 : /*** genus2charpoly ***/
1650 : /*** ***/
1651 : /****************************************************************************/
1652 :
1653 : /* Half stable reduction */
1654 :
1655 : static long
1656 588 : Zst_val(GEN P, GEN f, GEN p, long vt, GEN *pR)
1657 : {
1658 588 : pari_sp av = avma;
1659 588 : long v = varn(P);
1660 : while(1)
1661 1260 : {
1662 1848 : long i, j, dm = LONG_MAX;
1663 1848 : GEN Pm = NULL;
1664 1848 : long dP = degpol(P);
1665 7532 : for (i = 0; i <= minss(dP, dm); i++)
1666 : {
1667 5684 : GEN Py = gel(P, i+2);
1668 5684 : if (signe(Py))
1669 : {
1670 4186 : if (typ(Py)==t_POL)
1671 : {
1672 3864 : long dPy = degpol(Py);
1673 12502 : for (j = 0; j <= minss(dPy, dm-i); j++)
1674 : {
1675 8638 : GEN c = gel(Py, j+2);
1676 8638 : if (signe(c))
1677 : {
1678 3556 : if (i+j < dm)
1679 : {
1680 1848 : dm = i+j;
1681 1848 : Pm = monomial(gen_1, dm, v);
1682 1848 : gel(Pm,dm+2) = gen_0;
1683 : }
1684 3556 : gel(Pm,i+2) = c;
1685 : }
1686 : }
1687 : } else
1688 : {
1689 322 : if (i < dm)
1690 : {
1691 77 : dm = i;
1692 77 : Pm = monomial(Py, dm, v);
1693 : }
1694 : else
1695 245 : gel(Pm, i+2) = Py;
1696 : }
1697 : }
1698 : }
1699 1848 : Pm = RgX_renormalize(Pm);
1700 1848 : if (ZX_pval(Pm,p)==0)
1701 : {
1702 588 : *pR = gc_GEN(av, P);
1703 588 : return dm;
1704 : }
1705 1260 : Pm = RgX_homogenize_deg(Pm, dm, vt);
1706 1260 : P = gadd(gsub(P, Pm), gmul(f, ZXX_Z_divexact(Pm, p)));
1707 : }
1708 : }
1709 :
1710 : static long
1711 588 : Zst_normval(GEN P, GEN f, GEN p, long vt, GEN *pR)
1712 : {
1713 588 : long v = Zst_val(P, f, p, vt, pR);
1714 588 : long e = RgX_val(*pR)>>1;
1715 588 : if (e > 0)
1716 : {
1717 0 : v -= 2*e;
1718 0 : *pR = RgX_shift(*pR, -2*e);
1719 : }
1720 588 : return v;
1721 : }
1722 :
1723 : static GEN
1724 1176 : RgXY_swapsafe(GEN P, long v1, long v2)
1725 : {
1726 1176 : if (varn(P)==v2)
1727 : {
1728 77 : P = shallowcopy(P); setvarn(P,v1); return P;
1729 : } else
1730 1099 : return RgXY_swap(P, RgXY_degreex(P), v2);
1731 : }
1732 :
1733 : static GEN
1734 588 : Zst_red1(GEN P, GEN f, GEN p, long vt)
1735 : {
1736 588 : pari_sp av = avma;
1737 : GEN r, f1, f2, P1, P2;
1738 588 : long vs = varn(P);
1739 588 : long w = Zst_normval(P, f, p, vt, &r), ww = w-odd(w);
1740 588 : GEN st = monomial(pol_x(vt), 1, vs);
1741 588 : f1 = gsubst(f, vt, st);
1742 588 : P1 = gsubst(gdiv(r, monomial(gen_1,ww,vs)),vt,st);
1743 588 : f2 = gsubst(f, vs, st);
1744 588 : P2 = gsubst(gdiv(r, monomial(gen_1,ww,vt)),vs,st);
1745 588 : f2 = RgXY_swapsafe(f2, vs, vt);
1746 588 : P2 = RgXY_swapsafe(P2, vs, vt);
1747 588 : return gc_GEN(av, mkvec4(P1, f1, P2, f2));
1748 : }
1749 :
1750 : static GEN
1751 1176 : Zst_reduce(GEN P, GEN p, long vt, long *pv)
1752 : {
1753 : GEN C;
1754 1176 : long v = RgX_val(P);
1755 1176 : *pv = v + ZXX_pvalrem(RgX_shift(P, -v), p, &P);
1756 1176 : C = constant_coeff(P);
1757 1176 : C = typ(C) == t_POL ? C: scalarpol_shallow(C, vt);
1758 1176 : return FpX_red(C, p);
1759 : }
1760 :
1761 : static GEN
1762 588 : Zst_red3(GEN C, GEN p, long vt)
1763 : {
1764 : while(1)
1765 511 : {
1766 588 : GEN P1 = gel(C,1), f1 = gel(C,2), Poo = gel(C,3), foo= gel(C,4);
1767 : long e;
1768 588 : GEN Qoop = Zst_reduce(Poo, p, vt, &e), Qp, R;
1769 588 : if (RgX_val(Qoop) >= 3-e)
1770 : {
1771 0 : C = Zst_red1(Poo, foo, p, vt);
1772 511 : continue;
1773 : }
1774 588 : Qp = Zst_reduce(P1, p, vt, &e);
1775 588 : R = FpX_roots_mult(Qp, 3-e, p);
1776 588 : if (lg(R) > 1)
1777 511 : {
1778 511 : GEN xz = deg1pol_shallow(gen_1, gel(R,1), vt);
1779 511 : C = Zst_red1(gsubst(P1, vt, xz), gsubst(f1, vt, xz), p, vt);
1780 511 : continue;
1781 : }
1782 77 : return Qp;
1783 : }
1784 : }
1785 :
1786 : static GEN
1787 77 : genus2_halfstablemodel_i(GEN P, GEN p, long vt)
1788 : {
1789 : GEN Qp, R, Poo, Qoop;
1790 77 : long e = ZX_pvalrem(P, p, &Qp);
1791 77 : R = FpX_roots_mult(FpX_red(Qp,p), 4-e, p);
1792 77 : if (lg(R) > 1)
1793 : {
1794 77 : GEN C = Zst_red1(ZX_Z_translate(P, gel(R,1)), pol_x(vt), p, vt);
1795 77 : return Zst_red3(C, p, vt);
1796 : }
1797 0 : Poo = RgXn_recip_shallow(P, 7);
1798 0 : e = ZX_pvalrem(Poo, p, &Qoop);
1799 0 : Qoop = FpX_red(Qoop,p);
1800 0 : if (RgX_val(Qoop)>=4-e)
1801 : {
1802 0 : GEN C = Zst_red1(Poo, pol_x(vt), p, vt);
1803 0 : return Zst_red3(C, p, vt);
1804 : }
1805 0 : return gcopy(P);
1806 : }
1807 :
1808 : static GEN
1809 77 : genus2_halfstablemodel(GEN P, GEN p)
1810 : {
1811 77 : pari_sp av = avma;
1812 77 : long vt = fetch_var(), vs = varn(P);
1813 77 : GEN S = genus2_halfstablemodel_i(P, p, vt);
1814 77 : setvarn(S, vs); delete_var();
1815 77 : return gc_GEN(av, S);
1816 : }
1817 :
1818 : /* semi-stable reduction */
1819 :
1820 : static GEN
1821 1015 : genus2_redmodel(GEN P, GEN p)
1822 : {
1823 : GEN LP, U, F;
1824 : long i, k, r;
1825 1015 : if (degpol(P) < 0) return mkvec2(cgetg(1, t_COL), P);
1826 980 : F = FpX_factor_squarefree(P, p);
1827 980 : r = lg(F); U = NULL;
1828 3416 : for (i = k = 1; i < r; i++)
1829 : {
1830 2436 : GEN f = gel(F,i);
1831 2436 : long df = degpol(f);
1832 2436 : if (!df) continue;
1833 1687 : if (odd(i)) U = U? FpX_mul(U, f, p): f;
1834 1687 : if (i > 1) gel(F,k++) = df == 1? mkcol(f): gel(FpX_factor(f, p), 1);
1835 : }
1836 980 : LP = leading_coeff(P);
1837 980 : if (!U)
1838 154 : U = scalarpol_shallow(LP, varn(P));
1839 : else
1840 : {
1841 826 : GEN LU = leading_coeff(U);
1842 826 : if (!equalii(LU, LP)) U = FpX_Fp_mul(U, Fp_div(LP, LU, p), p);
1843 : }
1844 980 : setlg(F,k); if (k > 1) F = shallowconcat1(F);
1845 980 : return mkvec2(F, U);
1846 : }
1847 :
1848 : static GEN
1849 5852 : xdminusone(long d)
1850 : {
1851 5852 : return gsub(pol_xn(d, 0),gen_1);
1852 : }
1853 :
1854 : static GEN
1855 623 : ellfromeqncharpoly(GEN P, GEN Q, GEN p)
1856 : {
1857 : long v;
1858 : GEN E, F, t, y;
1859 623 : v = fetch_var();
1860 623 : y = pol_x(v);
1861 623 : F = gsub(gadd(ZX_sqr(y), gmul(y, Q)), P);
1862 623 : E = ellinit(ellfromeqn(F), p, DEFAULTPREC);
1863 623 : delete_var();
1864 623 : t = ellcharpoly(E, p);
1865 623 : obj_free(E);
1866 623 : return t;
1867 : }
1868 :
1869 : static GEN
1870 0 : nfellcharpoly(GEN e, GEN T, GEN p)
1871 : {
1872 : GEN nf, E, t;
1873 0 : e = shallowcopy(e);
1874 0 : nf = nfinit(mkvec2(T, mkvec(p)), DEFAULTPREC);
1875 : while(1)
1876 : {
1877 0 : E = ellinit(e, nf, DEFAULTPREC);
1878 0 : if (lg(E)!=1) break;
1879 0 : gel(e,5) = gadd(gel(e,5), p);
1880 : }
1881 0 : t = elleulerf(E, p);
1882 0 : obj_free(E);
1883 0 : return RgX_recip(ginv(t));
1884 : }
1885 :
1886 : static GEN
1887 0 : genus2_red5(GEN P, GEN T, GEN p)
1888 : {
1889 0 : long vx = varn(P), vy = varn(T);
1890 0 : GEN f = shallowcopy(T), pi = shifti(p,-1);
1891 0 : setvarn(f, vx);
1892 : while(1)
1893 0 : {
1894 : GEN Pr, R, r, Rs;
1895 0 : long v = ZXX_pvalrem(P, p, &Pr);
1896 0 : R = FpXQX_roots_mult(Pr, 2-v, T, p);
1897 0 : if (lg(R)==1) return P;
1898 0 : r = FpX_center(gel(R,1), p, pi);
1899 0 : Pr = RgX_affine(P, p, r);
1900 0 : setvarn(r, vx);
1901 0 : f = RgX_Rg_div(gsub(f, r), p);
1902 0 : Rs = RgX_rem(RgXY_swap(Pr, 3, vy), gsub(f, pol_x(vy)));
1903 0 : Pr = RgXY_swap(Rs, 3, vy);
1904 0 : if (ZXX_pvalrem(Pr, sqri(p), &Pr)==0) return P;
1905 0 : P = Pr;
1906 : }
1907 : }
1908 :
1909 : static GEN
1910 1029 : genus2_type5(GEN P, GEN p)
1911 : {
1912 : GEN E, F, T, a, a2, Q;
1913 : long v;
1914 1029 : if (equaliu(p, 2))
1915 224 : (void) ZXX_pvalrem(P, sqri(p), &P);
1916 1029 : (void) ZX_pvalrem(P, p, &F);
1917 1029 : F = FpX_red(F, p);
1918 1029 : if (degpol(F) < 1) return NULL;
1919 1015 : F = FpX_factor(F, p);
1920 1015 : if (mael(F,2,1) != 3 || degpol(gmael(F,1,1)) != 2) return NULL;
1921 0 : T = gmael(F, 1, 1);
1922 0 : v = fetch_var_higher();
1923 0 : Q = RgV_to_RgX(ZX_digits(P, T), v);
1924 0 : Q = genus2_red5(Q, T, p);
1925 0 : a = gel(Q,5); a2 = ZX_sqr(a);
1926 0 : E = mkvec5(gen_0, gel(Q,4), gen_0, ZX_mul(gel(Q,3),a), ZX_mul(gel(Q,2),a2));
1927 0 : delete_var();
1928 0 : return nfellcharpoly(E, T, p);
1929 : }
1930 :
1931 : /* Assume P has semistable reduction at p */
1932 : static GEN
1933 1015 : genus2_eulerfact_semistable(GEN P, GEN p)
1934 : {
1935 1015 : GEN Pp = FpX_red(P, p);
1936 1015 : GEN GU = genus2_redmodel(Pp, p);
1937 1015 : long d = 6-degpol(Pp), v = d/2, w = odd(d);
1938 : GEN abe, tor;
1939 1015 : GEN ki, kp = pol_1(0), kq = pol_1(0);
1940 1015 : GEN F = gel(GU,1), Q = gel(GU,2);
1941 1015 : long dQ = degpol(Q), lF = lg(F)-1;
1942 :
1943 7 : abe = dQ >= 5 ? hyperellcharpoly(gmul(Q,gmodulo(gen_1,p)))
1944 2023 : : dQ >= 3 ? ellfromeqncharpoly(Q,gen_0,p)
1945 1008 : : pol_1(0);
1946 861 : ki = dQ != 0 ? xdminusone(1)
1947 1169 : : Fp_issquare(gel(Q,2),p) ? ZX_sqr(xdminusone(1))
1948 154 : : xdminusone(2);
1949 1015 : if (lF)
1950 : {
1951 : long i;
1952 2100 : for(i=1; i <= lF; i++)
1953 : {
1954 1183 : GEN Fi = gel(F, i);
1955 1183 : long d = degpol(Fi);
1956 1183 : GEN e = FpX_rem(Q, Fi, p);
1957 2100 : GEN kqf = lgpol(e)==0 ? xdminusone(d):
1958 1561 : FpXQ_issquare(e, Fi, p) ? ZX_sqr(xdminusone(d))
1959 917 : : xdminusone(2*d);
1960 1183 : kp = gmul(kp, xdminusone(d));
1961 1183 : kq = gmul(kq, kqf);
1962 : }
1963 : }
1964 1015 : if (v)
1965 : {
1966 273 : GEN kqoo = w==1 ? xdminusone(1):
1967 21 : Fp_issquare(leading_coeff(Q), p)? ZX_sqr(xdminusone(1))
1968 14 : : xdminusone(2);
1969 259 : kp = gmul(kp, xdminusone(1));
1970 259 : kq = gmul(kq, kqoo);
1971 : }
1972 1015 : tor = RgX_div(ZX_mul(xdminusone(1), kq), ZX_mul(ki, kp));
1973 1015 : return ZX_mul(abe, tor);
1974 : }
1975 :
1976 : GEN
1977 1428 : genus2_eulerfact(GEN P, GEN p, long ra, long rt)
1978 : {
1979 1428 : pari_sp av = avma;
1980 : GEN W, R, E;
1981 1428 : long d = 2*ra+rt;
1982 1428 : if (d == 0) return pol_1(0);
1983 805 : R = genus2_type5(P, p);
1984 805 : if (R) return R;
1985 805 : W = hyperellextremalmodels_i(P, 2, p);
1986 805 : if (lg(W) < 3)
1987 : {
1988 672 : GEN F = genus2_eulerfact_semistable(P,p);
1989 672 : if (degpol(F)!=d)
1990 : {
1991 77 : GEN S = genus2_halfstablemodel(P, p);
1992 77 : F = genus2_eulerfact_semistable(S, p);
1993 77 : if (degpol(F)!=d) pari_err_BUG("genus2charpoly");
1994 : }
1995 672 : return F;
1996 : }
1997 133 : E = gmul(genus2_eulerfact_semistable(gel(W,1),p),
1998 133 : genus2_eulerfact_semistable(gel(W,2),p));
1999 133 : return gc_upto(av, E);
2000 : }
2001 :
2002 : /* p = 2 */
2003 :
2004 : static GEN
2005 238 : F2x_genus2_find_trans(GEN P, GEN Q, GEN F)
2006 : {
2007 238 : pari_sp av = avma;
2008 238 : long i, d = F2x_degree(F), v = P[1];
2009 : GEN M, C, V;
2010 238 : M = cgetg(d+1, t_MAT);
2011 798 : for (i=1; i<=d; i++)
2012 : {
2013 560 : GEN Mi = F2x_rem(F2x_add(F2x_shift(Q,i-1), monomial_F2x(2*i-2,v)), F);
2014 560 : gel(M,i) = F2x_to_F2v(Mi, d);
2015 : }
2016 238 : C = F2x_to_F2v(F2x_rem(P, F), d);
2017 238 : V = F2m_F2c_invimage(M, C);
2018 238 : return gc_leaf(av, F2v_to_F2x(V, v));
2019 : }
2020 :
2021 : static GEN
2022 280 : F2x_genus2_trans(GEN P, GEN Q, GEN H)
2023 : {
2024 280 : return F2x_add(P,F2x_add(F2x_mul(H,Q), F2x_sqr(H)));
2025 : }
2026 :
2027 : static GEN
2028 693 : F2x_genus_redoo(GEN P, GEN Q, long k)
2029 : {
2030 693 : if (F2x_degree(P)==2*k)
2031 : {
2032 126 : long c = F2x_coeff(P,2*k-1), dQ = F2x_degree(Q);
2033 126 : if ((dQ==k-1 && c==1) || (dQ<k-1 && c==0))
2034 42 : return F2x_genus2_trans(P, Q, monomial_F2x(k, P[1]));
2035 : }
2036 651 : return P;
2037 : }
2038 :
2039 : static GEN
2040 434 : F2x_pseudodisc(GEN P, GEN Q)
2041 : {
2042 434 : GEN dP = F2x_deriv(P), dQ = F2x_deriv(Q);
2043 434 : return F2x_gcd(Q, F2x_add(F2x_mul(P, F2x_sqr(dQ)), F2x_sqr(dP)));
2044 : }
2045 :
2046 : static GEN
2047 231 : F2x_genus_red(GEN P, GEN Q)
2048 : {
2049 : long dP, dQ;
2050 : GEN F, FF;
2051 231 : P = F2x_genus_redoo(P, Q, 3);
2052 231 : P = F2x_genus_redoo(P, Q, 2);
2053 231 : P = F2x_genus_redoo(P, Q, 1);
2054 231 : dP = F2x_degree(P);
2055 231 : dQ = F2x_degree(Q);
2056 231 : FF = F = F2x_pseudodisc(P,Q);
2057 434 : while(F2x_degree(F)>0)
2058 : {
2059 203 : GEN M = gel(F2x_factor(F),1);
2060 203 : long i, l = lg(M);
2061 441 : for(i=1; i<l; i++)
2062 : {
2063 238 : GEN R = F2x_sqr(gel(M,i));
2064 238 : GEN H = F2x_genus2_find_trans(P, Q, R);
2065 238 : P = F2x_div(F2x_genus2_trans(P, Q, H), R);
2066 238 : Q = F2x_div(Q, gel(M,i));
2067 : }
2068 203 : F = F2x_pseudodisc(P, Q);
2069 : }
2070 231 : return mkvec4(P,Q,FF,mkvecsmall2(dP,dQ));
2071 : }
2072 :
2073 : /* Number of solutions of x^2+b*x+c */
2074 : static long
2075 210 : F2xqX_quad_nbroots(GEN b, GEN c, GEN T)
2076 : {
2077 210 : if (lgpol(b) > 0)
2078 : {
2079 147 : GEN d = F2xq_div(c, F2xq_sqr(b, T), T);
2080 147 : return F2xq_trace(d, T)? 0: 2;
2081 : }
2082 : else
2083 63 : return 1;
2084 : }
2085 :
2086 : static GEN
2087 231 : genus2_eulerfact2_semistable(GEN PQ)
2088 : {
2089 231 : GEN V = F2x_genus_red(ZX_to_F2x(gel(PQ, 1)), ZX_to_F2x(gel(PQ, 2)));
2090 231 : GEN P = gel(V, 1), Q = gel(V, 2);
2091 231 : GEN F = gel(V, 3), v = gel(V, 4);
2092 : GEN abe, tor;
2093 231 : GEN ki, kp = pol_1(0), kq = pol_1(0);
2094 231 : long dP = F2x_degree(P), dQ = F2x_degree(Q), d = maxss(dP, 2*dQ);
2095 231 : if (!lgpol(F)) return pol_1(0);
2096 259 : ki = dQ!=0 || dP>0 ? xdminusone(1):
2097 49 : dP==-1 ? ZX_sqr(xdminusone(1)): xdminusone(2);
2098 420 : abe = d>=5? hyperellcharpoly(gmul(PQ,gmodulss(1,2))):
2099 210 : d>=3? ellfromeqncharpoly(F2x_to_ZX(P), F2x_to_ZX(Q), gen_2):
2100 70 : pol_1(0);
2101 210 : if (lgpol(F))
2102 : {
2103 210 : GEN M = gel(F2x_factor(F), 1);
2104 210 : long i, lF = lg(M)-1;
2105 420 : for(i=1; i <= lF; i++)
2106 : {
2107 210 : GEN Fi = gel(M, i);
2108 210 : long d = F2x_degree(Fi);
2109 210 : long nb = F2xqX_quad_nbroots(F2x_rem(Q, Fi), F2x_rem(P, Fi), Fi);
2110 357 : GEN kqf = nb==1 ? xdminusone(d):
2111 63 : nb==2 ? ZX_sqr(xdminusone(d))
2112 147 : : xdminusone(2*d);
2113 210 : kp = gmul(kp, xdminusone(d));
2114 210 : kq = gmul(kq, kqf);
2115 : }
2116 : }
2117 210 : if (maxss(v[1],2*v[2])<5)
2118 : {
2119 70 : GEN kqoo = v[1]>2*v[2] ? xdminusone(1):
2120 7 : v[1]<2*v[2] ? ZX_sqr(xdminusone(1))
2121 21 : : xdminusone(2);
2122 49 : kp = gmul(kp, xdminusone(1));
2123 49 : kq = gmul(kq, kqoo);
2124 : }
2125 210 : tor = RgX_div(ZX_mul(xdminusone(1),kq), ZX_mul(ki, kp));
2126 210 : return ZX_mul(abe, tor);
2127 : }
2128 :
2129 : GEN
2130 224 : genus2_eulerfact2(GEN F, GEN PQ)
2131 : {
2132 224 : pari_sp av = avma;
2133 224 : GEN W, R = genus2_type5(F, gen_2), E;
2134 224 : if (R) return R;
2135 224 : W = hyperellextremalmodels_i(PQ, 2, gen_2);
2136 224 : if (lg(W) < 3) return genus2_eulerfact2_semistable(PQ);
2137 7 : E = gmul(genus2_eulerfact2_semistable(gel(W,1)),
2138 7 : genus2_eulerfact2_semistable(gel(W,2)));
2139 7 : return gc_upto(av, E);
2140 : }
2141 :
2142 : GEN
2143 889 : genus2charpoly(GEN G, GEN p)
2144 : {
2145 889 : pari_sp av = avma;
2146 889 : GEN gr = genus2red(G, p), F;
2147 889 : GEN PQ = gel(gr, 3), L = gel(gr, 4), r = gel(L, 4);
2148 889 : GEN P = gadd(gsqr(gel(PQ, 2)), gmul2n(gel(PQ, 1), 2));
2149 889 : if (equaliu(p,2))
2150 7 : F = genus2_eulerfact2(P, PQ);
2151 : else
2152 882 : F = genus2_eulerfact(P,p, r[1],r[2]);
2153 889 : return gc_upto(av, F);
2154 : }
|
