Code coverage tests

This page documents the degree to which the PARI/GP source code is tested by our public test suite, distributed with the source distribution in directory src/test/. This is measured by the gcov utility; we then process gcov output using the lcov frond-end.

We test a few variants depending on Configure flags on the pari.math.u-bordeaux.fr machine (x86_64 architecture), and agregate them in the final report:

The target is to exceed 90% coverage for all mathematical modules (given that branches depending on DEBUGLEVEL or DEBUGMEM are not covered). This script is run to produce the results below.

LCOV - code coverage report
Current view: top level - basemath - hyperell.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.18.1 lcov report (development 30735-65f72320ad) Lines: 1169 1251 93.4 %
Date: 2026-03-09 09:24:55 Functions: 104 107 97.2 %
Legend: Lines: hit not hit

          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          84 : F2x_genus2charpoly_naive(GEN P, GEN Q)
     626             : {
     627          84 :   long a, b = 1, c = 0;
     628          84 :   GEN T = mkvecsmall2(P[1], 7);
     629          84 :   GEN PT = F2x_rem(P, T), QT = F2x_rem(Q, T);
     630          84 :   long q0 = F2x_eval(Q, 0), q1 = F2x_eval(Q, 1);
     631          84 :   long dP = F2x_degree(P), dQ = F2x_degree(Q);
     632          84 :   a= dQ<3 ? 0: dP<=5 ? 1: -1;
     633          84 :   a += (q0? F2x_eval(P, 0)? -1: 1: 0) + (q1? F2x_eval(P, 1)? -1: 1: 0);
     634          84 :   b += q0 + q1;
     635          84 :   if (lgpol(QT))
     636          70 :     c = (F2xq_trace(F2xq_div(PT, F2xq_sqr(QT, T), T), T)==0 ? 1: -1);
     637          84 :   return mkvecsmalln(6, 0UL, 4UL, 2*a, (b+2*c+a*a)>>1, a, 1UL);
     638             : }
     639             : 
     640             : static GEN
     641         273 : Flx_difftable(GEN P, ulong p)
     642             : {
     643         273 :   long i, n = degpol(P);
     644         273 :   GEN V = cgetg(n+2, t_VEC);
     645         273 :   gel(V, n+1) = P;
     646        1869 :   for(i = n; i >= 1; i--)
     647        1596 :     gel(V, i) = Flx_diff1(gel(V, i+1), p);
     648         273 :   return V;
     649             : }
     650             : 
     651             : static GEN
     652        1617 : FlxV_Fl2_eval_pre(GEN V, GEN x, ulong D, ulong p, ulong pi)
     653             : {
     654        1617 :   long i, n = lg(V)-1;
     655        1617 :   GEN r = cgetg(n+1, t_VEC);
     656       12390 :   for (i = 1; i <= n; i++)
     657       10773 :     gel(r, i) = Flx_Fl2_eval_pre(gel(V, i), x, D, p, pi);
     658        1617 :   return r;
     659             : }
     660             : 
     661             : static GEN
     662       44898 : Fl2V_next(GEN V, ulong p)
     663             : {
     664       44898 :   long i, n = lg(V)-1;
     665       44898 :   GEN r = cgetg(n+1, t_VEC);
     666       44898 :   gel(r, 1) = gel(V, 1);
     667      288330 :   for (i = 2; i <= n; i++)
     668      243432 :     gel(r, i) = Flv_add(gel(V, i), gel(V, i-1), p);
     669       44898 :   return r;
     670             : }
     671             : 
     672             : static GEN
     673         273 : FlxV_constant(GEN x)
     674        2142 : { pari_APPLY_long(Flx_constant(gel(x,i))) }
     675             : 
     676             : static GEN
     677         273 : Flx_genus2charpoly_naive(GEN H, ulong p)
     678             : {
     679         273 :   pari_sp av = avma, av2;
     680         273 :   ulong pi = get_Fl_red(p);
     681         273 :   ulong i, j, p2 = p>>1, D = 2, e = ((p&2UL) == 0) ? -1 : 1;
     682         273 :   long a, b, c = 0, n = degpol(H);
     683         273 :   GEN t, d, k = const_vecsmall(p, -1);
     684         273 :   k[1] = 0;
     685        1890 :   for (i=1, j=1; i < p; i += 2, j = Fl_add(j, i, p)) k[j+1] = 1;
     686         329 :   while (k[1+D] >= 0) D++;
     687         273 :   b = n == 5 ? 0 : 1;
     688         273 :   a = b ? k[1+Flx_lead(H)]: 0;
     689         273 :   t = Flx_difftable(H, p);
     690         273 :   d = FlxV_constant(t);
     691         273 :   av2 = avma;
     692        3780 :   for (i=0; i < p; i++)
     693             :   {
     694        3507 :     ulong v = uel(d,n+1);
     695        3507 :     a += k[1+v];
     696        3507 :     b += !!v;
     697        3507 :     if (n==6)
     698        2373 :       uel(d,7) = Fl_add(uel(d,7), uel(d,6), p);
     699        3507 :     uel(d,6) = Fl_add(uel(d,6), uel(d,5), p);
     700        3507 :     uel(d,5) = Fl_add(uel(d,5), uel(d,4), p);
     701        3507 :     uel(d,4) = Fl_add(uel(d,4), uel(d,3), p);
     702        3507 :     uel(d,3) = Fl_add(uel(d,3), uel(d,2), p);
     703        3507 :     uel(d,2) = Fl_add(uel(d,2), uel(d,1), p);
     704             :   }
     705        1890 :   for (j=1; j <= p2; j++)
     706             :   {
     707        1617 :     GEN V = FlxV_Fl2_eval_pre(t, mkvecsmall2(0, j), D, p, pi);
     708        1617 :     for (i=0;; i++)
     709       44898 :     {
     710       46515 :       GEN r2 = gel(V, n+1);
     711       93030 :       c += uel(r2,2) ?
     712       44135 :         (uel(r2,1) ? uel(k,1+Fl2_norm_pre(r2, D, p, pi)): e)
     713       90650 :          : !!uel(r2,1);
     714       46515 :       if (i == p-1) break;
     715       44898 :       V = Fl2V_next(V, p);
     716             :     }
     717        1617 :     set_avma(av2);
     718             :   }
     719         273 :   set_avma(av);
     720         273 :   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      289545 : check_hyperell(GEN PQ)
     745             : {
     746             :   GEN H;
     747      289545 :   if (is_vec_t(typ(PQ)) && lg(PQ)==3)
     748      225686 :     H = gadd(gsqr(gel(PQ, 2)), gmul2n(gel(PQ, 1), 2));
     749             :   else
     750       63859 :     H = gmul2n(PQ, 2);
     751      289545 :   return typ(H) == t_POL? H: NULL;
     752             : }
     753             : 
     754             : GEN
     755        1050 : hyperellcharpoly(GEN PQ)
     756             : {
     757        1050 :   pari_sp av = avma;
     758        1050 :   GEN M, R, T=NULL, pp=NULL, q;
     759        1050 :   long d, n, eps = 0;
     760             :   ulong p;
     761        1050 :   GEN H = check_hyperell(PQ);
     762        1050 :   if (!H || !RgX_is_FpXQX(H, &T, &pp) || !pp)
     763           0 :     pari_err_TYPE("hyperellcharpoly", PQ);
     764        1050 :   p = itou(pp);
     765        1050 :   if (!T)
     766             :   {
     767         903 :     if (p==2 && is_vec_t(typ(PQ)))
     768             :     {
     769          84 :       long dP, dQ, v = varn(H);
     770          84 :       GEN P = gel(PQ,1), Q = gel(PQ,2);
     771          84 :       if (typ(P)!=t_POL)  P = scalarpol(P, v);
     772          84 :       if (typ(Q)!=t_POL)  Q = scalarpol(Q, v);
     773          84 :       dP = degpol(P); dQ = degpol(Q);
     774          84 :       if (dP<=6 && dQ <=3 && (dQ==3 || dP>=5))
     775             :       {
     776          84 :         GEN P2 = RgX_to_F2x(P), Q2 = RgX_to_F2x(Q);
     777          84 :         GEN D = F2x_add(F2x_mul(P2, F2x_sqr(F2x_deriv(Q2))), F2x_sqr(F2x_deriv(P2)));
     778          84 :         if (F2x_degree(F2x_gcd(D, Q2))) is_sing(PQ, 2);
     779          84 :         if (dP==6 && dQ<3 && F2x_coeff(P2,5)==F2x_coeff(Q2,2))
     780           0 :           is_sing(PQ, 2); /* The curve is singular at infinity */
     781          84 :         R = zx_to_ZX(F2x_genus2charpoly_naive(P2, Q2));
     782          84 :         return gc_upto(av, R);
     783             :       }
     784             :     }
     785         819 :     H = RgX_to_FpX(H, pp);
     786         819 :     d = degpol(H);
     787         819 :     if (d <= 0) is_sing(H, p);
     788         819 :     if (p > 2 && ((d == 5 && p < 17500) || (d == 6 && p < 24500)))
     789             :     {
     790         280 :       GEN Hp = ZX_to_Flx(H, p);
     791         280 :       if (!Flx_is_squarefree(Hp, p)) is_sing(H, p);
     792         273 :       R = zx_to_ZX(Flx_genus2charpoly_naive(Hp, p));
     793         273 :       return gc_upto(av, R);
     794             :     }
     795         539 :     n = hyperell_Weil_bound(pp, (d-1)>>1, pp);
     796         539 :     eps = odd(d)? 0: Fp_issquare(leading_coeff(H), pp);
     797         539 :     M = hyperellpadicfrobenius(H, p, n);
     798         539 :     R = centerlift(carberkowitz(M, 0));
     799         539 :     q = pp;
     800             :   }
     801             :   else
     802             :   {
     803             :     int fixvar;
     804         147 :     T = typ(T)==t_FFELT? FF_mod(T): RgX_to_FpX(T, pp);
     805         147 :     q = powuu(p, degpol(T));
     806         147 :     fixvar = (varncmp(varn(T),varn(H)) <= 0);
     807         147 :     if (fixvar) setvarn(T, fetch_var());
     808         147 :     H = RgX_to_FpXQX(H, T, pp);
     809         147 :     d = degpol(H);
     810         147 :     if (d <= 0) is_sing(H, p);
     811         147 :     eps = odd(d)? 0: Fq_issquare(leading_coeff(H), T, pp);
     812         147 :     n = hyperell_Weil_bound(q, (d-1)>>1, pp);
     813         147 :     M = nfhyperellpadicfrobenius(H, T, p, n);
     814         140 :     R = simplify_shallow(centerlift(liftpol_shallow(carberkowitz(M, 0))));
     815         140 :     if (fixvar) (void)delete_var();
     816             :   }
     817         679 :   if (!odd(d))
     818             :   {
     819         301 :     GEN b = get_basis(p, d) == 3 ? gen_1 : q;
     820         301 :     GEN pn = powuu(p, n);
     821         301 :     R = FpX_div_by_X_x(R, eps? b: negi(b), pn, NULL);
     822         301 :     R = FpX_center_i(R, pn, shifti(pn,-1));
     823             :   }
     824         679 :   return gc_upto(av, charpoly_funceq(R, q));
     825             : }
     826             : 
     827             : int
     828        3493 : hyperellisoncurve(GEN W, GEN P)
     829             : {
     830        3493 :   pari_sp av = avma;
     831             :   long res;
     832             :   GEN x, y;
     833        3493 :   if (typ(P)!=t_VEC || lg(P)!=3) pari_err_TYPE("hyperellisoncurve",P);
     834        3493 :   x = gel(P,1); y = gel(P,2);
     835        3493 :   if (typ(W)==t_POL)
     836           0 :     res = gequal(gsqr(y), poleval(W,x));
     837             :   else
     838             :   {
     839        3493 :     if (typ(W)!=t_VEC || lg(W)!=3) pari_err_TYPE("hyperellisoncurve",W);
     840        3493 :     res = gequal(gmul(y, gadd(y,poleval(gel(W,2), x))), poleval(gel(W,1), x));
     841             :   }
     842        3493 :   return gc_int(av, res);
     843             : }
     844             : 
     845             : GEN
     846          35 : hyperellordinate(GEN W, GEN x)
     847             : {
     848          35 :   pari_sp av = avma;
     849          35 :   if (typ(W)==t_POL)
     850             :   {
     851          14 :     GEN d = poleval(W,x), y;
     852          14 :     if (gequal0(d)) { return gc_GEN(av, mkvec(d)); }
     853          14 :     if (!issquareall(d, &y)) retgc_const(av, cgetg(1, t_VEC));
     854           7 :     return gc_GEN(av, mkvec2(y, gneg(y)));
     855             :   }
     856             :   else
     857             :   {
     858             :     GEN b, c, d, rd, y;
     859          21 :     if (typ(W)!=t_VEC || lg(W)!=3) pari_err_TYPE("hyperellisoncurve",W);
     860          21 :     b = poleval(gel(W,2), x); c = poleval(gel(W,1), x);
     861          21 :     d = gadd(gsqr(b), gmul2n(c, 2));
     862          21 :     if (gequal0(d)) { return gc_GEN(av, mkvec(gmul2n(gneg(b),-1))); }
     863          14 :     if (!issquareall(d, &rd)) retgc_const(av, cgetg(1, t_VEC));
     864           7 :     y = gmul2n(gsub(rd, b), -1);
     865           7 :     return gc_GEN(av, mkvec2(y, gsub(y,rd)));
     866             :   }
     867             : }
     868             : 
     869             : static long
     870      355646 : hyperellgenus(GEN H)
     871      355646 : { long d = degpol(H); return ((d+1)>>1)-1; }
     872             : 
     873             : GEN
     874      118971 : hyperelldisc(GEN PQ)
     875             : {
     876      118971 :   pari_sp av = avma;
     877      118971 :   GEN D, H = check_hyperell(PQ);
     878             :   long g;
     879      118971 :   if (!H || signe(H)==0) pari_err_TYPE("hyperelldisc",PQ);
     880      118971 :   g = hyperellgenus(H);
     881      118971 :   D = gmul2n(RgX_disc(H),-4*(g+1));
     882      118971 :   if (odd(degpol(H))) D = gmul(D, gsqr(leading_coeff(H)));
     883      118971 :   return gc_upto(av, D);
     884             : }
     885             : 
     886             : static long
     887      126677 : get_ep(GEN W)
     888             : {
     889      126677 :   GEN P = gel(W,1), Q = gel(W,2);
     890      126677 :   if (signe(Q)==0) return ZX_lval(P,2);
     891       86412 :   return minss(ZX_lval(P,2), ZX_lval(Q,2));
     892             : }
     893             : 
     894             : static GEN
     895       50910 : algo51(GEN W, GEN M)
     896             : {
     897       50910 :   GEN P = gel(W,1), Q = gel(W,2);
     898             :   for(;;)
     899       10647 :   {
     900       61557 :     long vP = ZX_lval(P,2);
     901       61557 :     long vQ = signe(Q) ? ZX_lval(Q,2): vP+1;
     902             :     long r;
     903             :     /* 1 */
     904       61557 :     if (vQ==0) break;
     905             :     /* 2 */
     906       36204 :     if (vP==0)
     907             :     {
     908             :       GEN H, H1;
     909             :       /* a */
     910       29644 :       RgX_even_odd(FpX_red(P,gen_2),&H, &H1);
     911       29644 :       if (signe(H1)) break;
     912             :       /* b */
     913       14965 :       P = ZX_add(P, ZX_mul(H, ZX_sub(Q, H)));
     914       14965 :       Q = ZX_sub(Q, ZX_shifti(H, 1));
     915       14965 :       vP = ZX_lval(P,2);
     916       14965 :       vQ = signe(Q) ? ZX_lval(Q,2): vP+1;
     917             :     }
     918             :     /* 2c */
     919       21525 :     if (vP==1) break;
     920             :     /* 2d */
     921       10647 :     r = minss(2*vQ, vP)>>1;
     922       10647 :     if (M) gel(M,1) = shifti(gel(M,1), r);
     923       10647 :     P = ZX_shifti(P, -2*r);
     924       10647 :     Q = ZX_shifti(Q, -r);
     925             :   }
     926       50910 :   return mkvec2(P,Q);
     927             : }
     928             : 
     929             : static GEN
     930      103469 : algo52(GEN W, GEN c, long *pt_lambda)
     931             : {
     932             :   long lambda;
     933      103469 :   GEN P = gel(W,1), Q = gel(W,2);
     934             :   for(;;)
     935      117043 :   {
     936             :     GEN H, H1;
     937             :     /* 1 */
     938      220512 :     GEN Pc = ZX_affine(P,gen_2,c), Qc = ZX_affine(Q,gen_2,c);
     939      220512 :     long mP = ZX_lval(Pc,2), mQ = signe(Qc) ? ZX_lval(Qc,2): mP+1;
     940             :     /* 2 */
     941      220512 :     if (2*mQ <= mP) { lambda = 2*mQ; break; }
     942             :     /* 3 */
     943      188220 :     if (odd(mP)) { lambda = mP; break; }
     944             :     /* 4 */
     945      127816 :     RgX_even_odd(FpX_red(ZX_shifti(Pc, -mP),gen_2),&H, &H1);
     946      127816 :     if (signe(H1)) { lambda = mP; break; }
     947             :     /* 5 */
     948      117043 :     P = ZX_add(P, ZX_mul(H, ZX_sub(Q, H)));
     949      117043 :     Q = ZX_sub(Q, ZX_shifti(H, 1));
     950             :   }
     951      103469 :   *pt_lambda = lambda;
     952      103469 :   return mkvec2(P,Q);
     953             : }
     954             : 
     955             : static long
     956      147407 : test53(long lambda, long ep, long g)
     957             : {
     958      147407 :   return (lambda <= g+1) || (odd(g) && lambda<g+3 && ep==1);
     959             : }
     960             : 
     961             : static long
     962      189648 : test55(GEN W, long ep, long g)
     963             : {
     964      189648 :   GEN P = gel(W,1), Q = gel(W,2);
     965      189648 :   GEN Pe = FpX_red(ep ? ZX_shifti(P,-1): P, gen_2);
     966      189648 :   GEN Qe = FpX_red(ep ? ZX_shifti(Q,-1): Q, gen_2);
     967      189648 :   if (ep==0)
     968             :   {
     969      149255 :     if (signe(Qe)!=0) return ZX_val(Qe) >= (g + 3)>>1;
     970       90814 :     else return ZX_val(FpX_deriv(Pe, gen_2)) >= g+1;
     971             :   }
     972             :   else
     973       40393 :     return ZX_val(Qe) >= (g+1)>>1 && ZX_val(Pe) >= g + 1;
     974             : }
     975             : 
     976             : static GEN
     977       50840 : hyperell_reverse(GEN W, long g)
     978             : {
     979       50840 :   return mkvec2(RgXn_recip_shallow(gel(W,1),2*g+3),
     980       50840 :                 RgXn_recip_shallow(gel(W,2),g+2));
     981             : }
     982             : 
     983             : /* [P,Q] -> [P(2x)/4^r, Q(2x)/2^r] */
     984             : static GEN
     985      169774 : ZX2_unscale(GEN W, long r)
     986             : {
     987      169774 :   GEN P = ZX_unscale2n(gel(W,1), 1);
     988      169774 :   GEN Q = ZX_unscale2n(gel(W,2), 1);
     989      169774 :   if (r)
     990             :   {
     991       30931 :     P = ZX_shifti(P, -2*r);
     992       30931 :     Q = ZX_shifti(Q, -r);
     993             :   }
     994      169774 :   return mkvec2(P,Q);
     995             : }
     996             : /* [P,Q] -> [P(2x+c)/4^r, Q(2x+c)/2^r] */
     997             : static GEN
     998      163756 : ZX2_affine_unscale(GEN W, long c, long r)
     999             : {
    1000      238886 :   if (c) W = mkvec2(ZX_Z_translate(gel(W,1), gen_1),
    1001       75130 :                     ZX_Z_translate(gel(W,2), gen_1));
    1002      163756 :   return ZX2_unscale(W, r);
    1003             : }
    1004             : 
    1005             : static GEN
    1006       50805 : algo56(GEN W, long g)
    1007             : {
    1008             :   long ep;
    1009       50805 :   GEN M = mkvec2(gen_1, matid(2)), Woo;
    1010       50805 :   W = algo51(W, M);
    1011       50805 :   Woo = hyperell_reverse(W, g);
    1012       50805 :   ep = get_ep(Woo);
    1013       50805 :   if (test55(Woo,ep,g))
    1014             :   {
    1015             :     long lambda;
    1016       11744 :     Woo = algo52(Woo, gen_0, &lambda);
    1017       11744 :     if (!test53(lambda,ep,g))
    1018             :     {
    1019        5969 :       long r = lambda>>1;
    1020        5969 :       gel(M,1) = shifti(gel(M,1), r);
    1021        5969 :       gel(M,2) = ZM2_mul(gel(M,2), mkmat22(gen_0, gen_1, gen_2, gen_0));
    1022        5969 :       W = ZX2_unscale(Woo, r);
    1023             :     }
    1024             :   }
    1025             :   for(;;)
    1026       24892 :   {
    1027       75697 :     long j, ep = get_ep(W);
    1028      189403 :     for (j = 0; j < 2; j++)
    1029      138598 :       if (test55(ZX2_affine_unscale(W, j, 0), ep, g))
    1030             :       {
    1031             :         long lambda;
    1032       91578 :         GEN c = utoi(j), Wc = algo52(W, c, &lambda);
    1033       91578 :         if (!test53(lambda,ep,g))
    1034             :         {
    1035       24892 :           long r = lambda>>1;
    1036       24892 :           gel(M,1) = shifti(gel(M,1), r);
    1037       24892 :           gel(M,2) = ZM2_mul(gel(M,2), mkmat22(gen_2, c, gen_0, gen_1));
    1038       24892 :           W = ZX2_affine_unscale(Wc, j, r);
    1039       24892 :           break;
    1040             :         }
    1041             :       }
    1042       75697 :     if (j==2) break;
    1043             :   }
    1044       50805 :   return mkvec2(W, M);
    1045             : }
    1046             : 
    1047             : static GEN
    1048         105 : algo56bis(GEN W, long g, long inf, long thr)
    1049             : {
    1050         105 :   pari_sp av = avma;
    1051         105 :   GEN vl = cgetg(3,t_VEC);
    1052         105 :   long nl = 1;
    1053         105 :   W = algo51(W, NULL);
    1054         105 :   if (inf)
    1055             :   {
    1056          35 :     GEN Woo = hyperell_reverse(W, g);
    1057          35 :     long ep = get_ep(Woo);
    1058          35 :     if (test55(ZX2_unscale(Woo, 0), ep, g))
    1059             :     {
    1060             :       long lambda;
    1061          28 :       Woo = algo52(Woo, gen_0, &lambda);
    1062          28 :       if (lambda == thr) gel(vl,nl++) = ZX2_unscale(Woo, lambda>>1);
    1063             :     }
    1064             :   }
    1065             :   {
    1066         105 :     long j, ep = get_ep(W);
    1067         315 :     for (j = 0; j < 2; j++)
    1068         210 :       if (test55(ZX2_affine_unscale(W, j, 0), ep, g))
    1069             :       {
    1070             :         long lambda;
    1071         119 :         GEN Wc = algo52(W, utoi(j), &lambda);
    1072         119 :         if (lambda == thr) gel(vl,nl++) = ZX2_affine_unscale(Wc, j, lambda>>1);
    1073             :       }
    1074             :   }
    1075         105 :   setlg(vl, nl);
    1076         105 :   return gc_GEN(av,vl);
    1077             : }
    1078             : 
    1079             : /* return the (degree 2) apolar invariant (the nth transvectant of P and P) */
    1080             : static GEN
    1081          91 : ZX_apolar(GEN P, long n)
    1082             : {
    1083          91 :   pari_sp av = avma;
    1084          91 :   long d = degpol(P), i;
    1085          91 :   GEN s = gen_0, g = cgetg(n+2,t_VEC);
    1086          91 :   gel(g,1) = gen_1;
    1087         637 :   for (i = 1; i <= n; i++) gel(g,i+1) = muliu(gel(g,i),i); /* g[i+1] = i! */
    1088         714 :   for (i = n-d; i <= d; i++)
    1089             :   {
    1090         623 :      GEN a = mulii(mulii(gel(g,i+1),gel(g,n-i+1)),
    1091         623 :                    mulii(gel(P,i+2),gel(P,n-i+2)));
    1092         623 :      s = odd(i)? subii(s, a): addii(s, a);
    1093             :   }
    1094          91 :   return gc_INT(av,s);
    1095             : }
    1096             : 
    1097             : static GEN
    1098       53038 : algo57(GEN F, long g, GEN pr)
    1099             : {
    1100             :   long i, l;
    1101       53038 :   GEN D, C = content(F);
    1102       53038 :   GEN e = gel(core2(shifti(C,-vali(C))),2);
    1103       53038 :   GEN M = mkvec2(e, matid(2));
    1104       53038 :   long minvd = (2*g+1)>>(odd(g) ? 4:2);
    1105       53038 :   F = ZX_Z_divexact(F, sqri(e));
    1106       53038 :   D = absi(hyperelldisc(F));
    1107       53038 :   if (!pr)
    1108             :   {
    1109          91 :     GEN A = gcdii(D, ZX_apolar(F, 2*g+2));
    1110          91 :     pr = gel(factor(shifti(A, -vali(A))),1);
    1111             :   }
    1112       53038 :   l = lg(pr);
    1113      312821 :   for (i = 1; i < l; i++)
    1114             :   {
    1115             :     long ep;
    1116      259783 :     GEN p = gel(pr, i), ps2 = shifti(p,-1), Fe;
    1117      259783 :     if (equaliu(p,2) || Z_pval(D,p) < minvd) continue;
    1118      197547 :     ep = ZX_pvalrem(F,p, &Fe); Fe = FpX_red(Fe, p);
    1119      197547 :     if (degpol(Fe) < g+1+ep)
    1120             :     {
    1121        6406 :       GEN Fi = ZX_unscale(RgXn_recip_shallow(F,2*g+3), p);
    1122        6406 :       long lambda = ZX_pval(Fi,p);
    1123        6406 :       if (!test53(lambda,ep,g))
    1124             :       {
    1125        3815 :         GEN ppr = powiu(p,lambda>>1);
    1126        3815 :         F = ZX_Z_divexact(Fi,sqri(ppr));
    1127        3815 :         gel(M,1) = mulii(gel(M,1), ppr);
    1128        3815 :         gel(M,2) = ZM2_mul(gel(M,2), mkmat22(gen_0,gen_1,p,gen_0));
    1129             :       }
    1130             :     }
    1131             :     for(;;)
    1132       25186 :     {
    1133             :       GEN Fe, R;
    1134      222733 :       long j, lR, ep = ZX_pvalrem(F,p, &Fe);
    1135      222733 :       R = FpX_roots_mult(FpX_red(Fe, p), g+2-ep, p); lR = lg(R);
    1136      235226 :       for (j = 1; j<lR; j++)
    1137             :       {
    1138       37679 :         GEN c = Fp_center(gel(R,j), p, ps2);
    1139       37679 :         GEN Fi = ZX_affine(F,p,c);
    1140       37679 :         long lambda = ZX_pval(Fi,p);
    1141       37679 :         if (!test53(lambda,ep,g))
    1142             :         {
    1143       25186 :           GEN ppr = powiu(p,lambda>>1);
    1144       25186 :           F = ZX_Z_divexact(Fi, sqri(ppr));
    1145       25186 :           gel(M,1) = mulii(gel(M,1), ppr);
    1146       25186 :           gel(M,2) = ZM2_mul(gel(M,2), mkmat22(p,c,gen_0,gen_1));
    1147       25186 :           break;
    1148             :         }
    1149             :       }
    1150      222733 :       if (j==lR) break;
    1151             :     }
    1152             :   }
    1153       53038 :   return mkvec2(F, M);
    1154             : }
    1155             : 
    1156             : /* if inf=0, ignore point at infinity */
    1157             : static GEN
    1158        3080 : algo57bis(GEN F, long g, GEN p, long inf, long thr)
    1159             : {
    1160        3080 :   pari_sp av = avma;
    1161        3080 :   GEN vl = cgetg(3,t_VEC), Fe;
    1162        3080 :   long nl = 1, ep = ZX_pvalrem(F,p, &Fe);
    1163        3080 :   Fe = FpX_red(Fe, p);
    1164             :   {
    1165        3080 :     GEN R = FpX_roots_mult(Fe, thr-ep, p);
    1166        3080 :     long j, lR = lg(R);
    1167        5999 :     for (j = 1; j<lR; j++)
    1168             :     {
    1169        2919 :       GEN Fj = ZX_affine(F, p, gel(R,j));
    1170        2919 :       long lambda = ZX_pvalrem(Fj, p, &Fj);
    1171        2919 :       if (lambda == thr) gel(vl,nl++) = odd(lambda)? ZX_Z_mul(Fj, p): Fj;
    1172             :     }
    1173             :   }
    1174        3080 :   if (inf==1 && 2*g+2-degpol(Fe) >= thr-ep)
    1175             :   {
    1176           0 :     GEN Fj = ZX_unscale(RgXn_recip_shallow(F,2*g+3), p);
    1177           0 :     long lambda = ZX_pvalrem(Fj, p, &Fj);
    1178           0 :     if (lambda == thr) gel(vl,nl++) = odd(lambda)? ZX_Z_mul(Fj, p): Fj;
    1179             :   }
    1180        3080 :   setlg(vl, nl);
    1181        3080 :   return gc_GEN(av,vl);
    1182             : }
    1183             : 
    1184             : static GEN
    1185        3185 : next_model(GEN G, long g, GEN p, long inf, long thr)
    1186             : {
    1187        3290 :   return equaliu(p,2) ? algo56bis(G, g,    inf, thr)
    1188        3290 :                       : algo57bis(G, g, p, inf, thr);
    1189             : }
    1190             : 
    1191             : static GEN
    1192        1498 : get_extremal_even(GEN F, GEN G, long g, GEN p, long *nb)
    1193             : {
    1194             :   while (1)
    1195        1274 :   {
    1196        1498 :     GEN Wi = next_model(G, g, p, 0, g+2);
    1197        1498 :     if (lg(Wi)==1) return F;
    1198        1365 :     F = gel(Wi,1); ++*nb;
    1199        1365 :     if (DEBUGLEVEL>1) err_printf("model %ld: %Ps\n", *nb, F);
    1200        1365 :     Wi = next_model(F, g, p, 0, g+1);
    1201        1365 :     if (lg(Wi)==1) return F;
    1202        1274 :     G = gel(Wi,1);
    1203             :   }
    1204             : }
    1205             : 
    1206             : static GEN
    1207           0 : get_extremal_odd(GEN F, long g, GEN p, long *nb)
    1208             : {
    1209             :   while (1)
    1210           0 :   {
    1211           0 :     GEN Wi = next_model(F, g, p, 0, g+2);
    1212           0 :     if (lg(Wi)==1) return F;
    1213           0 :     F = gel(Wi,1); ++*nb;
    1214           0 :     if (DEBUGLEVEL>1) err_printf("model %ld: %Ps\n", *nb, F);
    1215             :   }
    1216             : }
    1217             : 
    1218             : static GEN
    1219         357 : hyperellextremalmodels_nb(GEN F, long g, GEN p, long *nb)
    1220             : {
    1221         357 :   pari_sp av = avma;
    1222             :   GEN W, A, B;
    1223             :   long l;
    1224             : 
    1225         357 :   *nb = 1;
    1226         357 :   if (equaliu(p,2))
    1227             :   {
    1228          35 :     if (get_ep(F) > 0) retmkvec(gcopy(F));
    1229             :   } else
    1230             :   {
    1231         322 :     F = check_hyperell(F);
    1232         322 :     if (ZX_pval(F, p) > 0) return gc_GEN(av, mkvec(F));
    1233             :   }
    1234         322 :   if (DEBUGLEVEL>1) err_printf("model %ld: %Ps\n", *nb, F);
    1235         322 :   W = next_model(F, g, p, 1, odd(g)? g+2: g+1);
    1236         322 :   l = lg(W); if (l==1) return gc_GEN(av, mkvec(F));
    1237         210 :   if (odd(g))
    1238             :   {
    1239           0 :     *nb = l-1;
    1240           0 :     A = get_extremal_odd(gel(W,1), g, p, nb);
    1241           0 :     B = l==3 ? get_extremal_odd(gel(W,2), g, p, nb) : F;
    1242             :   }
    1243             :   else
    1244             :   {
    1245         210 :     A = get_extremal_even(F, gel(W,1), g, p, nb);
    1246         210 :     B = l==3 ? get_extremal_even(F, gel(W,2), g, p, nb) : F;
    1247             :   }
    1248         210 :   return gc_GEN(av, A == B? mkvec(A): mkvec2(A, B));
    1249             : }
    1250             : 
    1251             : static GEN
    1252         350 : hyperellextremalmodels_i(GEN F, long g, GEN p)
    1253             : {
    1254             :   long nb;
    1255         350 :   return hyperellextremalmodels_nb(F, g, p, &nb);
    1256             : }
    1257             : 
    1258             : GEN
    1259           7 : hyperellextremalmodels(GEN PQ, GEN p)
    1260             : {
    1261           7 :   pari_sp av = avma;
    1262           7 :   GEN H = check_hyperell(PQ), W, v;
    1263             :   long g, nb;
    1264           7 :   if (!H || signe(H)==0) pari_err_TYPE("hyperellextremalmodels",PQ);
    1265           7 :   if (typ(p)!=t_INT || signe(p)<=0) pari_err_TYPE("hyperellextremalmodels",p);
    1266           7 :   g = hyperellgenus(H);
    1267           7 :   W = hyperellminimalmodel(H,NULL,mkvec(p));
    1268           7 :   v = cgetg(3, t_VEC);
    1269           7 :   gel(v, 2) = hyperellextremalmodels_nb(W, g, p, &nb);
    1270           7 :   gel(v, 1) = stoi(nb);
    1271           7 :   return gc_upto(av, v);
    1272             : }
    1273             : 
    1274             : static GEN
    1275      303761 : RgX_RgM2_eval(GEN P, GEN A, GEN Bp, long d)
    1276             : {
    1277      303761 :   if (signe(P)==0)
    1278       79841 :     return P;
    1279             :   else
    1280             :   {
    1281      223920 :     long dP = degpol(P);
    1282      223920 :     GEN R = RgX_homogenous_evalpow(P, A, Bp);
    1283      223920 :     if (d > dP)
    1284       17298 :       R = gmul(R, gel(Bp,1+d-dP));
    1285      223920 :     return R;
    1286             :   }
    1287             : }
    1288             : 
    1289             : static GEN
    1290       53052 : minimalmodel_merge(GEN W2, GEN Modd, long g, long v)
    1291             : {
    1292       53052 :   GEN P = gel(W2,1), Q = gel(W2,2);
    1293       53052 :   GEN e = gel(Modd,1), M = gel(Modd,2);
    1294       53052 :   GEN A = deg1pol_shallow(gcoeff(M,1,1), gcoeff(M,1,2), v);
    1295       53052 :   GEN B = deg1pol_shallow(gcoeff(M,2,1), gcoeff(M,2,2), v);
    1296       53052 :   GEN Bp = gpowers(B, 2*g+2);
    1297       53052 :   long f = mod4(e)==1 ? 1: -1;
    1298       53052 :   GEN m = shifti(f > 0 ? subui(1,e): addui(1,e), -2);
    1299       53052 :   GEN  m24 = subii(shifti(m,1), shifti(sqri(m),2));
    1300       53052 :   P = RgX_RgM2_eval(P, A, Bp, 2*g+2);
    1301       53052 :   Q = RgX_RgM2_eval(Q, A, Bp, g+1);
    1302       53052 :   P = ZX_Z_divexact(ZX_add(P, ZX_Z_mul(ZX_sqr(Q), m24)),sqri(e));
    1303       53052 :   if (f < 0) Q = ZX_neg(Q);
    1304       53052 :   return mkvec2(P,Q);
    1305             : }
    1306             : 
    1307             : static GEN
    1308      106090 : hyperell_redQ(GEN W)
    1309             : {
    1310      106090 :   GEN P = gel(W,1), Q = gel(W,2);
    1311      106090 :   GEN Pr, Qr = FpX_red(Q, gen_2);
    1312      106090 :   Pr = ZX_add(P, ZX_shifti(ZX_mul(ZX_sub(Q, Qr),ZX_add(Q, Qr)),-2));
    1313      106090 :   return mkvec2(Pr, Qr);
    1314             : }
    1315             : 
    1316             : static GEN
    1317       50735 : minimalmodel_getH(GEN W, GEN Qn, GEN e, GEN M, long g, long v)
    1318             : {
    1319       50735 :   GEN Q = gel(W,2);
    1320       50735 :   GEN A = deg1pol_shallow(gcoeff(M,1,1), gcoeff(M,1,2), v);
    1321       50735 :   GEN B = deg1pol_shallow(gcoeff(M,2,1), gcoeff(M,2,2), v);
    1322       50735 :   GEN Bp = gpowers(B, g+1);
    1323       50735 :   return ZX_shifti(ZX_sub(ZX_Z_mul(Qn,e),RgX_RgM2_eval(Q, A, Bp, g+1)), -1);
    1324             : }
    1325             : 
    1326             : static void
    1327       53094 : check_hyperell_Q(const char *fun, GEN *pW, GEN *pF)
    1328             : {
    1329       53094 :   GEN W = *pW, F = check_hyperell(W);
    1330             :   long v, g;
    1331       53094 :   if (!F || !signe(F) || !RgX_is_ZX(F)) pari_err_TYPE(fun, W);
    1332       53087 :   if (!signe(ZX_disc(F))) pari_err_DOMAIN(fun,"disc(W)","==",gen_0,W);
    1333       53080 :   v = varn(F); g = hyperellgenus(F);
    1334       53080 :   if (g == 0) pari_err_DOMAIN(fun, "genus", "=", gen_0, gen_0);
    1335       53066 :   if (typ(W)==t_POL) W = mkvec2(W, pol_0(v));
    1336             :   else
    1337             :   {
    1338       43610 :     GEN P = gel(W, 1), Q = gel(W, 2);
    1339       43610 :     if (typ(P)!=t_POL) P = scalarpol_shallow(P, v);
    1340       43610 :     if (typ(Q)!=t_POL) Q = scalarpol_shallow(Q, v);
    1341       43610 :     if (!RgX_is_ZX(P) || !RgX_is_ZX(Q)) pari_err_TYPE(fun,W);
    1342       43610 :     if (degpol(P) > 2*g+2) pari_err_DOMAIN(fun, "deg(P)", ">", utoi(2*g+2), P);
    1343       43610 :     if (degpol(Q) > g+1) pari_err_DOMAIN(fun, "deg(Q)", ">", utoi(g+1), Q);
    1344       43610 :     W = mkvec2(P, Q);
    1345             :   }
    1346       53066 :   *pW = W; *pF = F;
    1347       53066 : }
    1348             : 
    1349             : GEN
    1350       53052 : hyperellminimalmodel(GEN W, GEN *pM, GEN pr)
    1351             : {
    1352       53052 :   pari_sp av = avma;
    1353             :   GEN Wr, F, WM2, F2, W2, M2, Modd, Wf, ef, Mf, Hf;
    1354             :   long g, v;
    1355       53052 :   check_hyperell_Q("hyperellminimalmodel",&W, &F);
    1356       53052 :   if (pr && (!is_vec_t(typ(pr)) || !RgV_is_ZV(pr)))
    1357          14 :     pari_err_TYPE("hyperellminimalmodel",pr);
    1358       53038 :   g = hyperellgenus(F); v = varn(F);
    1359       53038 :   Wr = hyperell_redQ(W);
    1360       53038 :   if (!pr || RgV_isin(pr, gen_2))
    1361             :   {
    1362       50805 :     WM2 = algo56(Wr,g); W2 = gel(WM2, 1); M2 = gel(WM2, 2);
    1363       50805 :     F2 = check_hyperell(W2);
    1364             :   }
    1365             :   else
    1366             :   {
    1367        2233 :     W2 = Wr; F2 = F; M2 = mkvec2(gen_1, matid(2));
    1368             :   }
    1369       53038 :   Modd = gel(algo57(F2, g, pr), 2);
    1370       53038 :   Wf = hyperell_redQ(minimalmodel_merge(W2, Modd, g, v));
    1371       53038 :   if (!pM) return gc_GEN(av, Wf);
    1372       50721 :   ef = mulii(gel(M2,1), gel(Modd,1));
    1373       50721 :   Mf = ZM2_mul(gel(M2,2), gel(Modd,2));
    1374       50721 :   Hf = minimalmodel_getH(W, gel(Wf,2), ef, Mf, g, v);
    1375       50721 :   *pM =  mkvec3(ef, Mf, Hf);
    1376       50721 :   return gc_all(av, 2, &Wf, pM);
    1377             : }
    1378             : 
    1379             : GEN
    1380          14 : hyperellminimaldisc(GEN W, GEN pr)
    1381             : {
    1382          14 :   pari_sp av = avma;
    1383          14 :   GEN C = hyperellminimalmodel(W, NULL, pr);
    1384          14 :   return gc_INT(av, hyperelldisc(C));
    1385             : }
    1386             : 
    1387             : static GEN
    1388          35 : redqfbsplit(GEN a, GEN b, GEN c, GEN d)
    1389             : {
    1390          35 :   GEN p = subii(d,b), q = shifti(a,1);
    1391          35 :   GEN U, Q, u, v, w = bezout(p, q, &u, &v);
    1392             : 
    1393          35 :   if (!equali1(w)) { p = diviiexact(p, w); q = diviiexact(q, w); }
    1394          35 :   U = mkmat22(p, negi(v), q, u);
    1395          35 :   Q = qfb3_SL2_apply(mkvec3(a,b,c), U);
    1396          35 :   b = gel(Q, 2); c = gel(Q,3);
    1397          35 :   if (signe(b) < 0) gel(U,2) = mkcol2(v, negi(u));
    1398          35 :   gel(U,2) = ZC_lincomb(gen_1, truedivii(negi(c), d), gel(U,2), gel(U,1));
    1399          35 :   return U;
    1400             : }
    1401             : 
    1402             : static GEN
    1403       16386 : polreduce(GEN P, GEN M)
    1404             : {
    1405       16386 :   long v = varn(P), dP = degpol(P), d = odd(dP) ? dP+1: dP;
    1406       16386 :   GEN A = deg1pol_shallow(gcoeff(M,1,1), gcoeff(M,1,2), v);
    1407       16386 :   GEN B = deg1pol_shallow(gcoeff(M,2,1), gcoeff(M,2,2), v);
    1408       16386 :   return RgX_RgM2_eval(P, A, gpowers(B, d), d);
    1409             : }
    1410             : 
    1411             : /* assume deg(P) > 2 */
    1412             : static GEN
    1413        8193 : red_Cremona_Stoll(GEN P, GEN *pM)
    1414             : {
    1415             :   GEN q1, q2, q3, M, R;
    1416        8193 :   long i, prec = nbits2prec(2*gexpo(P)) + EXTRAPRECWORD, d = degpol(P);
    1417        8193 :   GEN dP = ZX_deriv(P);
    1418             :   for (;;)
    1419           0 :   {
    1420        8193 :     GEN r = QX_complex_roots(P, prec);
    1421        8193 :     q1 = gen_0; q2 = gen_0; q3 = gen_0;
    1422       41000 :     for (i = 1; i <= d; i++)
    1423             :     {
    1424       32807 :       GEN ri = gel(r,i);
    1425       32807 :       GEN s = ginv(gabs(RgX_cxeval(dP,ri,NULL), prec));
    1426       32807 :       if (d!=4) s = gpow(s, gdivgs(gen_2,d-2), prec);
    1427       32807 :       q1 = gadd(q1, s);
    1428       32807 :       q2 = gsub(q2, gmul(real_i(ri), s));
    1429       32807 :       q3 = gadd(q3, gmul(gnorm(ri), s));
    1430             :     }
    1431        8193 :     M = lllgram(mkmat22(q1,q2,q2,q3));
    1432        8193 :     if (M && lg(M) == 3) break;
    1433           0 :     prec = precdbl(prec);
    1434             :   }
    1435        8193 :   R = polreduce(P, M);
    1436        8193 :   *pM = M;
    1437        8193 :   return R;
    1438             : }
    1439             : 
    1440             : /* assume deg(P) > 2 */
    1441             : GEN
    1442        8193 : ZX_hyperellred(GEN P, GEN *pM)
    1443             : {
    1444        8193 :   pari_sp av = avma;
    1445        8193 :   long d = degpol(P);
    1446             :   GEN q1, q2, q3, D, vD;
    1447        8193 :   GEN a = gel(P,d+2), b = gel(P,d+1), c = gel(P, d);
    1448             :   GEN M, R, M2;
    1449             : 
    1450        8193 :   q1 = muliu(sqri(a), d);
    1451        8193 :   q2 = shifti(mulii(a,b), 1);
    1452        8193 :   q3 = subii(sqri(b), shifti(mulii(a,c), 1));
    1453        8193 :   D = gcdii(gcdii(q1, q2), q3);
    1454        8193 :   if (!equali1(D))
    1455             :   {
    1456        8172 :     q1 = diviiexact(q1, D);
    1457        8172 :     q2 = diviiexact(q2, D);
    1458        8172 :     q3 = diviiexact(q3, D);
    1459             :   }
    1460        8193 :   D = qfb_disc3(q1, q2, q3);
    1461        8193 :   if (!signe(D))
    1462          49 :     M = mkmat22(gen_1, truedivii(negi(q2),shifti(q1,1)), gen_0, gen_1);
    1463        8144 :   else if (issquareall(D,&vD))
    1464          35 :     M = redqfbsplit(q1, q2, q3, vD);
    1465             :   else
    1466        8109 :     M = gel(qfbredsl2(mkqfb(q1,q2,q3,D), NULL), 2);
    1467        8193 :   R = red_Cremona_Stoll(polreduce(P, M), &M2);
    1468        8193 :   if (pM) *pM = gmul(M, M2);
    1469        8193 :   return gc_all(av, pM ? 2: 1, &R, pM);
    1470             : }
    1471             : 
    1472             : GEN
    1473          42 : hyperellred(GEN W, GEN *pM)
    1474             : {
    1475          42 :   pari_sp av = avma;
    1476             :   long g, v;
    1477             :   GEN F, M, Wf, Hf;
    1478          42 :   check_hyperell_Q("hyperellred", &W, &F);
    1479          14 :   g = hyperellgenus(F); v = varn(F);
    1480          14 :   (void) ZX_hyperellred(F, &M);
    1481          14 :   Wf = hyperell_redQ(minimalmodel_merge(W, mkvec2(gen_1, M), g, v));
    1482          14 :   Hf = minimalmodel_getH(W, gel(Wf,2), gen_1, M, g, v);
    1483          14 :   if (pM) *pM = mkvec3(gen_1, M, Hf);
    1484          14 :   return gc_all(av, pM ? 2: 1, &Wf, pM);
    1485             : }
    1486             : 
    1487             : static void
    1488       65296 : check_hyperell_Rg(const char *fun, GEN *pW, GEN *pF)
    1489             : {
    1490       65296 :   GEN W = *pW, F = check_hyperell(W);
    1491             :   long v;
    1492       65296 :   if (!F)
    1493           7 :     pari_err_TYPE(fun, W);
    1494       65289 :   if (degpol(F) <= 0) pari_err_CONSTPOL(fun);
    1495       65282 :   v = varn(F);
    1496       65282 :   if (typ(W)==t_POL) W = mkvec2(W, pol_0(v));
    1497             :   else
    1498             :   {
    1499       65254 :     GEN P = gel(W, 1), Q = gel(W, 2);
    1500       65254 :     long g = hyperellgenus(F);
    1501       65254 :     if( typ(P)!=t_POL) P = scalarpol(P, v);
    1502       65254 :     if( typ(Q)!=t_POL) Q = scalarpol(Q, v);
    1503       65254 :     if (degpol(P) > 2*g+2)
    1504           0 :       pari_err_DOMAIN(fun, "poldegree(P)", ">", utoi(2*g+2), P);
    1505       65254 :     if (degpol(Q) > g+1)
    1506           0 :       pari_err_DOMAIN(fun, "poldegree(Q)", ">", utoi(g+1), Q);
    1507             : 
    1508       65254 :     W = mkvec2(P, Q);
    1509             :   }
    1510       65282 :   if (pF) *pF = F;
    1511       65282 :   *pW = W;
    1512       65282 : }
    1513             : 
    1514             : static void
    1515       65282 : check_hyperell_vc(const char *fun, GEN C, long v, GEN *e, GEN *M, GEN *H)
    1516             : {
    1517       65282 :   if (typ(C) != t_VEC || lg(C) != 4) pari_err_TYPE(fun,C);
    1518       65275 :   *e = gel(C,1); *M = gel(C,2); *H = gel(C,3);
    1519       65275 :   if (typ(*M) != t_MAT || lg(*M) != 3 || lgcols(*M) != 3) pari_err_TYPE(fun,C);
    1520       65268 :   if (typ(*H)!=t_POL || varncmp(varn(*H),v) > 0) *H = scalarpol_shallow(*H,v);
    1521       65268 : }
    1522             : 
    1523             : GEN
    1524       65296 : hyperellchangecurve(GEN W, GEN C)
    1525             : {
    1526       65296 :   pari_sp av = avma;
    1527             :   GEN F, P, Q, A, B, Bp, e, M, H;
    1528             :   long g, v;
    1529       65296 :   check_hyperell_Rg("hyperellchangecurve",&W,&F);
    1530       65282 :   P = gel(W,1); Q = gel(W,2);
    1531       65282 :   g = hyperellgenus(F); v = varn(F);
    1532       65282 :   check_hyperell_vc("hyperellchangecurve", C, v, &e, &M, &H);
    1533       65268 :   if (varncmp(gvar(M),v) <= 0)
    1534           0 :     pari_err_PRIORITY("hyperellchangecurve",M,"<=",v);
    1535       65268 :   A = deg1pol_shallow(gcoeff(M,1,1), gcoeff(M,1,2), v);
    1536       65268 :   B = deg1pol_shallow(gcoeff(M,2,1), gcoeff(M,2,2), v);
    1537       65268 :   Bp = gpowers(B, 2*g+2);
    1538       65268 :   P = RgX_RgM2_eval(P, A, Bp, 2*g+2);
    1539       65268 :   Q = RgX_RgM2_eval(Q, A, Bp, g+1);
    1540       65268 :   P = RgX_Rg_div(RgX_sub(P, RgX_mul(H,RgX_add(Q,H))), gsqr(e));
    1541       65268 :   Q = RgX_Rg_div(RgX_add(Q, RgX_mul2n(H,1)), e);
    1542       65268 :   return gc_GEN(av, mkvec2(P,Q));
    1543             : }
    1544             : 
    1545             : /****************************************************************************/
    1546             : /***                                                                      ***/
    1547             : /***                        genus2charpoly                                ***/
    1548             : /***                                                                      ***/
    1549             : /****************************************************************************/
    1550             : 
    1551             : /* Half stable reduction */
    1552             : 
    1553             : static long
    1554         588 : Zst_val(GEN P, GEN f, GEN p, long vt, GEN *pR)
    1555             : {
    1556         588 :   pari_sp av = avma;
    1557         588 :   long v = varn(P);
    1558             :   while(1)
    1559        1260 :   {
    1560        1848 :     long i, j, dm = LONG_MAX;
    1561        1848 :     GEN Pm = NULL;
    1562        1848 :     long dP = degpol(P);
    1563        7532 :     for (i = 0; i <= minss(dP, dm); i++)
    1564             :     {
    1565        5684 :       GEN Py = gel(P, i+2);
    1566        5684 :       if (signe(Py))
    1567             :       {
    1568        4186 :         if (typ(Py)==t_POL)
    1569             :         {
    1570        3864 :           long dPy = degpol(Py);
    1571       12502 :           for (j = 0; j <= minss(dPy, dm-i); j++)
    1572             :           {
    1573        8638 :             GEN c = gel(Py, j+2);
    1574        8638 :             if (signe(c))
    1575             :             {
    1576        3556 :                 if (i+j < dm)
    1577             :                 {
    1578        1848 :                   dm = i+j;
    1579        1848 :                   Pm = monomial(gen_1, dm, v);
    1580        1848 :                   gel(Pm,dm+2) = gen_0;
    1581             :                 }
    1582        3556 :                 gel(Pm,i+2) = c;
    1583             :             }
    1584             :           }
    1585             :         } else
    1586             :         {
    1587         322 :           if (i < dm)
    1588             :           {
    1589          77 :             dm = i;
    1590          77 :             Pm = monomial(Py, dm, v);
    1591             :           }
    1592             :           else
    1593         245 :             gel(Pm, i+2) = Py;
    1594             :         }
    1595             :       }
    1596             :     }
    1597        1848 :     Pm = RgX_renormalize(Pm);
    1598        1848 :     if (ZX_pval(Pm,p)==0)
    1599             :     {
    1600         588 :       *pR = gc_GEN(av, P);
    1601         588 :       return dm;
    1602             :     }
    1603        1260 :     Pm = RgX_homogenize_deg(Pm, dm, vt);
    1604        1260 :     P = gadd(gsub(P, Pm), gmul(f, ZXX_Z_divexact(Pm, p)));
    1605             :   }
    1606             : }
    1607             : 
    1608             : static long
    1609         588 : Zst_normval(GEN P, GEN f, GEN p, long vt, GEN *pR)
    1610             : {
    1611         588 :   long v = Zst_val(P, f, p, vt, pR);
    1612         588 :   long e = RgX_val(*pR)>>1;
    1613         588 :   if (e > 0)
    1614             :   {
    1615           0 :     v -= 2*e;
    1616           0 :     *pR = RgX_shift(*pR, -2*e);
    1617             :   }
    1618         588 :   return v;
    1619             : }
    1620             : 
    1621             : static GEN
    1622        1176 : RgXY_swapsafe(GEN P, long v1, long v2)
    1623             : {
    1624        1176 :   if (varn(P)==v2)
    1625             :   {
    1626          77 :     P = shallowcopy(P); setvarn(P,v1); return P;
    1627             :   } else
    1628        1099 :     return RgXY_swap(P, RgXY_degreex(P), v2);
    1629             : }
    1630             : 
    1631             : static GEN
    1632         588 : Zst_red1(GEN P, GEN f, GEN p, long vt)
    1633             : {
    1634         588 :   pari_sp av = avma;
    1635             :   GEN r, f1, f2, P1, P2;
    1636         588 :   long vs = varn(P);
    1637         588 :   long w = Zst_normval(P, f, p, vt, &r), ww = w-odd(w);
    1638         588 :   GEN st = monomial(pol_x(vt), 1, vs);
    1639         588 :   f1 = gsubst(f, vt, st);
    1640         588 :   P1 = gsubst(gdiv(r, monomial(gen_1,ww,vs)),vt,st);
    1641         588 :   f2 = gsubst(f, vs, st);
    1642         588 :   P2 = gsubst(gdiv(r, monomial(gen_1,ww,vt)),vs,st);
    1643         588 :   f2 = RgXY_swapsafe(f2, vs, vt);
    1644         588 :   P2 = RgXY_swapsafe(P2, vs, vt);
    1645         588 :   return gc_GEN(av, mkvec4(P1, f1, P2, f2));
    1646             : }
    1647             : 
    1648             : static GEN
    1649        1176 : Zst_reduce(GEN P, GEN p, long vt, long *pv)
    1650             : {
    1651             :   GEN C;
    1652        1176 :   long v = RgX_val(P);
    1653        1176 :   *pv = v + ZXX_pvalrem(RgX_shift(P, -v), p, &P);
    1654        1176 :   C = constant_coeff(P);
    1655        1176 :   C = typ(C) == t_POL ? C: scalarpol_shallow(C, vt);
    1656        1176 :   return FpX_red(C, p);
    1657             : }
    1658             : 
    1659             : static GEN
    1660         588 : Zst_red3(GEN C, GEN p, long vt)
    1661             : {
    1662             :   while(1)
    1663         511 :   {
    1664         588 :     GEN P1 = gel(C,1) ,f1 = gel(C,2), Poo = gel(C,3), foo= gel(C,4);
    1665             :     long e;
    1666         588 :     GEN Qoop = Zst_reduce(Poo, p, vt, &e), Qp, R;
    1667         588 :     if (RgX_val(Qoop) >= 3-e)
    1668             :     {
    1669           0 :       C = Zst_red1(Poo, foo, p, vt);
    1670         511 :       continue;
    1671             :     }
    1672         588 :     Qp = Zst_reduce(P1, p, vt, &e);
    1673         588 :     R = FpX_roots_mult(Qp, 3-e, p);
    1674         588 :     if (lg(R) > 1)
    1675         511 :     {
    1676         511 :       GEN xz = deg1pol_shallow(gen_1, gel(R,1), vt);
    1677         511 :       C = Zst_red1(gsubst(P1, vt, xz), gsubst(f1, vt, xz), p, vt);
    1678         511 :       continue;
    1679             :     }
    1680          77 :     return Qp;
    1681             :   }
    1682             : }
    1683             : 
    1684             : static GEN
    1685          77 : genus2_halfstablemodel_i(GEN P, GEN p, long vt)
    1686             : {
    1687             :   GEN Qp, R, Poo, Qoop;
    1688          77 :   long e = ZX_pvalrem(P, p, &Qp);
    1689          77 :   R = FpX_roots_mult(FpX_red(Qp,p), 4-e, p);
    1690          77 :   if (lg(R) > 1)
    1691             :   {
    1692          77 :     GEN C = Zst_red1(ZX_Z_translate(P, gel(R,1)), pol_x(vt), p, vt);
    1693          77 :     return Zst_red3(C, p, vt);
    1694             :   }
    1695           0 :   Poo = RgXn_recip_shallow(P, 7);
    1696           0 :   e = ZX_pvalrem(Poo, p, &Qoop);
    1697           0 :   Qoop = FpX_red(Qoop,p);
    1698           0 :   if (RgX_val(Qoop)>=4-e)
    1699             :   {
    1700           0 :     GEN C = Zst_red1(Poo, pol_x(vt), p, vt);
    1701           0 :     return Zst_red3(C, p, vt);
    1702             :   }
    1703           0 :   return gcopy(P);
    1704             : }
    1705             : 
    1706             : static GEN
    1707          77 : genus2_halfstablemodel(GEN P, GEN p)
    1708             : {
    1709          77 :   pari_sp av = avma;
    1710          77 :   long vt = fetch_var(), vs = varn(P);
    1711          77 :   GEN S = genus2_halfstablemodel_i(P, p, vt);
    1712          77 :   setvarn(S, vs); delete_var();
    1713          77 :   return gc_GEN(av, S);
    1714             : }
    1715             : 
    1716             : /* semi-stable reduction */
    1717             : 
    1718             : static GEN
    1719         532 : genus2_redmodel(GEN P, GEN p)
    1720             : {
    1721             :   GEN LP, U, F;
    1722             :   long i, k, r;
    1723         532 :   if (degpol(P) < 0) return mkvec2(cgetg(1, t_COL), P);
    1724         497 :   F = FpX_factor_squarefree(P, p);
    1725         497 :   r = lg(F); U = NULL;
    1726        1967 :   for (i = k = 1; i < r; i++)
    1727             :   {
    1728        1470 :     GEN f = gel(F,i);
    1729        1470 :     long df = degpol(f);
    1730        1470 :     if (!df) continue;
    1731         777 :     if (odd(i)) U = U? FpX_mul(U, f, p): f;
    1732         777 :     if (i > 1) gel(F,k++) = df == 1? mkcol(f): gel(FpX_factor(f, p), 1);
    1733             :   }
    1734         497 :   LP = leading_coeff(P);
    1735         497 :   if (!U)
    1736         112 :     U = scalarpol_shallow(LP, varn(P));
    1737             :   else
    1738             :   {
    1739         385 :     GEN LU = leading_coeff(U);
    1740         385 :     if (!equalii(LU, LP)) U = FpX_Fp_mul(U, Fp_div(LP, LU, p), p);
    1741             :   }
    1742         497 :   setlg(F,k); if (k > 1) F = shallowconcat1(F);
    1743         497 :   return mkvec2(F, U);
    1744             : }
    1745             : 
    1746             : static GEN
    1747        2926 : xdminusone(long d)
    1748             : {
    1749        2926 :   return gsub(pol_xn(d, 0),gen_1);
    1750             : }
    1751             : 
    1752             : static GEN
    1753         140 : ellfromeqncharpoly(GEN P, GEN Q, GEN p)
    1754             : {
    1755             :   long v;
    1756             :   GEN E, F, t, y;
    1757         140 :   v = fetch_var();
    1758         140 :   y = pol_x(v);
    1759         140 :   F = gsub(gadd(ZX_sqr(y), gmul(y, Q)), P);
    1760         140 :   E = ellinit(ellfromeqn(F), p, DEFAULTPREC);
    1761         140 :   delete_var();
    1762         140 :   t = ellcharpoly(E, p);
    1763         140 :   obj_free(E);
    1764         140 :   return t;
    1765             : }
    1766             : 
    1767             : static GEN
    1768           0 : nfellcharpoly(GEN e, GEN T, GEN p)
    1769             : {
    1770             :   GEN nf, E, t;
    1771           0 :   e = shallowcopy(e);
    1772           0 :   nf = nfinit(mkvec2(T, mkvec(p)), DEFAULTPREC);
    1773             :   while(1)
    1774             :   {
    1775           0 :     E = ellinit(e, nf, DEFAULTPREC);
    1776           0 :     if (lg(E)!=1) break;
    1777           0 :     gel(e,5) = gadd(gel(e,5), p);
    1778             :   }
    1779           0 :   t = elleulerf(E, p);
    1780           0 :   obj_free(E);
    1781           0 :   return RgX_recip(ginv(t));
    1782             : }
    1783             : 
    1784             : static GEN
    1785           0 : genus2_red5(GEN P, GEN T, GEN p)
    1786             : {
    1787           0 :   long vx = varn(P), vy = varn(T);
    1788           0 :   GEN f = shallowcopy(T), pi = shifti(p,-1);
    1789           0 :   setvarn(f, vx);
    1790             :   while(1)
    1791           0 :   {
    1792             :     GEN Pr, R, r, Rs;
    1793           0 :     long v = ZXX_pvalrem(P, p, &Pr);
    1794           0 :     R = FpXQX_roots_mult(Pr, 2-v, T, p);
    1795           0 :     if (lg(R)==1) return P;
    1796           0 :     r = FpX_center(gel(R,1), p, pi);
    1797           0 :     Pr = RgX_affine(P, p, r);
    1798           0 :     setvarn(r, vx);
    1799           0 :     f = RgX_Rg_div(gsub(f, r), p);
    1800           0 :     Rs = RgX_rem(RgXY_swap(Pr, 3, vy), gsub(f, pol_x(vy)));
    1801           0 :     Pr = RgXY_swap(Rs, 3, vy);
    1802           0 :     if (ZXX_pvalrem(Pr, sqri(p), &Pr)==0) return P;
    1803           0 :     P = Pr;
    1804             :   }
    1805             : }
    1806             : 
    1807             : static GEN
    1808         350 : genus2_type5(GEN P, GEN p)
    1809             : {
    1810             :   GEN E, F, T, a, a2, Q;
    1811             :   long v;
    1812         350 :   if (equaliu(p, 2))
    1813          28 :     (void) ZXX_pvalrem(P, sqri(p), &P);
    1814         350 :   (void) ZX_pvalrem(P, p, &F);
    1815         350 :   F = FpX_red(F, p);
    1816         350 :   if (degpol(F) < 1) return NULL;
    1817         343 :   F = FpX_factor(F, p);
    1818         343 :   if (mael(F,2,1) != 3 || degpol(gmael(F,1,1)) != 2) return NULL;
    1819           0 :   T = gmael(F, 1, 1);
    1820           0 :   v = fetch_var_higher();
    1821           0 :   Q = RgV_to_RgX(ZX_digits(P, T), v);
    1822           0 :   Q = genus2_red5(Q, T, p);
    1823           0 :   a = gel(Q,5); a2 = ZX_sqr(a);
    1824           0 :   E = mkvec5(gen_0, gel(Q,4), gen_0, ZX_mul(gel(Q,3),a), ZX_mul(gel(Q,2),a2));
    1825           0 :   delete_var();
    1826           0 :   return nfellcharpoly(E, T, p);
    1827             : }
    1828             : 
    1829             : /* Assume P has semistable reduction at p */
    1830             : static GEN
    1831         532 : genus2_eulerfact_semistable(GEN P, GEN p)
    1832             : {
    1833         532 :   GEN Pp = FpX_red(P, p);
    1834         532 :   GEN GU = genus2_redmodel(Pp, p);
    1835         532 :   long d = 6-degpol(Pp), v = d/2, w = odd(d);
    1836             :   GEN abe, tor;
    1837         532 :   GEN ki, kp = pol_1(0), kq = pol_1(0);
    1838         532 :   GEN F = gel(GU,1), Q = gel(GU,2);
    1839         532 :   long dQ = degpol(Q), lF = lg(F)-1;
    1840             : 
    1841           7 :   abe = dQ >= 5 ? hyperellcharpoly(gmul(Q,gmodulo(gen_1,p)))
    1842        1057 :       : dQ >= 3 ? ellfromeqncharpoly(Q,gen_0,p)
    1843         525 :                 : pol_1(0);
    1844         420 :   ki = dQ != 0 ? xdminusone(1)
    1845         644 :               : Fp_issquare(gel(Q,2),p) ? ZX_sqr(xdminusone(1))
    1846         112 :                                         : xdminusone(2);
    1847         532 :   if (lF)
    1848             :   {
    1849             :     long i;
    1850        1064 :     for(i=1; i <= lF; i++)
    1851             :     {
    1852         616 :       GEN Fi = gel(F, i);
    1853         616 :       long d = degpol(Fi);
    1854         616 :       GEN e = FpX_rem(Q, Fi, p);
    1855         980 :       GEN kqf = lgpol(e)==0 ? xdminusone(d):
    1856         714 :                 FpXQ_issquare(e, Fi, p) ? ZX_sqr(xdminusone(d))
    1857         364 :                                         : xdminusone(2*d);
    1858         616 :       kp = gmul(kp, xdminusone(d));
    1859         616 :       kq = gmul(kq, kqf);
    1860             :     }
    1861             :   }
    1862         532 :   if (v)
    1863             :   {
    1864         245 :     GEN kqoo = w==1 ? xdminusone(1):
    1865           0 :                Fp_issquare(leading_coeff(Q), p)? ZX_sqr(xdminusone(1))
    1866           0 :                                               : xdminusone(2);
    1867         245 :     kp = gmul(kp, xdminusone(1));
    1868         245 :     kq = gmul(kq, kqoo);
    1869             :   }
    1870         532 :   tor = RgX_div(ZX_mul(xdminusone(1), kq), ZX_mul(ki, kp));
    1871         532 :   return ZX_mul(abe, tor);
    1872             : }
    1873             : 
    1874             : GEN
    1875         945 : genus2_eulerfact(GEN P, GEN p, long ra, long rt)
    1876             : {
    1877         945 :   pari_sp av = avma;
    1878             :   GEN W, R, E;
    1879         945 :   long d = 2*ra+rt;
    1880         945 :   if (d == 0) return pol_1(0);
    1881         322 :   R = genus2_type5(P, p);
    1882         322 :   if (R) return R;
    1883         322 :   W = hyperellextremalmodels_i(P, 2, p);
    1884         322 :   if (lg(W) < 3)
    1885             :   {
    1886         189 :     GEN F = genus2_eulerfact_semistable(P,p);
    1887         189 :     if (degpol(F)!=d)
    1888             :     {
    1889          77 :       GEN S = genus2_halfstablemodel(P, p);
    1890          77 :       F = genus2_eulerfact_semistable(S, p);
    1891          77 :       if (degpol(F)!=d) pari_err_BUG("genus2charpoly");
    1892             :     }
    1893         189 :     return F;
    1894             :   }
    1895         133 :   E =  gmul(genus2_eulerfact_semistable(gel(W,1),p),
    1896         133 :             genus2_eulerfact_semistable(gel(W,2),p));
    1897         133 :   return gc_upto(av, E);
    1898             : }
    1899             : 
    1900             : /*   p = 2  */
    1901             : 
    1902             : static GEN
    1903          28 : F2x_genus2_find_trans(GEN P, GEN Q, GEN F)
    1904             : {
    1905          28 :   pari_sp av = avma;
    1906          28 :   long i, d = F2x_degree(F), v = P[1];
    1907             :   GEN M, C, V;
    1908          28 :   M = cgetg(d+1, t_MAT);
    1909          84 :   for (i=1; i<=d; i++)
    1910             :   {
    1911          56 :     GEN Mi = F2x_rem(F2x_add(F2x_shift(Q,i-1), monomial_F2x(2*i-2,v)), F);
    1912          56 :     gel(M,i) = F2x_to_F2v(Mi, d);
    1913             :   }
    1914          28 :   C = F2x_to_F2v(F2x_rem(P, F), d);
    1915          28 :   V = F2m_F2c_invimage(M, C);
    1916          28 :   return gc_leaf(av, F2v_to_F2x(V, v));
    1917             : }
    1918             : 
    1919             : static GEN
    1920          42 : F2x_genus2_trans(GEN P, GEN Q, GEN H)
    1921             : {
    1922          42 :   return F2x_add(P,F2x_add(F2x_mul(H,Q), F2x_sqr(H)));
    1923             : }
    1924             : 
    1925             : static GEN
    1926         105 : F2x_genus_redoo(GEN P, GEN Q, long k)
    1927             : {
    1928         105 :   if (F2x_degree(P)==2*k)
    1929             :   {
    1930          21 :     long c = F2x_coeff(P,2*k-1), dQ = F2x_degree(Q);
    1931          21 :     if ((dQ==k-1 && c==1) || (dQ<k-1 && c==0))
    1932          14 :      return F2x_genus2_trans(P, Q, monomial_F2x(k, P[1]));
    1933             :   }
    1934          91 :   return P;
    1935             : }
    1936             : 
    1937             : static GEN
    1938          56 : F2x_pseudodisc(GEN P, GEN Q)
    1939             : {
    1940          56 :   GEN dP = F2x_deriv(P), dQ = F2x_deriv(Q);
    1941          56 :   return F2x_gcd(Q, F2x_add(F2x_mul(P, F2x_sqr(dQ)), F2x_sqr(dP)));
    1942             : }
    1943             : 
    1944             : static GEN
    1945          35 : F2x_genus_red(GEN P, GEN Q)
    1946             : {
    1947             :   long dP, dQ;
    1948             :   GEN F, FF;
    1949          35 :   P = F2x_genus_redoo(P, Q, 3);
    1950          35 :   P = F2x_genus_redoo(P, Q, 2);
    1951          35 :   P = F2x_genus_redoo(P, Q, 1);
    1952          35 :   dP = F2x_degree(P);
    1953          35 :   dQ = F2x_degree(Q);
    1954          35 :   FF = F = F2x_pseudodisc(P,Q);
    1955          56 :   while(F2x_degree(F)>0)
    1956             :   {
    1957          21 :     GEN M = gel(F2x_factor(F),1);
    1958          21 :     long i, l = lg(M);
    1959          49 :     for(i=1; i<l; i++)
    1960             :     {
    1961          28 :       GEN R = F2x_sqr(gel(M,i));
    1962          28 :       GEN H = F2x_genus2_find_trans(P, Q, R);
    1963          28 :       P = F2x_div(F2x_genus2_trans(P, Q, H), R);
    1964          28 :       Q = F2x_div(Q, gel(M,i));
    1965             :     }
    1966          21 :     F = F2x_pseudodisc(P, Q);
    1967             :   }
    1968          35 :   return mkvec4(P,Q,FF,mkvecsmall2(dP,dQ));
    1969             : }
    1970             : 
    1971             : /* Number of solutions of x^2+b*x+c */
    1972             : static long
    1973          21 : F2xqX_quad_nbroots(GEN b, GEN c, GEN T)
    1974             : {
    1975          21 :   if (lgpol(b) > 0)
    1976             :   {
    1977          14 :     GEN d = F2xq_div(c, F2xq_sqr(b, T), T);
    1978          14 :     return F2xq_trace(d, T)? 0: 2;
    1979             :   }
    1980             :   else
    1981           7 :     return 1;
    1982             : }
    1983             : 
    1984             : static GEN
    1985          35 : genus2_eulerfact2_semistable(GEN PQ)
    1986             : {
    1987          35 :   GEN V = F2x_genus_red(ZX_to_F2x(gel(PQ, 1)), ZX_to_F2x(gel(PQ, 2)));
    1988          35 :   GEN P = gel(V, 1), Q = gel(V, 2);
    1989          35 :   GEN F = gel(V, 3), v = gel(V, 4);
    1990             :   GEN abe, tor;
    1991          35 :   GEN ki, kp = pol_1(0), kq = pol_1(0);
    1992          35 :   long dP = F2x_degree(P), dQ = F2x_degree(Q), d = maxss(dP, 2*dQ);
    1993          35 :   if (!lgpol(F)) return pol_1(0);
    1994          35 :   ki = dQ!=0 || dP>0 ? xdminusone(1):
    1995           7 :       dP==-1 ? ZX_sqr(xdminusone(1)): xdminusone(2);
    1996          56 :   abe = d>=5? hyperellcharpoly(gmul(PQ,gmodulss(1,2))):
    1997          28 :         d>=3? ellfromeqncharpoly(F2x_to_ZX(P), F2x_to_ZX(Q), gen_2):
    1998          14 :         pol_1(0);
    1999          28 :   if (lgpol(F))
    2000             :   {
    2001          28 :     GEN M = gel(F2x_factor(F), 1);
    2002          28 :     long i, lF = lg(M)-1;
    2003          49 :     for(i=1; i <= lF; i++)
    2004             :     {
    2005          21 :       GEN Fi = gel(M, i);
    2006          21 :       long d = F2x_degree(Fi);
    2007          21 :       long nb  = F2xqX_quad_nbroots(F2x_rem(Q, Fi), F2x_rem(P, Fi), Fi);
    2008          35 :       GEN kqf = nb==1 ? xdminusone(d):
    2009           0 :                 nb==2 ? ZX_sqr(xdminusone(d))
    2010          14 :                       : xdminusone(2*d);
    2011          21 :       kp = gmul(kp, xdminusone(d));
    2012          21 :       kq = gmul(kq, kqf);
    2013             :     }
    2014             :   }
    2015          28 :   if (maxss(v[1],2*v[2])<5)
    2016             :   {
    2017          28 :     GEN kqoo = v[1]>2*v[2] ? xdminusone(1):
    2018           0 :                v[1]<2*v[2] ? ZX_sqr(xdminusone(1))
    2019           7 :                            : xdminusone(2);
    2020          21 :     kp = gmul(kp, xdminusone(1));
    2021          21 :     kq = gmul(kq, kqoo);
    2022             :   }
    2023          28 :   tor = RgX_div(ZX_mul(xdminusone(1),kq), ZX_mul(ki, kp));
    2024          28 :   return ZX_mul(abe, tor);
    2025             : }
    2026             : 
    2027             : GEN
    2028          28 : genus2_eulerfact2(GEN F, GEN PQ)
    2029             : {
    2030          28 :   pari_sp av = avma;
    2031          28 :   GEN W, R = genus2_type5(F, gen_2), E;
    2032          28 :   if (R) return R;
    2033          28 :   W = hyperellextremalmodels_i(PQ, 2, gen_2);
    2034          28 :   if (lg(W) < 3) return genus2_eulerfact2_semistable(PQ);
    2035           7 :   E = gmul(genus2_eulerfact2_semistable(gel(W,1)),
    2036           7 :            genus2_eulerfact2_semistable(gel(W,2)));
    2037           7 :   return gc_upto(av, E);
    2038             : }
    2039             : 
    2040             : GEN
    2041         889 : genus2charpoly(GEN G, GEN p)
    2042             : {
    2043         889 :   pari_sp av = avma;
    2044         889 :   GEN gr = genus2red(G, p), F;
    2045         889 :   GEN PQ = gel(gr, 3), L = gel(gr, 4), r = gel(L, 4);
    2046         889 :   GEN P = gadd(gsqr(gel(PQ, 2)), gmul2n(gel(PQ, 1), 2));
    2047         889 :   if (equaliu(p,2))
    2048           7 :     F = genus2_eulerfact2(P, PQ);
    2049             :   else
    2050         882 :     F = genus2_eulerfact(P,p, r[1],r[2]);
    2051         889 :   return gc_upto(av, F);
    2052             : }

Generated by: LCOV version 1.16