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 - base2.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.12.1 lcov report (development 24215-a79de5b25) Lines: 2167 2287 94.8 %
Date: 2019-08-24 05:50:50 Functions: 171 175 97.7 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2000  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. It is distributed in the hope that it will be useful, but WITHOUT
       8             : ANY WARRANTY WHATSOEVER.
       9             : Check the License for details. You should have received a copy of it, along
      10             : with the package; see the file 'COPYING'. If not, write to the Free Software
      11             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      12             : 
      13             : /*******************************************************************/
      14             : /*                                                                 */
      15             : /*                       MAXIMAL ORDERS                            */
      16             : /*                                                                 */
      17             : /*******************************************************************/
      18             : #include "pari.h"
      19             : #include "paripriv.h"
      20             : 
      21             : /* allow p = -1 from factorizations, avoid oo loop on p = 1 */
      22             : static long
      23        2198 : safe_Z_pvalrem(GEN x, GEN p, GEN *z)
      24             : {
      25        2198 :   if (is_pm1(p))
      26             :   {
      27          28 :     if (signe(p) > 0) return gvaluation(x,p); /*error*/
      28          21 :     *z = absi(x); return 1;
      29             :   }
      30        2170 :   return Z_pvalrem(x, p, z);
      31             : }
      32             : /* D an integer, P a ZV, return a factorization matrix for D over P, removing
      33             :  * entries with 0 exponent. */
      34             : static GEN
      35         609 : fact_from_factors(GEN D, GEN P, long flag)
      36             : {
      37         609 :   long i, l = lg(P), iq = 1;
      38         609 :   GEN Q = cgetg(l+1,t_COL);
      39         609 :   GEN E = cgetg(l+1,t_COL);
      40        2800 :   for (i=1; i<l; i++)
      41             :   {
      42        2198 :     GEN p = gel(P,i);
      43             :     long k;
      44        2198 :     if (flag && !equalim1(p))
      45             :     {
      46          14 :       p = gcdii(p, D);
      47          14 :       if (is_pm1(p)) continue;
      48             :     }
      49        2198 :     k = safe_Z_pvalrem(D, p, &D);
      50        2191 :     if (k) { gel(Q,iq) = p; gel(E,iq) = utoipos(k); iq++; }
      51             :   }
      52         602 :   D = absi_shallow(D);
      53         602 :   if (!equali1(D))
      54             :   {
      55         140 :     long k = Z_isanypower(D, &D);
      56         140 :     if (!k) k = 1;
      57         140 :     gel(Q,iq) = D; gel(E,iq) = utoipos(k); iq++;
      58             :   }
      59         602 :   setlg(Q,iq);
      60         602 :   setlg(E,iq); return mkmat2(Q,E);
      61             : }
      62             : 
      63             : /* d a t_INT; f a t_MAT factorisation of some t_INT sharing some divisors
      64             :  * with d, or a prime (t_INT). Return a factorization F of d: "primes"
      65             :  * entries in f _may_ be composite, and are included as is in d. */
      66             : static GEN
      67         616 : update_fact(GEN d, GEN f)
      68             : {
      69             :   GEN P;
      70         616 :   switch (typ(f))
      71             :   {
      72         602 :     case t_INT: case t_VEC: case t_COL: return f;
      73             :     case t_MAT:
      74          14 :       if (lg(f) == 3) { P = gel(f,1); break; }
      75             :     /*fall through*/
      76             :     default:
      77           7 :       pari_err_TYPE("nfbasis [factorization expected]",f);
      78             :       return NULL;/*LCOV_EXCL_LINE*/
      79             :   }
      80           7 :   return fact_from_factors(d, P, 1);
      81             : }
      82             : 
      83             : /* T = C T0(X/L); C = L^d / lt(T0), d = deg(T)
      84             :  * disc T = C^2(d - 1) L^-(d(d-1)) disc T0 = (L^d / lt(T0)^2)^(d-1) disc T0 */
      85             : static GEN
      86       21826 : set_disc(nfmaxord_t *S)
      87             : {
      88             :   GEN l0, L, dT;
      89             :   long d;
      90       21826 :   if (S->T0 == S->T) return ZX_disc(S->T);
      91        6909 :   d = degpol(S->T0);
      92        6909 :   l0 = leading_coeff(S->T0);
      93        6909 :   L = S->unscale;
      94        6909 :   if (typ(L) == t_FRAC && abscmpii(gel(L,1), gel(L,2)) < 0)
      95         826 :     dT = ZX_disc(S->T); /* more efficient */
      96             :   else
      97             :   {
      98        6083 :     GEN a = gpowgs(gdiv(gpowgs(L, d), sqri(l0)), d-1);
      99        6083 :     dT = gmul(a, ZX_disc(S->T0)); /* more efficient */
     100             :   }
     101        6909 :   return S->dT = dT;
     102             : }
     103             : 
     104             : /* dT != 0 */
     105             : static GEN
     106       16324 : poldiscfactors_i(GEN T, GEN dT, long flag)
     107             : {
     108       16324 :   GEN fa = absZ_factor_limit(dT, minuu(tridiv_bound(dT), maxprime()));
     109       16324 :   GEN Z, E, P = gel(fa,1), Tp = NULL;
     110       16324 :   long i, first = 1, l = lg(P);
     111       16324 :   if (l == 1) return fa;
     112       14966 :   setlg(P, l-1); Z = mkcol(gel(P, l-1));
     113       45024 :   while (lg(Z) != 1)
     114             :   { /* handle last element of Z */
     115       15092 :     GEN p = gel(Z, lg(Z)-1), r;
     116       15092 :     setlg(Z, lg(Z)-1);
     117       15092 :     if (!first) (void)Z_isanypower(p, &p);
     118       15092 :     first = 0;
     119       15092 :     if ((flag || lgefint(p)==3) && BPSW_psp(p))
     120       14936 :     { P = vec_append(P, p); continue; }
     121         156 :     if (!Tp) Tp = ZX_deriv(T);
     122         156 :     r = FpX_gcd_check(T, Tp, p);
     123         156 :     if (r)
     124          63 :       Z = shallowconcat(Z, Z_cba(r, diviiexact(p,r)));
     125          93 :     else if (flag)
     126          49 :       P = shallowconcat(P, gel(Z_factor(p),1));
     127             :     else
     128          44 :       P = vec_append(P, p);
     129             :   }
     130       14966 :   ZV_sort_inplace(P); l = lg(P); E = cgetg(l, t_COL);
     131       43869 :   for (i = 1; i < l; i++)
     132       28903 :     gel(E,i) = utoipos(Z_pvalrem(dT, gel(P,i), &dT));
     133       14966 :   return mkmat2(P,E);
     134             : }
     135             : 
     136             : GEN
     137          42 : poldiscfactors(GEN T, long flag)
     138             : {
     139          42 :   pari_sp av = avma;
     140             :   GEN dT;
     141          42 :   if (typ(T) != t_POL || !RgX_is_ZX(T)) pari_err_TYPE("poldiscfactors",T);
     142          42 :   if (flag < 0 || flag > 1) pari_err_FLAG("poldiscfactors");
     143          42 :   dT = ZX_disc(T);
     144          42 :   if (!signe(dT)) retmkvec2(gen_0, Z_factor(gen_0));
     145          35 :   return gerepilecopy(av, mkvec2(dT, poldiscfactors_i(T, dT, flag)));
     146             : }
     147             : 
     148             : static void
     149       21826 : nfmaxord_check_args(nfmaxord_t *S, GEN T, long flag)
     150             : {
     151       21826 :   GEN dT, L, E, P, fa = NULL;
     152             :   pari_timer t;
     153       21826 :   long l, ty = typ(T);
     154             : 
     155       21826 :   if (DEBUGLEVEL) timer_start(&t);
     156       21826 :   if (ty == t_VEC) {
     157        5537 :     if (lg(T) != 3) pari_err_TYPE("nfmaxord",T);
     158        5537 :     fa = gel(T,2); T = gel(T,1); ty = typ(T);
     159             :   }
     160       21826 :   if (ty != t_POL) pari_err_TYPE("nfmaxord",T);
     161       21826 :   T = Q_primpart(T);
     162       21826 :   if (degpol(T) <= 0) pari_err_CONSTPOL("nfmaxord");
     163       21826 :   RgX_check_ZX(T, "nfmaxord");
     164       21826 :   S->T0 = T;
     165       21826 :   T = ZX_Q_normalize(T, &L);
     166       21826 :   S->unscale = L;
     167       21826 :   S->T = T;
     168       21826 :   S->dT = dT = set_disc(S);
     169       21826 :   if (fa)
     170             :   {
     171        5537 :     const long MIN = 100; /* include at least all p < 101 */
     172             :     long tf;
     173        5537 :     if (!isint1(L)) fa = update_fact(dT, fa);
     174        5530 :     tf = typ(fa);
     175        5530 :     switch(tf)
     176             :     {
     177             :       case t_MAT:
     178         238 :         if (!is_Z_factornon0(fa)) pari_err_TYPE("nfmaxord",fa);
     179         231 :         fa = gel(fa,1); tf = t_COL; /* fall through */
     180             :       case t_VEC: case t_COL:
     181         602 :         P = gel(absZ_factor_limit(dT, MIN), 1); l = lg(P);
     182         602 :         if (l > 1 && abscmpiu(gel(P,1), MIN) <= 0)
     183             :         {
     184         504 :           if (abscmpiu(gel(P,l-1), MIN) > 0) setlg(P,l-1);
     185         504 :           settyp(P,tf); fa = ZV_sort_uniq(shallowconcat(fa,P));
     186             :         }
     187         602 :         fa = fact_from_factors(dT, fa, 0);
     188         595 :         break;
     189             :       case t_INT:
     190        4914 :         fa = absZ_factor_limit(dT, (signe(fa) <= 0)? 1: maxuu(itou(fa), MIN));
     191        4914 :         break;
     192             :       default:
     193           7 :         pari_err_TYPE("nfmaxord",fa);
     194             :     }
     195        5509 :     if (!signe(dT)) pari_err_IRREDPOL("nfmaxord",mkvec2(T,fa));
     196             :   }
     197             :   else
     198       16289 :     fa = poldiscfactors_i(T, dT, !(flag & nf_PARTIALFACT));
     199       21798 :   P = gel(fa,1); l = lg(P);
     200       21798 :   E = gel(fa,2);
     201       21798 :   if (l > 1 && is_pm1(gel(P,1)))
     202             :   {
     203          21 :     l--;
     204          21 :     P = vecslice(P, 2, l);
     205          21 :     E = vecslice(E, 2, l);
     206             :   }
     207       21798 :   S->dTP = P;
     208       21798 :   S->dTE = vec_to_vecsmall(E);
     209       21798 :   if (DEBUGLEVEL>2) timer_printf(&t, "disc. factorisation");
     210       21798 : }
     211             : 
     212             : static int
     213       57802 : fnz(GEN x,long j)
     214             : {
     215             :   long i;
     216      268745 :   for (i=1; i<j; i++)
     217      218098 :     if (signe(gel(x,i))) return 0;
     218       50647 :   return 1;
     219             : }
     220             : /* return list u[i], 2 by 2 coprime with the same prime divisors as ab */
     221             : static GEN
     222         105 : get_coprimes(GEN a, GEN b)
     223             : {
     224         105 :   long i, k = 1;
     225         105 :   GEN u = cgetg(3, t_COL);
     226         105 :   gel(u,1) = a;
     227         105 :   gel(u,2) = b;
     228             :   /* u1,..., uk 2 by 2 coprime */
     229         497 :   while (k+1 < lg(u))
     230             :   {
     231         287 :     GEN d, c = gel(u,k+1);
     232         287 :     if (is_pm1(c)) { k++; continue; }
     233         539 :     for (i=1; i<=k; i++)
     234             :     {
     235         364 :       GEN ui = gel(u,i);
     236         364 :       if (is_pm1(ui)) continue;
     237         182 :       d = gcdii(c, ui);
     238         182 :       if (d == gen_1) continue;
     239         182 :       c = diviiexact(c, d);
     240         182 :       gel(u,i) = diviiexact(ui, d);
     241         182 :       u = shallowconcat(u, d);
     242             :     }
     243         175 :     gel(u,++k) = c;
     244             :   }
     245         497 :   for (i = k = 1; i < lg(u); i++)
     246         392 :     if (!is_pm1(gel(u,i))) gel(u,k++) = gel(u,i);
     247         105 :   setlg(u, k); return u;
     248             : }
     249             : 
     250             : /*******************************************************************/
     251             : /*                                                                 */
     252             : /*                            ROUND 4                              */
     253             : /*                                                                 */
     254             : /*******************************************************************/
     255             : typedef struct {
     256             :   /* constants */
     257             :   long pisprime; /* -1: unknown, 1: prime,  0: composite */
     258             :   GEN p, f; /* goal: factor f p-adically */
     259             :   long df;
     260             :   GEN pdf; /* p^df = reduced discriminant of f */
     261             :   long mf; /* */
     262             :   GEN psf, pmf; /* stability precision for f, wanted precision for f */
     263             :   long vpsf; /* v_p(p_f) */
     264             :   /* these are updated along the way */
     265             :   GEN phi; /* a p-integer, in Q[X] */
     266             :   GEN phi0; /* a p-integer, in Q[X] from testb2 / testc2, to be composed with
     267             :              * phi when correct precision is known */
     268             :   GEN chi; /* characteristic polynomial of phi (mod psc) in Z[X] */
     269             :   GEN nu; /* irreducible divisor of chi mod p, in Z[X] */
     270             :   GEN invnu; /* numerator ( 1/ Mod(nu, chi) mod pmr ) */
     271             :   GEN Dinvnu;/* denominator ( ... ) */
     272             :   long vDinvnu; /* v_p(Dinvnu) */
     273             :   GEN prc, psc; /* reduced discriminant of chi, stability precision for chi */
     274             :   long vpsc; /* v_p(p_c) */
     275             :   GEN ns, nsf, precns; /* cached Newton sums for nsf and their precision */
     276             : } decomp_t;
     277             : static GEN maxord_i(decomp_t *S, GEN p, GEN f, long mf, GEN w, long flag);
     278             : static GEN dbasis(GEN p, GEN f, long mf, GEN alpha, GEN U);
     279             : static GEN maxord(GEN p,GEN f,long mf);
     280             : static GEN ZX_Dedekind(GEN F, GEN *pg, GEN p);
     281             : 
     282             : /* Warning: data computed for T = ZX_Q_normalize(T0). If S.unscale !=
     283             :  * gen_1, caller must take steps to correct the components if it wishes
     284             :  * to stick to the original T0. Return a vector of p-maximal orders, for
     285             :  * those p s.t p^2 | disc(T) [ = S->dTP ]*/
     286             : static GEN
     287       21826 : get_maxord(nfmaxord_t *S, GEN T0, long flag)
     288             : {
     289             :   VOLATILE GEN P, E, O;
     290             :   VOLATILE long lP, i, k;
     291             : 
     292       21826 :   nfmaxord_check_args(S, T0, flag);
     293       21798 :   P = S->dTP; lP = lg(P);
     294       21798 :   E = S->dTE;
     295       21798 :   O = cgetg(1, t_VEC);
     296       99464 :   for (i=1; i<lP; i++)
     297             :   {
     298             :     VOLATILE pari_sp av;
     299             :     /* includes the silly case where P[i] = -1 */
     300       77666 :     if (E[i] <= 1) { O = shallowconcat(O, gen_1); continue; }
     301       61006 :     av = avma;
     302       61006 :     pari_CATCH(CATCH_ALL) {
     303         105 :       GEN N, u, err = pari_err_last();
     304             :       long l;
     305         105 :       switch(err_get_num(err))
     306             :       {
     307             :         case e_INV:
     308             :         {
     309         105 :           GEN p, x = err_get_compo(err, 2);
     310         105 :           if (typ(x) == t_INTMOD)
     311             :           { /* caught false prime, update factorization */
     312         105 :             p = gcdii(gel(x,1), gel(x,2));
     313         105 :             u = diviiexact(gel(x,1),p);
     314         105 :             if (DEBUGLEVEL) pari_warn(warner,"impossible inverse: %Ps", x);
     315         105 :             gerepileall(av, 2, &p, &u);
     316             : 
     317         105 :             u = get_coprimes(p, u); l = lg(u);
     318             :             /* no small factors, but often a prime power */
     319         105 :             for (k = 1; k < l; k++) (void)Z_isanypower(gel(u,k), &gel(u,k));
     320         105 :             break;
     321             :           }
     322             :           /* fall through */
     323             :         }
     324             :         case e_PRIME: case e_IRREDPOL:
     325             :         { /* we're here because we failed BPSW_isprime(), no point in
     326             :            * reporting a possible counter-example to the BPSW test */
     327           0 :           GEN p = gel(P,i);
     328           0 :           set_avma(av);
     329           0 :           if (DEBUGLEVEL)
     330           0 :             pari_warn(warner,"large composite in nfmaxord:loop(), %Ps", p);
     331           0 :           if (expi(p) < 100) /* factor should require ~20ms for this */
     332           0 :             u = gel(Z_factor(p), 1);
     333             :           else
     334             :           { /* give up, probably not maximal */
     335           0 :             GEN B, g, k = ZX_Dedekind(S->T, &g, p);
     336           0 :             k = FpX_normalize(k, p);
     337           0 :             B = dbasis(p, S->T, E[i], NULL, FpX_div(S->T,k,p));
     338           0 :             O = shallowconcat(O, mkvec(B));
     339           0 :             pari_CATCH_reset(); continue;
     340             :           }
     341           0 :           break;
     342             :         }
     343           0 :         default: pari_err(0, err);
     344             :           return NULL;/*LCOV_EXCL_LINE*/
     345             :       }
     346         105 :       l = lg(u);
     347         105 :       gel(P,i) = gel(u,1);
     348         105 :       P = shallowconcat(P, vecslice(u, 2, l-1));
     349         105 :       av = avma;
     350         105 :       N = S->dT; E[i] = Z_pvalrem(N, gel(P,i), &N);
     351         105 :       for (k=lP, lP=lg(P); k < lP; k++) E[k] = Z_pvalrem(N, gel(P,k), &N);
     352       61111 :     } pari_RETRY {
     353       61111 :       if (DEBUGLEVEL>2) err_printf("Treating p^k = %Ps^%ld\n",P[i],E[i]);
     354       61111 :       O = shallowconcat(O, mkvec( maxord(gel(P,i),S->T,E[i]) ));
     355       61006 :     } pari_ENDCATCH;
     356             :   }
     357       21798 :   S->dTP = P; return O;
     358             : }
     359             : 
     360             : /* M a QM, return denominator of diagonal. All denominators are powers of
     361             :  * a given integer */
     362             : static GEN
     363       11432 : diag_denom(GEN M)
     364             : {
     365       11432 :   GEN d = gen_1;
     366       11432 :   long j, l = lg(M);
     367      121876 :   for (j=1; j<l; j++)
     368             :   {
     369      110444 :     GEN t = gcoeff(M,j,j);
     370      110444 :     if (typ(t) == t_INT) continue;
     371       28631 :     t = gel(t,2);
     372       28631 :     if (abscmpii(t,d) > 0) d = t;
     373             :   }
     374       11432 :   return d;
     375             : }
     376             : static void
     377       18655 : setPE(GEN D, GEN P, GEN *pP, GEN *pE)
     378             : {
     379       18655 :   long k, j, l = lg(P);
     380             :   GEN P2, E2;
     381       18655 :   *pP = P2 = cgetg(l, t_COL);
     382       18655 :   *pE = E2 = cgetg(l, t_VECSMALL);
     383       57226 :   for (k = j = 1; j < l; j++)
     384             :   {
     385       38571 :     long v = Z_pvalrem(D, gel(P,j), &D);
     386       38571 :     if (v) { gel(P2,k) = gel(P,j); E2[k] = v; k++; }
     387             :   }
     388       18655 :   setlg(P2, k);
     389       18655 :   setlg(E2, k);
     390       18655 : }
     391             : void
     392       18641 : nfmaxord(nfmaxord_t *S, GEN T0, long flag)
     393             : {
     394       18641 :   GEN O = get_maxord(S, T0, flag);
     395       18634 :   GEN f = S->T, P = S->dTP, a = NULL, da = NULL;
     396       18634 :   long n = degpol(f), lP = lg(P), i, j, k;
     397       18634 :   int centered = 0;
     398       18634 :   pari_sp av = avma;
     399             :   /* r1 & basden not initialized here */
     400       18634 :   S->r1 = -1;
     401       18634 :   S->basden = NULL;
     402       57156 :   for (i=1; i<lP; i++)
     403             :   {
     404       38522 :     GEN M, db, b = gel(O,i);
     405       38522 :     if (b == gen_1) continue;
     406       11432 :     db = diag_denom(b);
     407       11432 :     if (db == gen_1) continue;
     408             : 
     409             :     /* db = denom(b), (da,db) = 1. Compute da Im(b) + db Im(a) */
     410       11432 :     b = Q_muli_to_int(b,db);
     411       11432 :     if (!da) { da = db; a = b; }
     412             :     else
     413             :     { /* optimization: easy as long as both matrix are diagonal */
     414        7155 :       j=2; while (j<=n && fnz(gel(a,j),j) && fnz(gel(b,j),j)) j++;
     415        7155 :       k = j-1; M = cgetg(2*n-k+1,t_MAT);
     416       39469 :       for (j=1; j<=k; j++)
     417             :       {
     418       32314 :         gel(M,j) = gel(a,j);
     419       32314 :         gcoeff(M,j,j) = mulii(gcoeff(a,j,j),gcoeff(b,j,j));
     420             :       }
     421             :       /* could reduce mod M(j,j) but not worth it: usually close to da*db */
     422        7155 :       for (  ; j<=n;     j++) gel(M,j) = ZC_Z_mul(gel(a,j), db);
     423        7155 :       for (  ; j<=2*n-k; j++) gel(M,j) = ZC_Z_mul(gel(b,j+k-n), da);
     424        7155 :       da = mulii(da,db);
     425        7155 :       a = ZM_hnfmodall_i(M, da, hnf_MODID|hnf_CENTER);
     426        7155 :       gerepileall(av, 2, &a, &da);
     427        7155 :       centered = 1;
     428             :     }
     429             :   }
     430       18634 :   if (da)
     431             :   {
     432        4277 :     GEN index = diviiexact(da, gcoeff(a,1,1));
     433        4277 :     for (j=2; j<=n; j++) index = mulii(index, diviiexact(da, gcoeff(a,j,j)));
     434        4277 :     if (!centered) a = ZM_hnfcenter(a);
     435        4277 :     a = RgM_Rg_div(a, da);
     436        4277 :     S->index = index;
     437        4277 :     S->dK = diviiexact(S->dT, sqri(index));
     438             :   }
     439             :   else
     440             :   {
     441       14357 :     S->index = gen_1;
     442       14357 :     S->dK = S->dT;
     443       14357 :     a = matid(n);
     444             :   }
     445       18634 :   setPE(S->dK, P, &S->dKP, &S->dKE);
     446       18634 :   S->basis = RgM_to_RgXV(a, varn(f));
     447       18634 : }
     448             : GEN
     449          70 : nfbasis(GEN x, GEN *pdK)
     450             : {
     451          70 :   pari_sp av = avma;
     452             :   nfmaxord_t S;
     453             :   GEN B;
     454          70 :   nfmaxord(&S, x, 0);
     455          70 :   B = RgXV_unscale(S.basis, S.unscale);
     456          70 :   if (pdK)  *pdK = S.dK;
     457          70 :   gerepileall(av, pdK? 2: 1, &B, pdK); return B;
     458             : }
     459             : /* field discriminant: faster than nfmaxord, use local data only */
     460             : static GEN
     461        3185 : maxord_disc(nfmaxord_t *S, GEN x)
     462             : {
     463        3185 :   GEN O = get_maxord(S, x, 0), I = gen_1;
     464        3164 :   long n = degpol(S->T), lP = lg(O), i, j;
     465       42308 :   for (i = 1; i < lP; i++)
     466             :   {
     467       39144 :     GEN b = gel(O,i);
     468       39144 :     if (b == gen_1) continue;
     469      399595 :     for (j = 1; j <= n; j++)
     470             :     {
     471      365736 :       GEN c = gcoeff(b,j,j);
     472      365736 :       if (typ(c) == t_FRAC) I = mulii(I, gel(c,2)) ;
     473             :     }
     474             :   }
     475        3164 :   return diviiexact(S->dT, sqri(I));
     476             : }
     477             : GEN
     478        3171 : nfdisc(GEN x)
     479             : {
     480        3171 :   pari_sp av = avma;
     481             :   nfmaxord_t S;
     482        3171 :   return gerepileuptoint(av, maxord_disc(&S, x));
     483             : }
     484             : GEN
     485          21 : nfdiscfactors(GEN x)
     486             : {
     487          21 :   pari_sp av = avma;
     488          21 :   GEN E, P, D, nf = checknf_i(x);
     489          21 :   if (nf)
     490             :   {
     491           7 :     D = nf_get_disc(nf);
     492           7 :     P = nf_get_ramified_primes(nf);
     493             :   }
     494             :   else
     495             :   {
     496             :     nfmaxord_t S;
     497          14 :     D = maxord_disc(&S, x);
     498          14 :     P = S.dTP;
     499             :   }
     500          21 :   setPE(D, P, &P, &E); settyp(P, t_COL);
     501          21 :   return gerepilecopy(av, mkvec2(D, mkmat2(P, zc_to_ZC(E))));
     502             : }
     503             : 
     504             : static ulong
     505      121051 : Flx_checkdeflate(GEN x)
     506             : {
     507      121051 :   ulong d = 0, i, lx = (ulong)lg(x);
     508      267744 :   for (i=3; i<lx; i++)
     509      230365 :     if (x[i]) { d = ugcd(d,i-2); if (d == 1) break; }
     510      121051 :   return d;
     511             : }
     512             : 
     513             : /* product of (monic) irreducible factors of f over Fp[X]
     514             :  * Assume f reduced mod p, otherwise valuation at x may be wrong */
     515             : static GEN
     516      121051 : Flx_radical(GEN f, ulong p)
     517             : {
     518      121051 :   long v0 = Flx_valrem(f, &f);
     519             :   ulong du, d, e;
     520             :   GEN u;
     521             : 
     522      121051 :   d = Flx_checkdeflate(f);
     523      121051 :   if (!d) return v0? polx_Flx(f[1]): pol1_Flx(f[1]);
     524      103029 :   if (u_lvalrem(d,p, &e)) f = Flx_deflate(f, d/e); /* f(x^p^i) -> f(x) */
     525      103029 :   u = Flx_gcd(f, Flx_deriv(f, p), p); /* (f,f') */
     526      103023 :   du = degpol(u);
     527      103023 :   if (du)
     528             :   {
     529       70256 :     if (du == (ulong)degpol(f))
     530           0 :       f = Flx_radical(Flx_deflate(f,p), p);
     531             :     else
     532             :     {
     533       70256 :       u = Flx_normalize(u, p);
     534       70256 :       f = Flx_div(f, u, p);
     535       70256 :       if (p <= du)
     536             :       {
     537        8029 :         GEN w = (degpol(f) >= degpol(u))? Flx_rem(f, u, p): f;
     538        8029 :         w = Flxq_powu(w, du, u, p);
     539        8029 :         w = Flx_div(u, Flx_gcd(w,u,p), p); /* u / gcd(u, v^(deg u-1)) */
     540        8029 :         f = Flx_mul(f, Flx_radical(Flx_deflate(w,p), p), p);
     541             :       }
     542             :     }
     543             :   }
     544      103023 :   if (v0) f = Flx_shift(f, 1);
     545      103023 :   return f;
     546             : }
     547             : /* Assume f reduced mod p, otherwise valuation at x may be wrong */
     548             : static GEN
     549        3365 : FpX_radical(GEN f, GEN p)
     550             : {
     551             :   GEN u;
     552             :   long v0;
     553        3365 :   if (lgefint(p) == 3)
     554             :   {
     555         572 :     ulong q = p[2];
     556         572 :     return Flx_to_ZX( Flx_radical(ZX_to_Flx(f, q), q) );
     557             :   }
     558        2793 :   v0 = ZX_valrem(f, &f);
     559        2793 :   u = FpX_gcd(f,FpX_deriv(f, p), p);
     560        2694 :   if (degpol(u)) f = FpX_div(f, u, p);
     561        2694 :   if (v0) f = RgX_shift(f, 1);
     562        2694 :   return f;
     563             : }
     564             : /* f / a */
     565             : static GEN
     566      112450 : zx_z_div(GEN f, ulong a)
     567             : {
     568      112450 :   long i, l = lg(f);
     569      112450 :   GEN g = cgetg(l, t_VECSMALL);
     570      112450 :   g[1] = f[1];
     571      112450 :   for (i = 2; i < l; i++) g[i] = f[i] / a;
     572      112450 :   return g;
     573             : }
     574             : /* Dedekind criterion; return k = gcd(g,h, (f-gh)/p), where
     575             :  *   f = \prod f_i^e_i, g = \prod f_i, h = \prod f_i^{e_i-1}
     576             :  * k = 1 iff Z[X]/(f) is p-maximal */
     577             : static GEN
     578      115815 : ZX_Dedekind(GEN F, GEN *pg, GEN p)
     579             : {
     580             :   GEN k, h, g, f, f2;
     581      115815 :   ulong q = p[2];
     582      115815 :   if (lgefint(p) == 3 && q < (1UL << BITS_IN_HALFULONG))
     583      112450 :   {
     584      112450 :     ulong q2 = q*q;
     585      112450 :     f2 = ZX_to_Flx(F, q2);
     586      112450 :     f = Flx_red(f2, q);
     587      112450 :     g = Flx_radical(f, q);
     588      112450 :     h = Flx_div(f, g, q);
     589      112450 :     k = zx_z_div(Flx_sub(f2, Flx_mul(g,h,q2), q2), q);
     590      112450 :     k = Flx_gcd(k, Flx_gcd(g,h,q), q);
     591      112450 :     k = Flx_to_ZX(k);
     592      112450 :     g = Flx_to_ZX(g);
     593             :   }
     594             :   else
     595             :   {
     596        3365 :     f2 = FpX_red(F, sqri(p));
     597        3365 :     f = FpX_red(f2, p);
     598        3365 :     g = FpX_radical(f, p);
     599        3260 :     h = FpX_div(f, g, p);
     600        3260 :     k = ZX_Z_divexact(ZX_sub(f2, ZX_mul(g,h)), p);
     601        3260 :     k = FpX_gcd(FpX_red(k, p), FpX_gcd(g,h,p), p);
     602             :   }
     603      115710 :   *pg = g; return k;
     604             : }
     605             : 
     606             : /* p-maximal order of Z[x]/f; mf = v_p(Disc(f)) or < 0 [unknown].
     607             :  * Return gen_1 if p-maximal */
     608             : static GEN
     609      115815 : maxord(GEN p, GEN f, long mf)
     610             : {
     611      115815 :   const pari_sp av = avma;
     612      115815 :   GEN res, g, k = ZX_Dedekind(f, &g, p);
     613      115710 :   long dk = degpol(k);
     614      115710 :   if (DEBUGLEVEL>2) err_printf("  ZX_Dedekind: gcd has degree %ld\n", dk);
     615      115710 :   if (!dk) { set_avma(av); return gen_1; }
     616       76901 :   if (mf < 0) mf = ZpX_disc_val(f, p);
     617       76901 :   k = FpX_normalize(k, p);
     618       76901 :   if (2*dk >= mf-1)
     619       38668 :     res = dbasis(p, f, mf, NULL, FpX_div(f,k,p));
     620             :   else
     621             :   {
     622             :     GEN w, F1, F2;
     623             :     decomp_t S;
     624       38233 :     F1 = FpX_factor(k,p);
     625       38233 :     F2 = FpX_factor(FpX_div(g,k,p),p);
     626       38233 :     w = merge_sort_uniq(gel(F1,1),gel(F2,1),(void*)cmpii,&gen_cmp_RgX);
     627       38233 :     res = maxord_i(&S, p, f, mf, w, 0);
     628             :   }
     629       76901 :   return gerepilecopy(av,res);
     630             : }
     631             : /* T monic separable ZX, p prime */
     632             : GEN
     633           0 : ZpX_primedec(GEN T, GEN p)
     634             : {
     635           0 :   const pari_sp av = avma;
     636           0 :   GEN w, F1, F2, res, g, k = ZX_Dedekind(T, &g, p);
     637             :   decomp_t S;
     638           0 :   if (!degpol(k)) return zm_to_ZM(FpX_degfact(T, p));
     639           0 :   k = FpX_normalize(k, p);
     640           0 :   F1 = FpX_factor(k,p);
     641           0 :   F2 = FpX_factor(FpX_div(g,k,p),p);
     642           0 :   w = merge_sort_uniq(gel(F1,1),gel(F2,1),(void*)cmpii,&gen_cmp_RgX);
     643           0 :   res = maxord_i(&S, p, T, ZpX_disc_val(T, p), w, -1);
     644           0 :   if (!res)
     645             :   {
     646           0 :     long f = degpol(S.nu), e = degpol(T) / f;
     647           0 :     set_avma(av); retmkmat2(mkcols(f), mkcols(e));
     648             :   }
     649           0 :   return gerepilecopy(av,res);
     650             : }
     651             : 
     652             : static GEN
     653      935272 : Zlx_sylvester_echelon(GEN f1, GEN f2, long early_abort, ulong p, ulong pm)
     654             : {
     655      935272 :   long j, n = degpol(f1);
     656      935269 :   GEN h, a = cgetg(n+1,t_MAT);
     657      935264 :   f1 = Flx_get_red(f1, pm);
     658      935161 :   h = Flx_rem(f2,f1,pm);
     659     4050671 :   for (j=1;; j++)
     660             :   {
     661     7166070 :     gel(a,j) = Flx_to_Flv(h, n);
     662     4050688 :     if (j == n) break;
     663     3115441 :     h = Flx_rem(Flx_shift(h, 1), f1, pm);
     664             :   }
     665      935247 :   return zlm_echelon(a, early_abort, p, pm);
     666             : }
     667             : /* Sylvester's matrix, mod p^m (assumes f1 monic). If early_abort
     668             :  * is set, return NULL if one pivot is 0 mod p^m */
     669             : static GEN
     670       15626 : ZpX_sylvester_echelon(GEN f1, GEN f2, long early_abort, GEN p, GEN pm)
     671             : {
     672       15626 :   long j, n = degpol(f1);
     673       15626 :   GEN h, a = cgetg(n+1,t_MAT);
     674       15626 :   h = FpXQ_red(f2,f1,pm);
     675      155936 :   for (j=1;; j++)
     676             :   {
     677      296246 :     gel(a,j) = RgX_to_RgC(h, n);
     678      155936 :     if (j == n) break;
     679      140310 :     h = FpX_rem(RgX_shift_shallow(h, 1), f1, pm);
     680             :   }
     681       15626 :   return ZpM_echelon(a, early_abort, p, pm);
     682             : }
     683             : 
     684             : /* polynomial gcd mod p^m (assumes f1 monic). Return a QpX ! */
     685             : static GEN
     686       21554 : Zlx_gcd(GEN f1, GEN f2, ulong p, ulong pm)
     687             : {
     688       21554 :   pari_sp av = avma;
     689       21554 :   GEN a = Zlx_sylvester_echelon(f1,f2,0,p,pm);
     690       21554 :   long c, l = lg(a), sv = f1[1];
     691      126579 :   for (c = 1; c < l; c++)
     692             :   {
     693      126579 :     ulong t = ucoeff(a,c,c);
     694      126579 :     if (t)
     695             :     {
     696       21554 :       a = Flx_to_ZX(Flv_to_Flx(gel(a,c), sv));
     697       21554 :       if (t == 1) return gerepilecopy(av, a);
     698        2912 :       return gerepileupto(av, RgX_Rg_div(a, utoipos(t)));
     699             :     }
     700             :   }
     701           0 :   set_avma(av);
     702           0 :   a = cgetg(2,t_POL); a[1] = sv; return a;
     703             : }
     704             : GEN
     705       27597 : ZpX_gcd(GEN f1, GEN f2, GEN p, GEN pm)
     706             : {
     707       27597 :   pari_sp av = avma;
     708             :   GEN a;
     709             :   long c, l, v;
     710       27597 :   if (lgefint(pm) == 3)
     711             :   {
     712       21554 :     ulong q = pm[2];
     713       21554 :     return Zlx_gcd(ZX_to_Flx(f1, q), ZX_to_Flx(f2,q), p[2], q);
     714             :   }
     715        6043 :   a = ZpX_sylvester_echelon(f1,f2,0,p,pm);
     716        6043 :   l = lg(a); v = varn(f1);
     717       41036 :   for (c = 1; c < l; c++)
     718             :   {
     719       41036 :     GEN t = gcoeff(a,c,c);
     720       41036 :     if (signe(t))
     721             :     {
     722        6043 :       a = RgV_to_RgX(gel(a,c), v);
     723        6043 :       if (equali1(t)) return gerepilecopy(av, a);
     724        1715 :       return gerepileupto(av, RgX_Rg_div(a, t));
     725             :     }
     726             :   }
     727           0 :   set_avma(av); return pol_0(v);
     728             : }
     729             : 
     730             : /* Return m > 0, such that p^m ~ 2^16 for initial value of m; p > 1 */
     731             : static long
     732      865285 : init_m(GEN p)
     733             : {
     734      865285 :   if (lgefint(p) > 3) return 1;
     735      865086 :   return (long)(16 / log2(p[2]));
     736             : }
     737             : 
     738             : /* reduced resultant mod p^m (assumes x monic) */
     739             : GEN
     740       92491 : ZpX_reduced_resultant(GEN x, GEN y, GEN p, GEN pm)
     741             : {
     742       92491 :   pari_sp av = avma;
     743             :   GEN z;
     744       92491 :   if (lgefint(pm) == 3)
     745             :   {
     746       87070 :     ulong q = pm[2];
     747       87070 :     z = Zlx_sylvester_echelon(ZX_to_Flx(x,q), ZX_to_Flx(y,q),0,p[2],q);
     748       87070 :     if (lg(z) > 1)
     749             :     {
     750       87070 :       ulong c = ucoeff(z,1,1);
     751       87070 :       if (c) { set_avma(av); return utoipos(c); }
     752             :     }
     753             :   }
     754             :   else
     755             :   {
     756        5421 :     z = ZpX_sylvester_echelon(x,y,0,p,pm);
     757        5421 :     if (lg(z) > 1)
     758             :     {
     759        5421 :       GEN c = gcoeff(z,1,1);
     760        5421 :       if (signe(c)) return gerepileuptoint(av, c);
     761             :     }
     762             :   }
     763       36940 :   set_avma(av); return gen_0;
     764             : }
     765             : /* Assume Res(f,g) divides p^M. Return Res(f, g), using dynamic p-adic
     766             :  * precision (until result is non-zero or p^M). */
     767             : GEN
     768       63651 : ZpX_reduced_resultant_fast(GEN f, GEN g, GEN p, long M)
     769             : {
     770       63651 :   GEN R, q = NULL;
     771             :   long m;
     772       63651 :   m = init_m(p); if (m < 1) m = 1;
     773       28840 :   for(;; m <<= 1) {
     774      121331 :     if (M < 2*m) break;
     775       44104 :     q = q? sqri(q): powiu(p, m); /* p^m */
     776       44104 :     R = ZpX_reduced_resultant(f,g, p, q); if (signe(R)) return R;
     777             :   }
     778       48387 :   q = powiu(p, M);
     779       48387 :   R = ZpX_reduced_resultant(f,g, p, q); return signe(R)? R: q;
     780             : }
     781             : 
     782             : /* v_p(Res(x,y) mod p^m), assumes (lc(x),p) = 1 */
     783             : static long
     784      830851 : ZpX_resultant_val_i(GEN x, GEN y, GEN p, GEN pm)
     785             : {
     786      830851 :   pari_sp av = avma;
     787             :   GEN z;
     788             :   long i, l, v;
     789      830851 :   if (lgefint(pm) == 3)
     790             :   {
     791      826689 :     ulong q = pm[2], pp = p[2];
     792      826689 :     z = Zlx_sylvester_echelon(ZX_to_Flx(x,q), ZX_to_Flx(y,q), 1, pp, q);
     793      826692 :     if (!z) return gc_long(av,-1); /* failure */
     794      755392 :     v = 0; l = lg(z);
     795      755392 :     for (i = 1; i < l; i++) v += u_lval(ucoeff(z,i,i), pp);
     796             :   }
     797             :   else
     798             :   {
     799        4162 :     z = ZpX_sylvester_echelon(x, y, 1, p, pm);
     800        4162 :     if (!z) return gc_long(av,-1); /* failure */
     801        3599 :     v = 0; l = lg(z);
     802        3599 :     for (i = 1; i < l; i++) v += Z_pval(gcoeff(z,i,i), p);
     803             :   }
     804      758993 :   return v;
     805             : }
     806             : 
     807             : /* assume (lc(f),p) = 1; no assumption on g */
     808             : long
     809      801618 : ZpX_resultant_val(GEN f, GEN g, GEN p, long M)
     810             : {
     811      801618 :   pari_sp av = avma;
     812      801618 :   GEN q = NULL;
     813             :   long v, m;
     814      801618 :   m = init_m(p); if (m < 2) m = 2;
     815       29207 :   for(;; m <<= 1) {
     816      860046 :     if (m > M) m = M;
     817      830839 :     q = q? sqri(q): powiu(p, m); /* p^m */
     818      830853 :     v = ZpX_resultant_val_i(f,g, p, q); if (v >= 0) return gc_long(av,v);
     819       71863 :     if (m == M) return gc_long(av,M);
     820             :   }
     821             : }
     822             : 
     823             : /* assume f separable and (lc(f),p) = 1 */
     824             : long
     825       32520 : ZpX_disc_val(GEN f, GEN p)
     826             : {
     827       32520 :   pari_sp av = avma;
     828             :   long v;
     829       32520 :   if (degpol(f) == 1) return 0;
     830       32520 :   v = ZpX_resultant_val(f, ZX_deriv(f), p, LONG_MAX);
     831       32520 :   return gc_long(av,v);
     832             : }
     833             : 
     834             : /* *e a ZX, *d, *z in Z, *d = p^(*vd). Simplify e / d by cancelling a
     835             :  * common factor p^v; if z!=NULL, update it by cancelling the same power of p */
     836             : static void
     837      588750 : update_den(GEN p, GEN *e, GEN *d, long *vd, GEN *z)
     838             : {
     839             :   GEN newe;
     840      588750 :   long ve = ZX_pvalrem(*e, p, &newe);
     841      588750 :   if (ve) {
     842             :     GEN newd;
     843      332342 :     long v = minss(*vd, ve);
     844      332342 :     if (v) {
     845      332342 :       if (v == *vd)
     846             :       { /* rare, denominator cancelled */
     847       38593 :         if (ve != v) newe = ZX_Z_mul(newe, powiu(p, ve - v));
     848       38593 :         newd = gen_1;
     849       38593 :         *vd = 0;
     850       38593 :         if (z) *z =diviiexact(*z, powiu(p, v));
     851             :       }
     852             :       else
     853             :       { /* v = ve < vd, generic case */
     854      293749 :         GEN q = powiu(p, v);
     855      293749 :         newd = diviiexact(*d, q);
     856      293749 :         *vd -= v;
     857      293749 :         if (z) *z = diviiexact(*z, q);
     858             :       }
     859      332342 :       *e = newe;
     860      332342 :       *d = newd;
     861             :     }
     862             :   }
     863      588750 : }
     864             : 
     865             : /* return denominator, a power of p */
     866             : static GEN
     867      386377 : QpX_denom(GEN x)
     868             : {
     869      386377 :   long i, l = lg(x);
     870      386377 :   GEN maxd = gen_1;
     871     1965003 :   for (i=2; i<l; i++)
     872             :   {
     873     1578626 :     GEN d = gel(x,i);
     874     1578626 :     if (typ(d) == t_FRAC && cmpii(gel(d,2), maxd) > 0) maxd = gel(d,2);
     875             :   }
     876      386377 :   return maxd;
     877             : }
     878             : static GEN
     879       54704 : QpXV_denom(GEN x)
     880             : {
     881       54704 :   long l = lg(x), i;
     882       54704 :   GEN maxd = gen_1;
     883      277975 :   for (i = 1; i < l; i++)
     884             :   {
     885      223271 :     GEN d = QpX_denom(gel(x,i));
     886      223271 :     if (cmpii(d, maxd) > 0) maxd = d;
     887             :   }
     888       54704 :   return maxd;
     889             : }
     890             : 
     891             : static GEN
     892      163106 : QpX_remove_denom(GEN x, GEN p, GEN *pdx, long *pv)
     893             : {
     894      163106 :   *pdx = QpX_denom(x);
     895      163106 :   if (*pdx == gen_1) { *pv = 0; *pdx = NULL; }
     896             :   else {
     897      122798 :     x = Q_muli_to_int(x,*pdx);
     898      122798 :     *pv = Z_pval(*pdx, p);
     899             :   }
     900      163106 :   return x;
     901             : }
     902             : 
     903             : /* p^v * f o g mod (T,q). q = p^vq  */
     904             : static GEN
     905       23565 : compmod(GEN p, GEN f, GEN g, GEN T, GEN q, long v)
     906             : {
     907       23565 :   GEN D = NULL, z, df, dg, qD;
     908       23565 :   long vD = 0, vdf, vdg;
     909             : 
     910       23565 :   f = QpX_remove_denom(f, p, &df, &vdf);
     911       23565 :   if (typ(g) == t_VEC) /* [num,den,v_p(den)] */
     912           0 :   { vdg = itos(gel(g,3)); dg = gel(g,2); g = gel(g,1); }
     913             :   else
     914       23565 :     g = QpX_remove_denom(g, p, &dg, &vdg);
     915       23565 :   if (df) { D = df; vD = vdf; }
     916       23565 :   if (dg) {
     917        4242 :     long degf = degpol(f);
     918        4242 :     D = mul_content(D, powiu(dg, degf));
     919        4242 :     vD += degf * vdg;
     920             :   }
     921       23565 :   qD = D ? mulii(q, D): q;
     922       23565 :   if (dg) f = FpX_rescale(f, dg, qD);
     923       23565 :   z = FpX_FpXQ_eval(f, g, T, qD);
     924       23565 :   if (!D) {
     925           0 :     if (v) {
     926           0 :       if (v > 0)
     927           0 :         z = ZX_Z_mul(z, powiu(p, v));
     928             :       else
     929           0 :         z = RgX_Rg_div(z, powiu(p, -v));
     930             :     }
     931           0 :     return z;
     932             :   }
     933       23565 :   update_den(p, &z, &D, &vD, NULL);
     934       23565 :   qD = mulii(D,q);
     935       23565 :   if (v) vD -= v;
     936       23565 :   z = FpX_center_i(z, qD, shifti(qD,-1));
     937       23565 :   if (vD > 0)
     938       23565 :     z = RgX_Rg_div(z, powiu(p, vD));
     939           0 :   else if (vD < 0)
     940           0 :     z = ZX_Z_mul(z, powiu(p, -vD));
     941       23565 :   return z;
     942             : }
     943             : 
     944             : /* fast implementation of ZM_hnfmodid(M, D) / D, D = p^k */
     945             : static GEN
     946       38233 : ZpM_hnfmodid(GEN M, GEN p, GEN D)
     947             : {
     948       38233 :   long i, l = lg(M);
     949       38233 :   M = RgM_Rg_div(ZpM_echelon(M,0,p,D), D);
     950      308650 :   for (i = 1; i < l; i++)
     951      270417 :     if (gequal0(gcoeff(M,i,i))) gcoeff(M,i,i) = gen_1;
     952       38233 :   return M;
     953             : }
     954             : 
     955             : /* Return Z-basis for Z[a] + U(a)/p Z[a] in Z[t]/(f), mf = v_p(disc f), U
     956             :  * a ZX. Special cases: a = t is coded as NULL, U = 0 is coded as NULL */
     957             : static GEN
     958       49549 : dbasis(GEN p, GEN f, long mf, GEN a, GEN U)
     959             : {
     960       49549 :   long n = degpol(f), i, dU;
     961             :   GEN b, h;
     962             : 
     963       49549 :   if (n == 1) return matid(1);
     964       49549 :   if (a && gequalX(a)) a = NULL;
     965       49549 :   if (DEBUGLEVEL>5)
     966             :   {
     967           0 :     err_printf("  entering Dedekind Basis with parameters p=%Ps\n",p);
     968           0 :     err_printf("  f = %Ps,\n  a = %Ps\n",f, a? a: pol_x(varn(f)));
     969             :   }
     970       49549 :   if (a)
     971             :   {
     972       10881 :     GEN pd = powiu(p, mf >> 1);
     973       10881 :     GEN da, pdp = mulii(pd,p), D = pdp;
     974             :     long vda;
     975       10881 :     dU = U ? degpol(U): 0;
     976       10881 :     b = cgetg(n+1, t_MAT);
     977       10881 :     h = scalarpol(pd, varn(f));
     978       10881 :     a = QpX_remove_denom(a, p, &da, &vda);
     979       10881 :     if (da) D = mulii(D, da);
     980       10881 :     gel(b,1) = scalarcol_shallow(pd, n);
     981       47146 :     for (i=2; i<=n; i++)
     982             :     {
     983       36265 :       if (i == dU+1)
     984           0 :         h = compmod(p, U, mkvec3(a,da,stoi(vda)), f, pdp, (mf>>1) - 1);
     985             :       else
     986             :       {
     987       36265 :         h = FpXQ_mul(h, a, f, D);
     988       36265 :         if (da) h = ZX_Z_divexact(h, da);
     989             :       }
     990       36265 :       gel(b,i) = RgX_to_RgC(h,n);
     991             :     }
     992       10881 :     return ZpM_hnfmodid(b, p, pd);
     993             :   }
     994             :   else
     995             :   {
     996       38668 :     if (!U) return matid(n);
     997       38668 :     dU = degpol(U);
     998       38668 :     if (dU == n) return matid(n);
     999       38668 :     U = FpX_normalize(U, p);
    1000       38668 :     b = cgetg(n+1, t_MAT);
    1001       38668 :     for (i = 1; i <= dU; i++) gel(b,i) = vec_ei(n, i);
    1002       38668 :     h = RgX_Rg_div(U, p);
    1003       51513 :     for ( ; i <= n; i++)
    1004             :     {
    1005       51513 :       gel(b, i) = RgX_to_RgC(h,n);
    1006       51513 :       if (i == n) break;
    1007       12845 :       h = RgX_shift_shallow(h,1);
    1008             :     }
    1009       38668 :     return b;
    1010             :   }
    1011             : }
    1012             : 
    1013             : static GEN
    1014       54704 : get_partial_order_as_pols(GEN p, GEN f)
    1015             : {
    1016       54704 :   GEN O = maxord(p, f, -1);
    1017       54704 :   long v = varn(f);
    1018       54704 :   return O == gen_1? pol_x_powers(degpol(f), v): RgM_to_RgXV(O, v);
    1019             : }
    1020             : 
    1021             : static long
    1022        1099 : p_is_prime(decomp_t *S)
    1023             : {
    1024        1099 :   if (S->pisprime < 0) S->pisprime = BPSW_psp(S->p);
    1025        1099 :   return S->pisprime;
    1026             : }
    1027             : static GEN ZpX_monic_factor_squarefree(GEN f, GEN p, long prec);
    1028             : 
    1029             : /* if flag = 0, maximal order, else factorization to precision r = flag */
    1030             : static GEN
    1031       27597 : Decomp(decomp_t *S, long flag)
    1032             : {
    1033       27597 :   pari_sp av = avma;
    1034             :   GEN fred, pr2, pr, pk, ph2, ph, b1, b2, a, e, de, f1, f2, dt, th, chip;
    1035       27597 :   GEN p = S->p;
    1036       27597 :   long vde, vdt, k, r = maxss(flag, 2*S->df + 1);
    1037             : 
    1038       27597 :   if (DEBUGLEVEL>5) err_printf("  entering Decomp: %Ps^%ld\n  f = %Ps\n",
    1039             :                                p, r, S->f);
    1040       27597 :   else if (DEBUGLEVEL>2) err_printf("  entering Decomp\n");
    1041       27597 :   chip = FpX_red(S->chi, p);
    1042       27597 :   if (!FpX_valrem(chip, S->nu, p, &b1))
    1043             :   {
    1044           0 :     if (!p_is_prime(S)) pari_err_PRIME("Decomp",p);
    1045           0 :     pari_err_BUG("Decomp (not a factor)");
    1046             :   }
    1047       27597 :   b2 = FpX_div(chip, b1, p);
    1048       27597 :   a = FpX_mul(FpXQ_inv(b2, b1, p), b2, p);
    1049             :   /* E = e / de, e in Z[X], de in Z,  E = a(phi) mod (f, p) */
    1050       27597 :   th = QpX_remove_denom(S->phi, p, &dt, &vdt);
    1051       27597 :   if (dt)
    1052             :   {
    1053       11130 :     long dega = degpol(a);
    1054       11130 :     vde = dega * vdt;
    1055       11130 :     de = powiu(dt, dega);
    1056       11130 :     pr = mulii(p, de);
    1057       11130 :     a = FpX_rescale(a, dt, pr);
    1058             :   }
    1059             :   else
    1060             :   {
    1061       16467 :     vde = 0;
    1062       16467 :     de = gen_1;
    1063       16467 :     pr = p;
    1064             :   }
    1065       27597 :   e = FpX_FpXQ_eval(a, th, S->f, pr);
    1066       27597 :   update_den(p, &e, &de, &vde, NULL);
    1067             : 
    1068       27597 :   pk = p; k = 1;
    1069             :   /* E, (1 - E) tend to orthogonal idempotents in Zp[X]/(f) */
    1070      174793 :   while (k < r + vde)
    1071             :   { /* E <-- E^2(3-2E) mod p^2k, with E = e/de */
    1072             :     GEN D;
    1073      119599 :     pk = sqri(pk); k <<= 1;
    1074      119599 :     e = ZX_mul(ZX_sqr(e), Z_ZX_sub(mului(3,de), gmul2n(e,1)));
    1075      119599 :     de= mulii(de, sqri(de));
    1076      119599 :     vde *= 3;
    1077      119599 :     D = mulii(pk, de);
    1078      119599 :     e = FpX_rem(e, centermod(S->f, D), D); /* e/de defined mod pk */
    1079      119599 :     update_den(p, &e, &de, &vde, NULL);
    1080             :   }
    1081             :   /* required precision of the factors */
    1082       27597 :   pr = powiu(p, r); pr2 = shifti(pr, -1);
    1083       27597 :   ph = mulii(de,pr);ph2 = shifti(ph, -1);
    1084       27597 :   fred = FpX_center_i(FpX_red(S->f, ph), ph, ph2);
    1085       27597 :   e    = FpX_center_i(FpX_red(e, ph), ph, ph2);
    1086             : 
    1087       27597 :   f1 = ZpX_gcd(fred, Z_ZX_sub(de, e), p, ph); /* p-adic gcd(f, 1-e) */
    1088       27597 :   fred = FpX_center_i(fred, pr, pr2);
    1089       27597 :   f1   = FpX_center_i(f1,   pr, pr2);
    1090       27597 :   f2 = FpX_div(fred,f1, pr);
    1091       27597 :   f2 = FpX_center_i(f2, pr, pr2);
    1092             : 
    1093       27597 :   if (DEBUGLEVEL>5)
    1094           0 :     err_printf("  leaving Decomp: f1 = %Ps\nf2 = %Ps\ne = %Ps\nde= %Ps\n", f1,f2,e,de);
    1095             : 
    1096       27597 :   if (flag < 0)
    1097             :   {
    1098           0 :     GEN m = vconcat(ZpX_primedec(f1, p), ZpX_primedec(f2, p));
    1099           0 :     return sort_factor(m, (void*)&cmpii, &cmp_nodata);
    1100             :   }
    1101       27597 :   else if (flag)
    1102             :   {
    1103         245 :     gerepileall(av, 2, &f1, &f2);
    1104         245 :     return shallowconcat(ZpX_monic_factor_squarefree(f1, p, flag),
    1105             :                          ZpX_monic_factor_squarefree(f2, p, flag));
    1106             :   } else {
    1107             :     GEN D, d1, d2, B1, B2, M;
    1108             :     long n, n1, n2, i;
    1109       27352 :     gerepileall(av, 4, &f1, &f2, &e, &de);
    1110       27352 :     D = de;
    1111       27352 :     B1 = get_partial_order_as_pols(p,f1); n1 = lg(B1)-1;
    1112       27352 :     B2 = get_partial_order_as_pols(p,f2); n2 = lg(B2)-1; n = n1+n2;
    1113       27352 :     d1 = QpXV_denom(B1);
    1114       27352 :     d2 = QpXV_denom(B2); if (cmpii(d1, d2) < 0) d1 = d2;
    1115       27352 :     if (d1 != gen_1) {
    1116       23691 :       B1 = Q_muli_to_int(B1, d1);
    1117       23691 :       B2 = Q_muli_to_int(B2, d1);
    1118       23691 :       D = mulii(d1, D);
    1119             :     }
    1120       27352 :     fred = centermod_i(S->f, D, shifti(D,-1));
    1121       27352 :     M = cgetg(n+1, t_MAT);
    1122      166740 :     for (i=1; i<=n1; i++)
    1123      139388 :       gel(M,i) = RgX_to_RgC(FpX_rem(FpX_mul(gel(B1,i),e,D), fred, D), n);
    1124       27352 :     e = Z_ZX_sub(de, e); B2 -= n1;
    1125      111235 :     for (   ; i<=n; i++)
    1126       83883 :       gel(M,i) = RgX_to_RgC(FpX_rem(FpX_mul(gel(B2,i),e,D), fred, D), n);
    1127       27352 :     return ZpM_hnfmodid(M, p, D);
    1128             :   }
    1129             : }
    1130             : 
    1131             : /* minimum extension valuation: L/E */
    1132             : static void
    1133       55261 : vstar(GEN p,GEN h, long *L, long *E)
    1134             : {
    1135       55261 :   long first, j, k, v, w, m = degpol(h);
    1136             : 
    1137       55261 :   first = 1; k = 1; v = 0;
    1138      376006 :   for (j=1; j<=m; j++)
    1139             :   {
    1140      320745 :     GEN c = gel(h, m-j+2);
    1141      320745 :     if (signe(c))
    1142             :     {
    1143      308220 :       w = Z_pval(c,p);
    1144      308220 :       if (first || w*k < v*j) { v = w; k = j; }
    1145      308220 :       first = 0;
    1146             :     }
    1147             :   }
    1148             :   /* v/k = min_j ( v_p(h_{m-j}) / j ) */
    1149       55261 :   w = (long)ugcd(v,k);
    1150       55261 :   *L = v/w;
    1151       55261 :   *E = k/w;
    1152       55261 : }
    1153             : 
    1154             : static GEN
    1155        6159 : redelt_i(GEN a, GEN N, GEN p, GEN *pda, long *pvda)
    1156             : {
    1157             :   GEN z;
    1158        6159 :   a = Q_remove_denom(a, pda);
    1159        6159 :   *pvda = 0;
    1160        6159 :   if (*pda)
    1161             :   {
    1162        6159 :     long v = Z_pvalrem(*pda, p, &z);
    1163        6159 :     if (v) {
    1164        6159 :       *pda = powiu(p, v);
    1165        6159 :       *pvda = v;
    1166        6159 :       N  = mulii(*pda, N);
    1167             :     }
    1168             :     else
    1169           0 :       *pda = NULL;
    1170        6159 :     if (!is_pm1(z)) a = ZX_Z_mul(a, Fp_inv(z, N));
    1171             :   }
    1172        6159 :   return centermod(a, N);
    1173             : }
    1174             : /* reduce the element a modulo N [ a power of p ], taking first care of the
    1175             :  * denominators */
    1176             : static GEN
    1177        1744 : redelt(GEN a, GEN N, GEN p)
    1178             : {
    1179             :   GEN da;
    1180             :   long vda;
    1181        1744 :   a = redelt_i(a, N, p, &da, &vda);
    1182        1744 :   if (da) a = RgX_Rg_div(a, da);
    1183        1744 :   return a;
    1184             : }
    1185             : 
    1186             : /* compute the c first Newton sums modulo pp of the
    1187             :    characteristic polynomial of a/d mod chi, d > 0 power of p (NULL = gen_1),
    1188             :    a, chi in Zp[X], vda = v_p(da)
    1189             :    ns = Newton sums of chi */
    1190             : static GEN
    1191       74561 : newtonsums(GEN p, GEN a, GEN da, long vda, GEN chi, long c, GEN pp, GEN ns)
    1192             : {
    1193             :   GEN va, pa, dpa, s;
    1194       74561 :   long j, k, vdpa, lns = lg(ns);
    1195             :   pari_sp av;
    1196             : 
    1197       74561 :   a = centermod(a, pp); av = avma;
    1198       74561 :   dpa = pa = NULL; /* -Wall */
    1199       74561 :   vdpa = 0;
    1200       74561 :   va = zerovec(c);
    1201      492214 :   for (j = 1; j <= c; j++)
    1202             :   { /* pa/dpa = (a/d)^(j-1) mod (chi, pp), dpa = p^vdpa */
    1203             :     long l;
    1204      418787 :     pa = j == 1? a: FpXQ_mul(pa, a, chi, pp);
    1205      418787 :     l = lg(pa); if (l == 2) break;
    1206      418787 :     if (lns < l) l = lns;
    1207             : 
    1208      418787 :     if (da) {
    1209      409633 :       dpa = j == 1? da: mulii(dpa, da);
    1210      409633 :       vdpa += vda;
    1211      409633 :       update_den(p, &pa, &dpa, &vdpa, &pp);
    1212             :     }
    1213      418787 :     s = mulii(gel(pa,2), gel(ns,2)); /* k = 2 */
    1214      418787 :     for (k = 3; k < l; k++) s = addii(s, mulii(gel(pa,k), gel(ns,k)));
    1215      418787 :     if (da) {
    1216             :       GEN r;
    1217      409633 :       s = dvmdii(s, dpa, &r);
    1218      409633 :       if (r != gen_0) return NULL;
    1219             :     }
    1220      417653 :     gel(va,j) = centermodii(s, pp, shifti(pp,-1));
    1221             : 
    1222      417653 :     if (gc_needed(av, 1))
    1223             :     {
    1224           7 :       if(DEBUGMEM>1) pari_warn(warnmem, "newtonsums");
    1225           7 :       gerepileall(av, dpa?4:3, &pa, &va, &pp, &dpa);
    1226             :     }
    1227             :   }
    1228       73427 :   for (; j <= c; j++) gel(va,j) = gen_0;
    1229       73427 :   return va;
    1230             : }
    1231             : 
    1232             : /* compute the characteristic polynomial of a/da mod chi (a in Z[X]), given
    1233             :  * by its Newton sums to a precision of pp using Newton sums */
    1234             : static GEN
    1235       73427 : newtoncharpoly(GEN pp, GEN p, GEN NS)
    1236             : {
    1237       73427 :   long n = lg(NS)-1, j, k;
    1238       73427 :   GEN c = cgetg(n + 2, t_VEC), pp2 = shifti(pp,-1);
    1239             : 
    1240       73427 :   gel(c,1) = (n & 1 ? gen_m1: gen_1);
    1241      488559 :   for (k = 2; k <= n+1; k++)
    1242             :   {
    1243      415153 :     pari_sp av2 = avma;
    1244      415153 :     GEN s = gen_0;
    1245             :     ulong z;
    1246      415153 :     long v = u_pvalrem(k - 1, p, &z);
    1247     3187262 :     for (j = 1; j < k; j++)
    1248             :     {
    1249     2772109 :       GEN t = mulii(gel(NS,j), gel(c,k-j));
    1250     2772109 :       if (!odd(j)) t = negi(t);
    1251     2772109 :       s = addii(s, t);
    1252             :     }
    1253      415153 :     if (v) {
    1254      145203 :       s = gdiv(s, powiu(p, v));
    1255      145203 :       if (typ(s) != t_INT) return NULL;
    1256             :     }
    1257      415132 :     s = mulii(s, Fp_inv(utoipos(z), pp));
    1258      415132 :     gel(c,k) = gerepileuptoint(av2, Fp_center_i(s, pp, pp2));
    1259             :   }
    1260       73406 :   for (k = odd(n)? 1: 2; k <= n+1; k += 2) gel(c,k) = negi(gel(c,k));
    1261       73406 :   return gtopoly(c, 0);
    1262             : }
    1263             : 
    1264             : static void
    1265       74561 : manage_cache(decomp_t *S, GEN f, GEN pp)
    1266             : {
    1267       74561 :   GEN t = S->precns;
    1268             : 
    1269       74561 :   if (!t) t = mulii(S->pmf, powiu(S->p, S->df));
    1270       74561 :   if (cmpii(t, pp) < 0) t = pp;
    1271             : 
    1272       74561 :   if (!S->precns || !RgX_equal(f, S->nsf) || cmpii(S->precns, t) < 0)
    1273             :   {
    1274       44822 :     if (DEBUGLEVEL>4)
    1275           0 :       err_printf("  Precision for cached Newton sums for %Ps: %Ps -> %Ps\n",
    1276           0 :                  f, S->precns? S->precns: gen_0, t);
    1277       44822 :     S->nsf = f;
    1278       44822 :     S->ns = FpX_Newton(f, degpol(f), t);
    1279       44822 :     S->precns = t;
    1280             :   }
    1281       74561 : }
    1282             : 
    1283             : /* return NULL if a mod f is not an integer
    1284             :  * The denominator of any integer in Zp[X]/(f) divides pdr */
    1285             : static GEN
    1286       74561 : mycaract(decomp_t *S, GEN f, GEN a, GEN pp, GEN pdr)
    1287             : {
    1288             :   pari_sp av;
    1289             :   GEN d, chi, prec1, prec2, prec3, ns;
    1290       74561 :   long vd, n = degpol(f);
    1291             : 
    1292       74561 :   if (gequal0(a)) return pol_0(varn(f));
    1293             : 
    1294       74561 :   a = QpX_remove_denom(a, S->p, &d, &vd);
    1295       74561 :   prec1 = pp;
    1296       74561 :   if (lgefint(S->p) == 3)
    1297       74530 :     prec1 = mulii(prec1, powiu(S->p, factorial_lval(n, itou(S->p))));
    1298       74561 :   if (d)
    1299             :   {
    1300       71787 :     GEN p1 = powiu(d, n);
    1301       71787 :     prec2 = mulii(prec1, p1);
    1302       71787 :     prec3 = mulii(prec1, gmin_shallow(mulii(p1, d), pdr));
    1303             :   }
    1304             :   else
    1305        2774 :     prec2 = prec3 = prec1;
    1306       74561 :   manage_cache(S, f, prec3);
    1307             : 
    1308       74561 :   av = avma;
    1309       74561 :   ns = newtonsums(S->p, a, d, vd, f, n, prec2, S->ns);
    1310       74561 :   if (!ns) return NULL;
    1311       73427 :   chi = newtoncharpoly(prec1, S->p, ns);
    1312       73427 :   if (!chi) return NULL;
    1313       73406 :   setvarn(chi, varn(f));
    1314       73406 :   return gerepileupto(av, centermod(chi, pp));
    1315             : }
    1316             : 
    1317             : static GEN
    1318       64220 : get_nu(GEN chi, GEN p, long *ptl)
    1319             : { /* split off powers of x first for efficiency */
    1320       64220 :   long v = ZX_valrem(FpX_red(chi,p), &chi), n;
    1321             :   GEN P;
    1322       64220 :   if (!degpol(chi)) { *ptl = 1; return pol_x(varn(chi)); }
    1323       59099 :   P = gel(FpX_factor(chi,p), 1); n = lg(P)-1;
    1324       59099 :   *ptl = v? n+1: n; return gel(P,n);
    1325             : }
    1326             : 
    1327             : /* Factor characteristic polynomial chi of phi mod p. If it splits, update
    1328             :  * S->{phi, chi, nu} and return 1. In any case, set *nu to an irreducible
    1329             :  * factor mod p of chi */
    1330             : static int
    1331       57832 : split_char(decomp_t *S, GEN chi, GEN phi, GEN phi0, GEN *nu)
    1332             : {
    1333             :   long l;
    1334       57832 :   *nu  = get_nu(chi, S->p, &l);
    1335       57832 :   if (l == 1) return 0; /* single irreducible factor: doesn't split */
    1336             :   /* phi o phi0 mod (p, f) */
    1337       11130 :   S->phi = compmod(S->p, phi, phi0, S->f, S->p, 0);
    1338       11130 :   S->chi = chi;
    1339       11130 :   S->nu = *nu; return 1;
    1340             : }
    1341             : 
    1342             : /* Return the prime element in Zp[phi], a t_INT (iff *Ep = 1) or QX;
    1343             :  * nup, chip are ZX. phi = NULL codes X
    1344             :  * If *Ep < oE or Ep divides Ediv (!=0) return NULL (uninteresting) */
    1345             : static GEN
    1346       53273 : getprime(decomp_t *S, GEN phi, GEN chip, GEN nup, long *Lp, long *Ep,
    1347             :          long oE, long Ediv)
    1348             : {
    1349             :   GEN z, chin, q, qp;
    1350             :   long r, s;
    1351             : 
    1352       53273 :   if (phi && dvdii(constant_coeff(chip), S->psc))
    1353             :   {
    1354         196 :     chip = mycaract(S, S->chi, phi, S->pmf, S->prc);
    1355         196 :     if (dvdii(constant_coeff(chip), S->pmf))
    1356          14 :       chip = ZXQ_charpoly(phi, S->chi, varn(chip));
    1357             :   }
    1358       53273 :   if (degpol(nup) == 1)
    1359             :   {
    1360       44773 :     GEN c = gel(nup,2); /* nup = X + c */
    1361       44773 :     chin = signe(c)? RgX_translate(chip, negi(c)): chip;
    1362             :   }
    1363             :   else
    1364        8500 :     chin = ZXQ_charpoly(nup, chip, varn(chip));
    1365             : 
    1366       53273 :   vstar(S->p, chin, Lp, Ep);
    1367       53273 :   if (*Ep < oE || (Ediv && Ediv % *Ep == 0)) return NULL;
    1368             : 
    1369       30440 :   if (*Ep == 1) return S->p;
    1370       16905 :   (void)cbezout(*Lp, -*Ep, &r, &s); /* = 1 */
    1371       16905 :   if (r <= 0)
    1372             :   {
    1373        2317 :     long t = 1 + ((-r) / *Ep);
    1374        2317 :     r += t * *Ep;
    1375        2317 :     s += t * *Lp;
    1376             :   }
    1377             :   /* r > 0 minimal such that r L/E - s = 1/E
    1378             :    * pi = nu^r / p^s is an element of valuation 1/E,
    1379             :    * so is pi + O(p) since 1/E < 1. May compute nu^r mod p^(s+1) */
    1380       16905 :   q = powiu(S->p, s); qp = mulii(q, S->p);
    1381       16905 :   nup = FpXQ_powu(nup, r, S->chi, qp);
    1382       16905 :   if (!phi) return RgX_Rg_div(nup, q); /* phi = X : no composition */
    1383        1744 :   z = compmod(S->p, nup, phi, S->chi, qp, -s);
    1384        1744 :   return signe(z)? z: NULL;
    1385             : }
    1386             : 
    1387             : static int
    1388       17079 : update_phi(decomp_t *S)
    1389             : {
    1390       17079 :   GEN PHI = NULL, prc, psc, X = pol_x(varn(S->f));
    1391             :   long k;
    1392       17207 :   for (k = 1;; k++)
    1393             :   {
    1394       17335 :     prc = ZpX_reduced_resultant_fast(S->chi, ZX_deriv(S->chi), S->p, S->vpsc);
    1395       17207 :     if (!equalii(prc, S->psc)) break;
    1396             : 
    1397             :     /* increase precision */
    1398         128 :     S->vpsc = maxss(S->vpsf, S->vpsc + 1);
    1399         128 :     S->psc = (S->vpsc == S->vpsf)? S->psf: mulii(S->psc, S->p);
    1400             : 
    1401         128 :     PHI = S->phi;
    1402         128 :     if (S->phi0) PHI = compmod(S->p, PHI, S->phi0, S->f, S->psc, 0);
    1403         128 :     PHI = gadd(PHI, ZX_Z_mul(X, mului(k, S->p)));
    1404         128 :     S->chi = mycaract(S, S->f, PHI, S->psc, S->pdf);
    1405             :   }
    1406       17079 :   psc = mulii(sqri(prc), S->p);
    1407             : 
    1408       17079 :   if (!PHI) /* ok above for k = 1 */
    1409             :   {
    1410       16951 :     PHI = S->phi;
    1411       16951 :     if (S->phi0)
    1412             :     {
    1413       10563 :       PHI = compmod(S->p, PHI, S->phi0, S->f, psc, 0);
    1414       10563 :       S->chi = mycaract(S, S->f, PHI, psc, S->pdf);
    1415             :     }
    1416             :   }
    1417       17079 :   S->phi = PHI;
    1418       17079 :   S->chi = FpX_red(S->chi, psc);
    1419             : 
    1420             :   /* may happen if p is unramified */
    1421       17079 :   if (is_pm1(prc)) return 0;
    1422       12604 :   S->psc = psc;
    1423       12604 :   S->vpsc = 2*Z_pval(prc, S->p) + 1;
    1424       12604 :   S->prc = mulii(prc, S->p); return 1;
    1425             : }
    1426             : 
    1427             : /* return 1 if at least 2 factors mod p ==> chi splits
    1428             :  * Replace S->phi such that F increases (to D) */
    1429             : static int
    1430        8954 : testb2(decomp_t *S, long D, GEN theta)
    1431             : {
    1432        8954 :   long v = varn(S->chi), dlim = degpol(S->chi)-1;
    1433        8954 :   GEN T0 = S->phi, chi, phi, nu;
    1434        8954 :   if (DEBUGLEVEL>4) err_printf("  Increasing Fa\n");
    1435             :   for (;;)
    1436             :   {
    1437        8954 :     phi = gadd(theta, random_FpX(dlim, v, S->p));
    1438        8954 :     chi = mycaract(S, S->chi, phi, S->psf, S->prc);
    1439             :     /* phi non-primary ? */
    1440        8954 :     if (split_char(S, chi, phi, T0, &nu)) return 1;
    1441        8947 :     if (degpol(nu) == D) break;
    1442             :   }
    1443             :   /* F_phi=lcm(F_alpha, F_theta)=D and E_phi=E_alpha */
    1444        8947 :   S->phi0 = T0;
    1445        8947 :   S->chi = chi;
    1446        8947 :   S->phi = phi;
    1447        8947 :   S->nu = nu; return 0;
    1448             : }
    1449             : 
    1450             : /* return 1 if at least 2 factors mod p ==> chi can be split.
    1451             :  * compute a new S->phi such that E = lcm(Ea, Et);
    1452             :  * A a ZX, T a t_INT (iff Et = 1, probably impossible ?) or QX */
    1453             : static int
    1454        1744 : testc2(decomp_t *S, GEN A, long Ea, GEN T, long Et)
    1455             : {
    1456        1744 :   GEN c, chi, phi, nu, T0 = S->phi;
    1457             : 
    1458        1744 :   if (DEBUGLEVEL>4) err_printf("  Increasing Ea\n");
    1459        1744 :   if (Et == 1) /* same as other branch, split for efficiency */
    1460           0 :     c = A; /* Et = 1 => s = 1, r = 0, t = 0 */
    1461             :   else
    1462             :   {
    1463             :     long r, s, t;
    1464        1744 :     (void)cbezout(Ea, Et, &r, &s); t = 0;
    1465        1744 :     while (r < 0) { r = r + Et; t++; }
    1466        1744 :     while (s < 0) { s = s + Ea; t++; }
    1467             : 
    1468             :     /* A^s T^r / p^t */
    1469        1744 :     c = RgXQ_mul(RgXQ_powu(A, s, S->chi), RgXQ_powu(T, r, S->chi), S->chi);
    1470        1744 :     c = RgX_Rg_div(c, powiu(S->p, t));
    1471        1744 :     c = redelt(c, S->psc, S->p);
    1472             :   }
    1473        1744 :   phi = RgX_add(c,  pol_x(varn(S->chi)));
    1474        1744 :   chi = mycaract(S, S->chi, phi, S->psf, S->prc);
    1475        1744 :   if (split_char(S, chi, phi, T0, &nu)) return 1;
    1476             :   /* E_phi = lcm(E_alpha,E_theta) */
    1477        1744 :   S->phi0 = T0;
    1478        1744 :   S->chi = chi;
    1479        1744 :   S->phi = phi;
    1480        1744 :   S->nu = nu; return 0;
    1481             : }
    1482             : 
    1483             : /* Return h^(-degpol(P)) P(x * h) if result is integral, NULL otherwise */
    1484             : static GEN
    1485        1701 : ZX_rescale_inv(GEN P, GEN h)
    1486             : {
    1487        1701 :   long i, l = lg(P);
    1488        1701 :   GEN Q = cgetg(l,t_POL), hi = h;
    1489        1701 :   gel(Q,l-1) = gel(P,l-1);
    1490        9527 :   for (i=l-2; i>=2; i--)
    1491             :   {
    1492             :     GEN r;
    1493        9527 :     gel(Q,i) = dvmdii(gel(P,i), hi, &r);
    1494        9527 :     if (signe(r)) return NULL;
    1495        9527 :     if (i == 2) break;
    1496        7826 :     hi = mulii(hi,h);
    1497             :   }
    1498        1701 :   Q[1] = P[1]; return Q;
    1499             : }
    1500             : 
    1501             : /* x p^-eq nu^-er mod p */
    1502             : static GEN
    1503       43651 : get_gamma(decomp_t *S, GEN x, long eq, long er)
    1504             : {
    1505       43651 :   GEN q, g = x, Dg = powiu(S->p, eq);
    1506       43651 :   long vDg = eq;
    1507       43651 :   if (er)
    1508             :   {
    1509        8356 :     if (!S->invnu)
    1510             :     {
    1511        4415 :       while (gdvd(S->chi, S->nu)) S->nu = RgX_Rg_add(S->nu, S->p);
    1512        4415 :       S->invnu = QXQ_inv(S->nu, S->chi);
    1513        4415 :       S->invnu = redelt_i(S->invnu, S->psc, S->p, &S->Dinvnu, &S->vDinvnu);
    1514             :     }
    1515        8356 :     if (S->Dinvnu) {
    1516        8356 :       Dg = mulii(Dg, powiu(S->Dinvnu, er));
    1517        8356 :       vDg += er * S->vDinvnu;
    1518             :     }
    1519        8356 :     q = mulii(S->p, Dg);
    1520        8356 :     g = ZX_mul(g, FpXQ_powu(S->invnu, er, S->chi, q));
    1521        8356 :     g = FpX_rem(g, S->chi, q);
    1522        8356 :     update_den(S->p, &g, &Dg, &vDg, NULL);
    1523        8356 :     g = centermod(g, mulii(S->p, Dg));
    1524             :   }
    1525       43651 :   if (!is_pm1(Dg)) g = RgX_Rg_div(g, Dg);
    1526       43651 :   return g;
    1527             : }
    1528             : static GEN
    1529       44197 : get_g(decomp_t *S, long Ea, long L, long E, GEN beta, GEN *pchig,
    1530             :       long *peq, long *per)
    1531             : {
    1532             :   long eq, er;
    1533       44197 :   GEN g, chig, chib = NULL;
    1534             :   for(;;) /* at most twice */
    1535             :   {
    1536       46507 :     if (L < 0)
    1537             :     {
    1538        1988 :       chib = ZXQ_charpoly(beta, S->chi, varn(S->chi));
    1539        1988 :       vstar(S->p, chib, &L, &E);
    1540             :     }
    1541       45352 :     eq = L / E; er = L*Ea / E - eq*Ea;
    1542             :     /* floor(L Ea/E) = eq Ea + er */
    1543       45352 :     if (er || !chib)
    1544             :     { /* g might not be an integer ==> chig = NULL */
    1545       43651 :       g = get_gamma(S, beta, eq, er);
    1546       43651 :       chig = mycaract(S, S->chi, g, S->psc, S->prc);
    1547             :     }
    1548             :     else
    1549             :     { /* g = beta/p^eq, special case of the above */
    1550        1701 :       GEN h = powiu(S->p, eq);
    1551        1701 :       g = RgX_Rg_div(beta, h);
    1552        1701 :       chig = ZX_rescale_inv(chib, h); /* chib(x h) / h^N */
    1553        1701 :       if (chig) chig = FpX_red(chig, S->pmf);
    1554             :     }
    1555             :     /* either success or second consecutive failure */
    1556       45352 :     if (chig || chib) break;
    1557             :     /* if g fails the v*-test, v(beta) was wrong. Retry once */
    1558        1155 :     L = -1;
    1559             :   }
    1560       44197 :   *pchig = chig; *peq = eq; *per = er; return g;
    1561             : }
    1562             : 
    1563             : /* return 1 if at least 2 factors mod p ==> chi can be split */
    1564             : static int
    1565       21821 : loop(decomp_t *S, long Ea)
    1566             : {
    1567       21821 :   pari_sp av = avma;
    1568       21821 :   GEN beta = FpXQ_powu(S->nu, Ea, S->chi, S->p);
    1569       21821 :   long N = degpol(S->f), v = varn(S->f);
    1570       21821 :   S->invnu = NULL;
    1571             :   for (;;)
    1572       22376 :   { /* beta tends to a factor of chi */
    1573             :     long L, i, Fg, eq, er;
    1574       44197 :     GEN chig = NULL, d, g, nug;
    1575             : 
    1576       44197 :     if (DEBUGLEVEL>4) err_printf("  beta = %Ps\n", beta);
    1577       44197 :     L = ZpX_resultant_val(S->chi, beta, S->p, S->mf+1);
    1578       44197 :     if (L > S->mf) L = -1; /* from scratch */
    1579       44197 :     g = get_g(S, Ea, L, N, beta, &chig, &eq, &er);
    1580       44197 :     if (DEBUGLEVEL>4) err_printf("  (eq,er) = (%ld,%ld)\n", eq,er);
    1581             :     /* g = beta p^-eq  nu^-er (a unit), chig = charpoly(g) */
    1582       66018 :     if (split_char(S, chig, g,S->phi, &nug)) return 1;
    1583             : 
    1584       33550 :     Fg = degpol(nug);
    1585       33550 :     if (Fg == 1)
    1586             :     { /* frequent special case nug = x - d */
    1587             :       long Le, Ee;
    1588             :       GEN chie, nue, e, pie;
    1589       21659 :       d = negi(gel(nug,2));
    1590       21659 :       chie = RgX_translate(chig, d);
    1591       21659 :       nue = pol_x(v);
    1592       21659 :       e = RgX_Rg_sub(g, d);
    1593       21659 :       pie = getprime(S, e, chie, nue, &Le, &Ee,  0,Ea);
    1594       21659 :       if (pie) return testc2(S, S->nu, Ea, pie, Ee);
    1595             :     }
    1596             :     else
    1597             :     {
    1598       11891 :       long Fa = degpol(S->nu), vdeng;
    1599             :       GEN deng, numg, nume;
    1600       21503 :       if (Fa % Fg) return testb2(S, ulcm(Fa,Fg), g);
    1601             :       /* nu & nug irreducible mod p, deg nug | deg nu. To improve beta, look
    1602             :        * for a root d of nug in Fp[phi] such that v_p(g - d) > 0 */
    1603        2937 :       if (ZX_equal(nug, S->nu))
    1604        1838 :         d = pol_x(v);
    1605             :       else
    1606             :       {
    1607        1099 :         if (!p_is_prime(S)) pari_err_PRIME("FpX_ffisom",S->p);
    1608        1099 :         d = FpX_ffisom(nug, S->nu, S->p);
    1609             :       }
    1610             :       /* write g = numg / deng, e = nume / deng */
    1611        2937 :       numg = QpX_remove_denom(g, S->p, &deng, &vdeng);
    1612        5150 :       for (i = 1; i <= Fg; i++)
    1613             :       {
    1614             :         GEN chie, nue, e;
    1615        5150 :         if (i != 1) d = FpXQ_pow(d, S->p, S->nu, S->p); /* next root */
    1616        5150 :         nume = ZX_sub(numg, ZX_Z_mul(d, deng));
    1617             :         /* test e = nume / deng */
    1618        5150 :         if (ZpX_resultant_val(S->chi, nume, S->p, vdeng*N+1) <= vdeng*N)
    1619        2213 :           continue;
    1620        2937 :         e = RgX_Rg_div(nume, deng);
    1621        2937 :         chie = mycaract(S, S->chi, e, S->psc, S->prc);
    1622        3595 :         if (split_char(S, chie, e,S->phi, &nue)) return 1;
    1623        2461 :         if (RgX_is_monomial(nue))
    1624             :         { /* v_p(e) = v_p(g - d) > 0 */
    1625             :           long Le, Ee;
    1626             :           GEN pie;
    1627        2461 :           pie = getprime(S, e, chie, nue, &Le, &Ee,  0,Ea);
    1628        2461 :           if (pie) return testc2(S, S->nu, Ea, pie, Ee);
    1629        2279 :           break;
    1630             :         }
    1631             :       }
    1632        2279 :       if (i > Fg)
    1633             :       {
    1634           0 :         if (!p_is_prime(S)) pari_err_PRIME("nilord",S->p);
    1635           0 :         pari_err_BUG("nilord (no root)");
    1636             :       }
    1637             :     }
    1638       22376 :     if (eq) d = gmul(d, powiu(S->p,  eq));
    1639       22376 :     if (er) d = gmul(d, gpowgs(S->nu, er));
    1640       22376 :     beta = gsub(beta, d);
    1641             : 
    1642       22376 :     if (gc_needed(av,1))
    1643             :     {
    1644           0 :       if (DEBUGMEM > 1) pari_warn(warnmem, "nilord");
    1645           0 :       gerepileall(av, S->invnu? 6: 4, &beta, &(S->precns), &(S->ns), &(S->nsf), &(S->invnu), &(S->Dinvnu));
    1646             :     }
    1647             :   }
    1648             : }
    1649             : 
    1650             : /* E and F cannot decrease; return 1 if O = Zp[phi], 2 if we can get a
    1651             :  * decomposition and 0 otherwise */
    1652             : static long
    1653       28696 : progress(decomp_t *S, GEN *ppa, long *pE)
    1654             : {
    1655       28696 :   long E = *pE, F;
    1656       28696 :   GEN pa = *ppa;
    1657       28696 :   S->phi0 = NULL; /* no delayed composition */
    1658             :   for(;;)
    1659         457 :   {
    1660             :     long l, La, Ea; /* N.B If E = 0, getprime cannot return NULL */
    1661       29153 :     GEN pia  = getprime(S, NULL, S->chi, S->nu, &La, &Ea, E,0);
    1662       29153 :     if (pia) { /* success, we break out in THIS loop */
    1663       28696 :       pa = (typ(pia) == t_POL)? RgX_RgXQ_eval(pia, S->phi, S->f): pia;
    1664       28696 :       E = Ea;
    1665       57392 :       if (La == 1) break; /* no need to change phi so that nu = pia */
    1666             :     }
    1667             :     /* phi += prime elt */
    1668       15613 :     S->phi = typ(pa) == t_INT? RgX_Rg_add_shallow(S->phi, pa)
    1669        9225 :                              : RgX_add(S->phi, pa);
    1670             :     /* recompute char. poly. chi from scratch */
    1671        6388 :     S->chi = mycaract(S, S->f, S->phi, S->psf, S->pdf);
    1672        6388 :     S->nu = get_nu(S->chi, S->p, &l);
    1673        6388 :     if (l > 1) return 2;
    1674        6388 :     if (!update_phi(S)) return 1; /* unramified */
    1675        6388 :     if (pia) break;
    1676             :   }
    1677       28696 :   *pE = E; *ppa = pa; F = degpol(S->nu);
    1678       28696 :   if (DEBUGLEVEL>4) err_printf("  (E, F) = (%ld,%ld)\n", E, F);
    1679       28696 :   if (E * F == degpol(S->f)) return 1;
    1680       21821 :   if (loop(S, E)) return 2;
    1681       10691 :   if (!update_phi(S)) return 1;
    1682        6216 :   return 0;
    1683             : }
    1684             : 
    1685             : /* flag != 0 iff we're looking for the p-adic factorization,
    1686             :    in which case it is the p-adic precision we want */
    1687             : static GEN
    1688       38947 : maxord_i(decomp_t *S, GEN p, GEN f, long mf, GEN w, long flag)
    1689             : {
    1690       38947 :   long oE, n = lg(w)-1; /* factor of largest degree */
    1691       38947 :   GEN opa, D = ZpX_reduced_resultant_fast(f, ZX_deriv(f), p, mf);
    1692       38947 :   S->pisprime = -1;
    1693       38947 :   S->p = p;
    1694       38947 :   S->mf = mf;
    1695       38947 :   S->nu = gel(w,n);
    1696       38947 :   S->df = Z_pval(D, p);
    1697       38947 :   S->pdf = powiu(p, S->df);
    1698       38947 :   S->phi = pol_x(varn(f));
    1699       38947 :   S->chi = S->f = f;
    1700       38947 :   if (n > 1) return Decomp(S, flag); /* FIXME: use bezout_lift_fact */
    1701             : 
    1702       22480 :   if (DEBUGLEVEL>4)
    1703           0 :     err_printf("  entering Nilord: %Ps^%ld\n  f = %Ps, nu = %Ps\n",
    1704             :                p, S->df, S->f, S->nu);
    1705       22480 :   else if (DEBUGLEVEL>2) err_printf("  entering Nilord\n");
    1706       22480 :   S->psf = S->psc = mulii(sqri(D), p);
    1707       22480 :   S->vpsf = S->vpsc = 2*S->df + 1;
    1708       22480 :   S->prc = mulii(D, p);
    1709       22480 :   S->chi = FpX_red(S->f, S->psc);
    1710       22480 :   S->pmf = powiu(p, S->mf+1);
    1711       22480 :   S->precns = NULL;
    1712       22480 :   for(opa = NULL, oE = 0;;)
    1713        6216 :   {
    1714       28696 :     long n = progress(S, &opa, &oE);
    1715       28696 :     if (n == 1) return flag? NULL: dbasis(p, S->f, S->mf, S->phi, S->chi);
    1716       17346 :     if (n == 2) return Decomp(S, flag);
    1717             :   }
    1718             : }
    1719             : 
    1720             : static int
    1721         763 : expo_is_squarefree(GEN e)
    1722             : {
    1723         763 :   long i, l = lg(e);
    1724        1134 :   for (i=1; i<l; i++)
    1725         903 :     if (e[i] != 1) return 0;
    1726         231 :   return 1;
    1727             : }
    1728             : /* pure round 4 */
    1729             : static GEN
    1730         714 : ZpX_round4(GEN f, GEN p, GEN w, long prec)
    1731             : {
    1732             :   decomp_t S;
    1733         714 :   GEN L = maxord_i(&S, p, f, ZpX_disc_val(f,p), w, prec);
    1734         714 :   return L? L: mkvec(f);
    1735             : }
    1736             : /* f a squarefree ZX with leading_coeff 1, degree > 0. Return list of
    1737             :  * irreducible factors in Zp[X] (computed mod p^prec) */
    1738             : static GEN
    1739         959 : ZpX_monic_factor_squarefree(GEN f, GEN p, long prec)
    1740             : {
    1741         959 :   pari_sp av = avma;
    1742             :   GEN L, fa, w, e;
    1743             :   long i, l;
    1744         959 :   if (degpol(f) == 1) return mkvec(f);
    1745         763 :   fa = FpX_factor(f,p); w = gel(fa,1); e = gel(fa,2);
    1746             :   /* no repeated factors: Hensel lift */
    1747         763 :   if (expo_is_squarefree(e)) return ZpX_liftfact(f, w, powiu(p,prec), p, prec);
    1748         532 :   l = lg(w);
    1749         532 :   if (l == 2)
    1750             :   {
    1751         336 :     L = ZpX_round4(f,p,w,prec);
    1752         336 :     if (lg(L) == 2) { set_avma(av); return mkvec(f); }
    1753             :   }
    1754             :   else
    1755             :   { /* >= 2 factors mod p: partial Hensel lift */
    1756         196 :     GEN D = ZpX_reduced_resultant_fast(f, ZX_deriv(f), p, ZpX_disc_val(f,p));
    1757         196 :     long r = maxss(2*Z_pval(D,p)+1, prec);
    1758         196 :     GEN W = cgetg(l, t_VEC);
    1759         637 :     for (i = 1; i < l; i++)
    1760         441 :       gel(W,i) = e[i] == 1? gel(w,i): FpX_powu(gel(w,i), e[i], p);
    1761         196 :     L = ZpX_liftfact(f, W, powiu(p,r), p, r);
    1762         637 :     for (i = 1; i < l; i++)
    1763         945 :       gel(L,i) = e[i] == 1? mkvec(gel(L,i))
    1764         504 :                           : ZpX_round4(gel(L,i), p, mkvec(gel(w,i)), prec);
    1765         196 :     L = shallowconcat1(L);
    1766             :   }
    1767         315 :   return gerepilecopy(av, L);
    1768             : }
    1769             : 
    1770             : /* assume f a ZX with leading_coeff 1, degree > 0 */
    1771             : GEN
    1772         462 : ZpX_monic_factor(GEN f, GEN p, long prec)
    1773             : {
    1774             :   GEN poly, ex, P, E;
    1775             :   long l, i;
    1776             : 
    1777         462 :   if (degpol(f) == 1) return mkmat2(mkcol(f), mkcol(gen_1));
    1778         462 :   poly = ZX_squff(f,&ex); l = lg(poly);
    1779         462 :   P = cgetg(l, t_VEC);
    1780         462 :   E = cgetg(l, t_VEC);
    1781         931 :   for (i = 1; i < l; i++)
    1782             :   {
    1783         469 :     GEN L = ZpX_monic_factor_squarefree(gel(poly,i), p, prec);
    1784         469 :     gel(P,i) = L; settyp(L, t_COL);
    1785         469 :     gel(E,i) = const_col(lg(L)-1, utoipos(ex[i]));
    1786             :   }
    1787         462 :   return mkmat2(shallowconcat1(P), shallowconcat1(E));
    1788             : }
    1789             : 
    1790             : /* DT = multiple of disc(T) or NULL
    1791             :  * Return a multiple of the denominator of an algebraic integer (in Q[X]/(T))
    1792             :  * when expressed in terms of the power basis */
    1793             : GEN
    1794        2604 : indexpartial(GEN T, GEN DT)
    1795             : {
    1796        2604 :   pari_sp av = avma;
    1797             :   long i, nb;
    1798        2604 :   GEN fa, E, P, res = gen_1, dT = ZX_deriv(T);
    1799             : 
    1800        2604 :   if (!DT) DT = ZX_disc(T);
    1801        2604 :   fa = absZ_factor_limit(DT, 0);
    1802        2604 :   P = gel(fa,1);
    1803        2604 :   E = gel(fa,2); nb = lg(P)-1;
    1804       17787 :   for (i = 1; i <= nb; i++)
    1805             :   {
    1806       15183 :     long e = itou(gel(E,i)), e2 = e >> 1;
    1807       15183 :     GEN p = gel(P,i), q = p;
    1808       15183 :     if (i == nb)
    1809        2590 :       q = powiu(p, (odd(e) && !BPSW_psp(p))? e2+1: e2);
    1810       12593 :     else if (e2 >= 2)
    1811        7301 :       q = ZpX_reduced_resultant_fast(T, dT, p, e2);
    1812       15183 :     res = mulii(res, q);
    1813             :   }
    1814        2604 :   return gerepileuptoint(av,res);
    1815             : }
    1816             : 
    1817             : /*******************************************************************/
    1818             : /*                                                                 */
    1819             : /*    2-ELT REPRESENTATION FOR PRIME IDEALS (dividing index)       */
    1820             : /*                                                                 */
    1821             : /*******************************************************************/
    1822             : /* to compute norm of elt in basis form */
    1823             : typedef struct {
    1824             :   long r1;
    1825             :   GEN M;  /* via embed_norm */
    1826             : 
    1827             :   GEN D, w, T; /* via resultant if M = NULL */
    1828             : } norm_S;
    1829             : 
    1830             : static GEN
    1831       72945 : get_norm(norm_S *S, GEN a)
    1832             : {
    1833       72945 :   if (S->M)
    1834             :   {
    1835             :     long e;
    1836       72279 :     GEN N = grndtoi( embed_norm(RgM_RgC_mul(S->M, a), S->r1), &e );
    1837       72279 :     if (e > -5) pari_err_PREC( "get_norm");
    1838       72279 :     return N;
    1839             :   }
    1840         666 :   if (S->w) a = RgV_RgC_mul(S->w, a);
    1841         666 :   return ZX_resultant_all(S->T, a, S->D, 0);
    1842             : }
    1843             : static void
    1844       21042 : init_norm(norm_S *S, GEN nf, GEN p)
    1845             : {
    1846       21042 :   GEN T = nf_get_pol(nf), M = nf_get_M(nf);
    1847       21042 :   long N = degpol(T), ex = gexpo(M) + gexpo(mului(8 * N, p));
    1848             : 
    1849       21042 :   S->r1 = nf_get_r1(nf);
    1850       21042 :   if (N * ex <= prec2nbits(gprecision(M)) - 20)
    1851             :   { /* enough prec to use embed_norm */
    1852       20983 :     S->M = M;
    1853       20983 :     S->D = NULL;
    1854       20983 :     S->w = NULL;
    1855       20983 :     S->T = NULL;
    1856             :   }
    1857             :   else
    1858             :   {
    1859          59 :     GEN w = leafcopy(nf_get_zkprimpart(nf)), D = nf_get_zkden(nf), Dp = sqri(p);
    1860             :     long i;
    1861          59 :     if (!equali1(D))
    1862             :     {
    1863          59 :       GEN w1 = D;
    1864          59 :       long v = Z_pval(D, p);
    1865          59 :       D = powiu(p, v);
    1866          59 :       Dp = mulii(D, Dp);
    1867          59 :       gel(w, 1) = remii(w1, Dp);
    1868             :     }
    1869          59 :     for (i=2; i<=N; i++) gel(w,i) = FpX_red(gel(w,i), Dp);
    1870          59 :     S->M = NULL;
    1871          59 :     S->D = D;
    1872          59 :     S->w = w;
    1873          59 :     S->T = T;
    1874             :   }
    1875       21042 : }
    1876             : /* f = f(pr/p), q = p^(f+1), a in pr.
    1877             :  * Return 1 if v_pr(a) = 1, and 0 otherwise */
    1878             : static int
    1879       72945 : is_uniformizer(GEN a, GEN q, norm_S *S) { return !dvdii(get_norm(S,a), q); }
    1880             : 
    1881             : /* Return x * y, x, y are t_MAT (Fp-basis of in O_K/p), assume (x,y)=1.
    1882             :  * Either x or y may be NULL (= O_K), not both */
    1883             : static GEN
    1884      174974 : mul_intersect(GEN x, GEN y, GEN p)
    1885             : {
    1886      174974 :   if (!x) return y;
    1887      121283 :   if (!y) return x;
    1888      103386 :   return FpM_intersect(x, y, p);
    1889             : }
    1890             : /* Fp-basis of (ZK/pr): applied to the primes found in primedec_aux()
    1891             :  * true nf */
    1892             : static GEN
    1893       70256 : Fp_basis(GEN nf, GEN pr)
    1894             : {
    1895             :   long i, j, l;
    1896             :   GEN x, y;
    1897             :   /* already in basis form (from Buchman-Lenstra) ? */
    1898       70256 :   if (typ(pr) == t_MAT) return pr;
    1899             :   /* ordinary prid (from Kummer) */
    1900       17547 :   x = pr_hnf(nf, pr);
    1901       17547 :   l = lg(x);
    1902       17547 :   y = cgetg(l, t_MAT);
    1903      205784 :   for (i=j=1; i<l; i++)
    1904      188237 :     if (gequal1(gcoeff(x,i,i))) gel(y,j++) = gel(x,i);
    1905       17547 :   setlg(y, j); return y;
    1906             : }
    1907             : /* Let Ip = prod_{ P | p } P be the p-radical. The list L contains the
    1908             :  * P (mod Ip) seen as sub-Fp-vector spaces of ZK/Ip.
    1909             :  * Return the list of (Ip / P) (mod Ip).
    1910             :  * N.B: All ideal multiplications are computed as intersections of Fp-vector
    1911             :  * spaces. true nf */
    1912             : static GEN
    1913       21042 : get_LV(GEN nf, GEN L, GEN p, long N)
    1914             : {
    1915       21042 :   long i, l = lg(L)-1;
    1916             :   GEN LV, LW, A, B;
    1917             : 
    1918       21042 :   LV = cgetg(l+1, t_VEC);
    1919       21042 :   if (l == 1) { gel(LV,1) = matid(N); return LV; }
    1920       17897 :   LW = cgetg(l+1, t_VEC);
    1921       17897 :   for (i=1; i<=l; i++) gel(LW,i) = Fp_basis(nf, gel(L,i));
    1922             : 
    1923             :   /* A[i] = L[1]...L[i-1], i = 2..l */
    1924       17897 :   A = cgetg(l+1, t_VEC); gel(A,1) = NULL;
    1925       17897 :   for (i=1; i < l; i++) gel(A,i+1) = mul_intersect(gel(A,i), gel(LW,i), p);
    1926             :   /* B[i] = L[i+1]...L[l], i = 1..(l-1) */
    1927       17897 :   B = cgetg(l+1, t_VEC); gel(B,l) = NULL;
    1928       17897 :   for (i=l; i>=2; i--) gel(B,i-1) = mul_intersect(gel(B,i), gel(LW,i), p);
    1929       17897 :   for (i=1; i<=l; i++) gel(LV,i) = mul_intersect(gel(A,i), gel(B,i), p);
    1930       17897 :   return LV;
    1931             : }
    1932             : 
    1933             : static void
    1934           0 : errprime(GEN p) { pari_err_PRIME("idealprimedec",p); }
    1935             : 
    1936             : /* P = Fp-basis (over O_K/p) for pr.
    1937             :  * V = Z-basis for I_p/pr. ramif != 0 iff some pr|p is ramified.
    1938             :  * Return a p-uniformizer for pr. Assume pr not inert, i.e. m > 0 */
    1939             : static GEN
    1940       53951 : uniformizer(GEN nf, norm_S *S, GEN P, GEN V, GEN p, int ramif)
    1941             : {
    1942       53951 :   long i, l, f, m = lg(P)-1, N = nf_get_degree(nf);
    1943             :   GEN u, Mv, x, q;
    1944             : 
    1945       53951 :   f = N - m; /* we want v_p(Norm(x)) = p^f */
    1946       53951 :   q = powiu(p,f+1);
    1947             : 
    1948       53951 :   u = FpM_FpC_invimage(shallowconcat(P, V), col_ei(N,1), p);
    1949       53951 :   setlg(u, lg(P));
    1950       53951 :   u = centermod(ZM_ZC_mul(P, u), p);
    1951       53951 :   if (is_uniformizer(u, q, S)) return u;
    1952       14909 :   if (signe(gel(u,1)) <= 0) /* make sure u[1] in ]-p,p] */
    1953       11660 :     gel(u,1) = addii(gel(u,1), p); /* try u + p */
    1954             :   else
    1955        3249 :     gel(u,1) = subii(gel(u,1), p); /* try u - p */
    1956       14909 :   if (!ramif || is_uniformizer(u, q, S)) return u;
    1957             : 
    1958             :   /* P/p ramified, u in P^2, not in Q for all other Q|p */
    1959        5535 :   Mv = zk_multable(nf, Z_ZC_sub(gen_1,u));
    1960        5535 :   l = lg(P);
    1961       13410 :   for (i=1; i<l; i++)
    1962             :   {
    1963       13410 :     x = centermod(ZC_add(u, ZM_ZC_mul(Mv, gel(P,i))), p);
    1964       13410 :     if (is_uniformizer(x, q, S)) return x;
    1965             :   }
    1966           0 :   errprime(p);
    1967             :   return NULL; /* LCOV_EXCL_LINE */
    1968             : }
    1969             : 
    1970             : /*******************************************************************/
    1971             : /*                                                                 */
    1972             : /*                   BUCHMANN-LENSTRA ALGORITHM                    */
    1973             : /*                                                                 */
    1974             : /*******************************************************************/
    1975             : static GEN
    1976      988602 : mk_pr(GEN p, GEN u, long e, long f, GEN t)
    1977      988602 : { return mkvec5(p, u, utoipos(e), utoipos(f), t); }
    1978             : 
    1979             : /* nf a true nf; pr = (p,u) of ramification index e */
    1980             : GEN
    1981      919004 : idealprimedec_kummer(GEN nf,GEN u,long e,GEN p)
    1982             : {
    1983      919004 :   GEN t, T = nf_get_pol(nf);
    1984      918993 :   long f = degpol(u), N = degpol(T);
    1985             : 
    1986      918972 :   if (f == N) /* inert */
    1987             :   {
    1988      136430 :     u = scalarcol_shallow(p,N);
    1989      136430 :     t = gen_1;
    1990             :   }
    1991             :   else
    1992             :   { /* make sure v_pr(u) = 1 (automatic if e>1) */
    1993      782542 :     t = poltobasis(nf, FpX_div(T,u,p));
    1994      782480 :     t = centermod(t, p);
    1995      782517 :     u = FpX_center_i(u, p, shifti(p,-1));
    1996      782562 :     if (e == 1 && ZpX_resultant_val(T, u, p, f+1) > f)
    1997       49404 :       gel(u,2) = addii(gel(u,2), p);
    1998      782563 :     u = poltobasis(nf,u);
    1999      782489 :     t = zk_multable(nf, t); /* t never a scalar here since pr is not inert */
    2000             :   }
    2001      918953 :   return mk_pr(p,u,e,f,t);
    2002             : }
    2003             : 
    2004             : typedef struct {
    2005             :   GEN nf, p;
    2006             :   long I;
    2007             : } eltmod_muldata;
    2008             : 
    2009             : static GEN
    2010      249177 : sqr_mod(void *data, GEN x)
    2011             : {
    2012      249177 :   eltmod_muldata *D = (eltmod_muldata*)data;
    2013      249177 :   return FpC_red(nfsqri(D->nf, x), D->p);
    2014             : }
    2015             : static GEN
    2016      123583 : ei_msqr_mod(void *data, GEN x)
    2017             : {
    2018      123583 :   GEN x2 = sqr_mod(data, x);
    2019      123583 :   eltmod_muldata *D = (eltmod_muldata*)data;
    2020      123583 :   return FpC_red(zk_ei_mul(D->nf, x2, D->I), D->p);
    2021             : }
    2022             : /* nf a true nf; compute lift(nf.zk[I]^p mod p) */
    2023             : static GEN
    2024      144763 : pow_ei_mod_p(GEN nf, long I, GEN p)
    2025             : {
    2026      144763 :   pari_sp av = avma;
    2027             :   eltmod_muldata D;
    2028      144763 :   long N = nf_get_degree(nf);
    2029      144763 :   GEN y = col_ei(N,I);
    2030      144763 :   if (I == 1) return y;
    2031      123427 :   D.nf = nf;
    2032      123427 :   D.p = p;
    2033      123427 :   D.I = I;
    2034      123427 :   y = gen_pow_fold(y, p, (void*)&D, &sqr_mod, &ei_msqr_mod);
    2035      123427 :   return gerepileupto(av,y);
    2036             : }
    2037             : 
    2038             : /* nf a true nf; return a Z basis of Z_K's p-radical, phi = x--> x^p-x */
    2039             : static GEN
    2040       21042 : pradical(GEN nf, GEN p, GEN *phi)
    2041             : {
    2042       21042 :   long i, N = nf_get_degree(nf);
    2043             :   GEN q,m,frob,rad;
    2044             : 
    2045             :   /* matrix of Frob: x->x^p over Z_K/p */
    2046       21042 :   frob = cgetg(N+1,t_MAT);
    2047       21042 :   for (i=1; i<=N; i++) gel(frob,i) = pow_ei_mod_p(nf,i,p);
    2048             : 
    2049       21042 :   m = frob; q = p;
    2050       21042 :   while (abscmpiu(q,N) < 0) { q = mulii(q,p); m = FpM_mul(m, frob, p); }
    2051       21042 :   rad = FpM_ker(m, p); /* m = Frob^k, s.t p^k >= N */
    2052       21042 :   for (i=1; i<=N; i++) gcoeff(frob,i,i) = subiu(gcoeff(frob,i,i), 1);
    2053       21042 :   *phi = frob; return rad;
    2054             : }
    2055             : 
    2056             : /* return powers of a: a^0, ... , a^d,  d = dim A */
    2057             : static GEN
    2058       32405 : get_powers(GEN mul, GEN p)
    2059             : {
    2060       32405 :   long i, d = lgcols(mul);
    2061       32405 :   GEN z, pow = cgetg(d+2,t_MAT), P = pow+1;
    2062             : 
    2063       32405 :   gel(P,0) = scalarcol_shallow(gen_1, d-1);
    2064       32405 :   z = gel(mul,1);
    2065      182652 :   for (i=1; i<=d; i++)
    2066             :   {
    2067      150247 :     gel(P,i) = z; /* a^i */
    2068      150247 :     if (i!=d) z = FpM_FpC_mul(mul, z, p);
    2069             :   }
    2070       32405 :   return pow;
    2071             : }
    2072             : 
    2073             : /* minimal polynomial of a in A (dim A = d).
    2074             :  * mul = multiplication table by a in A */
    2075             : static GEN
    2076       29458 : pol_min(GEN mul, GEN p)
    2077             : {
    2078       29458 :   pari_sp av = avma;
    2079       29458 :   GEN z = FpM_deplin(get_powers(mul, p), p);
    2080       29458 :   return gerepilecopy(av, RgV_to_RgX(z,0));
    2081             : }
    2082             : 
    2083             : static GEN
    2084       73121 : get_pr(GEN nf, norm_S *S, GEN p, GEN P, GEN V, int ramif, long N, long flim)
    2085             : {
    2086             :   GEN u, t;
    2087             :   long e, f;
    2088             : 
    2089       73121 :   if (typ(P) == t_VEC)
    2090             :   { /* already done (Kummer) */
    2091       17547 :     f = pr_get_f(P);
    2092       17547 :     if (flim > 0 && f > flim) return NULL;
    2093       16742 :     if (flim == -2) return (GEN)f;
    2094       16742 :     return P;
    2095             :   }
    2096       55574 :   f = N - (lg(P)-1);
    2097       55574 :   if (flim > 0 && f > flim) return NULL;
    2098       54903 :   if (flim == -2) return (GEN)f;
    2099             :   /* P = (p,u) prime. t is an anti-uniformizer: Z_K + t/p Z_K = P^(-1),
    2100             :    * so that v_P(t) = e(P/p)-1 */
    2101       54693 :   if (f == N) {
    2102         742 :     u = scalarcol_shallow(p,N);
    2103         742 :     t = gen_1;
    2104         742 :     e = 1;
    2105             :   } else {
    2106             :     GEN mt;
    2107       53951 :     u = uniformizer(nf, S, P, V, p, ramif);
    2108       53951 :     t = FpM_deplin(zk_multable(nf,u), p);
    2109       53951 :     mt = zk_multable(nf, t);
    2110       53951 :     e = ramif? 1 + ZC_nfval(t,mk_pr(p,u,0,0,mt)): 1;
    2111       53951 :     t = mt;
    2112             :   }
    2113       54693 :   return mk_pr(p,u,e,f,t);
    2114             : }
    2115             : 
    2116             : /* true nf */
    2117             : static GEN
    2118       21042 : primedec_end(GEN nf, GEN L, GEN p, long flim)
    2119             : {
    2120       21042 :   long i, j, l = lg(L), N = nf_get_degree(nf);
    2121       21042 :   GEN LV = get_LV(nf, L,p,N);
    2122       21042 :   int ramif = dvdii(nf_get_disc(nf), p);
    2123       21042 :   norm_S S; init_norm(&S, nf, p);
    2124       93806 :   for (i = j = 1; i < l; i++)
    2125             :   {
    2126       73121 :     GEN P = get_pr(nf, &S, p, gel(L,i), gel(LV,i), ramif, N, flim);
    2127       73121 :     if (!P) continue;
    2128       71645 :     gel(L,j++) = P;
    2129       71645 :     if (flim == -1) return P;
    2130             :   }
    2131       20685 :   setlg(L, j); return L;
    2132             : }
    2133             : 
    2134             : /* prime ideal decomposition of p; if flim>0, restrict to f(P,p) <= flim
    2135             :  * if flim = -1 return only the first P
    2136             :  * if flim = -2 return only the f(P/p) in a t_VECSMALL */
    2137             : static GEN
    2138      698221 : primedec_aux(GEN nf, GEN p, long flim)
    2139             : {
    2140      698221 :   const long TYP = (flim == -2)? t_VECSMALL: t_VEC;
    2141      698221 :   GEN E, F, L, Ip, phi, f, g, h, UN, T = nf_get_pol(nf);
    2142             :   long i, k, c, iL, N;
    2143             :   int kummer;
    2144             : 
    2145      698220 :   F = FpX_factor(T, p);
    2146      698260 :   E = gel(F,2);
    2147      698260 :   F = gel(F,1);
    2148             : 
    2149      698260 :   k = lg(F); if (k == 1) errprime(p);
    2150      698260 :   if ( !dvdii(nf_get_index(nf),p) ) /* p doesn't divide index */
    2151             :   {
    2152      675803 :     L = cgetg(k, TYP);
    2153     1569549 :     for (i=1; i<k; i++)
    2154             :     {
    2155     1113220 :       GEN t = gel(F,i);
    2156     1113220 :       long f = degpol(t);
    2157     1113298 :       if (flim > 0 && f > flim) { setlg(L, i); break; }
    2158      898168 :       if (flim == -2)
    2159           1 :         L[i] = f;
    2160             :       else
    2161      898167 :         gel(L,i) = idealprimedec_kummer(nf, t, E[i],p);
    2162      898122 :       if (flim == -1) return gel(L,1);
    2163             :     }
    2164      671455 :     return L;
    2165             :   }
    2166             : 
    2167       22351 :   kummer = 0;
    2168       22351 :   g = FpXV_prod(F, p);
    2169       22351 :   h = FpX_div(T,g,p);
    2170       22351 :   f = FpX_red(ZX_Z_divexact(ZX_sub(ZX_mul(g,h), T), p), p);
    2171             : 
    2172       22351 :   N = degpol(T);
    2173       22351 :   L = cgetg(N+1,TYP);
    2174       22351 :   iL = 1;
    2175       72915 :   for (i=1; i<k; i++)
    2176       51873 :     if (E[i] == 1 || signe(FpX_rem(f,gel(F,i),p)))
    2177       17547 :     {
    2178       18856 :       GEN t = gel(F,i);
    2179       18856 :       kummer = 1;
    2180       18856 :       gel(L,iL++) = idealprimedec_kummer(nf, t, E[i],p);
    2181       18856 :       if (flim == -1) return gel(L,1);
    2182             :     }
    2183             :     else /* F[i] | (f,g,h), happens at least once by Dedekind criterion */
    2184       33017 :       E[i] = 0;
    2185             : 
    2186             :   /* phi matrix of x -> x^p - x in algebra Z_K/p */
    2187       21042 :   Ip = pradical(nf,p,&phi);
    2188             : 
    2189             :   /* split etale algebra Z_K / (p,Ip) */
    2190       21042 :   h = cgetg(N+1,t_VEC);
    2191       21042 :   if (kummer)
    2192             :   { /* split off Kummer factors */
    2193        7428 :     GEN mb, b = NULL;
    2194       36907 :     for (i=1; i<k; i++)
    2195       29479 :       if (!E[i]) b = b? FpX_mul(b, gel(F,i), p): gel(F,i);
    2196        7428 :     if (!b) errprime(p);
    2197        7428 :     b = FpC_red(poltobasis(nf,b), p);
    2198        7428 :     mb = FpM_red(zk_multable(nf,b), p);
    2199             :     /* Fp-base of ideal (Ip, b) in ZK/p */
    2200        7428 :     gel(h,1) = FpM_image(shallowconcat(mb,Ip), p);
    2201             :   }
    2202             :   else
    2203       13614 :     gel(h,1) = Ip;
    2204             : 
    2205       21042 :   UN = col_ei(N, 1);
    2206       58530 :   for (c=1; c; c--)
    2207             :   { /* Let A:= (Z_K/p) / Ip etale; split A2 := A / Im H ~ Im M2
    2208             :        H * ? + M2 * Mi2 = Id_N ==> M2 * Mi2 projector A --> A2 */
    2209       37488 :     GEN M, Mi, M2, Mi2, phi2, mat1, H = gel(h,c); /* maximal rank */
    2210       37488 :     long dim, r = lg(H)-1;
    2211             : 
    2212       37488 :     M   = FpM_suppl(shallowconcat(H,UN), p);
    2213       37488 :     Mi  = FpM_inv(M, p);
    2214       37488 :     M2  = vecslice(M, r+1,N); /* M = (H|M2) invertible */
    2215       37488 :     Mi2 = rowslice(Mi,r+1,N);
    2216             :     /* FIXME: FpM_mul(,M2) could be done with vecpermute */
    2217       37488 :     phi2 = FpM_mul(Mi2, FpM_mul(phi,M2, p), p);
    2218       37488 :     mat1 = FpM_ker(phi2, p);
    2219       37488 :     dim = lg(mat1)-1; /* A2 product of 'dim' fields */
    2220       37488 :     if (dim > 1)
    2221             :     { /* phi2 v = 0 => a = M2 v in Ker phi, a not in Fp.1 + H */
    2222       29458 :       GEN R, a, mula, mul2, v = gel(mat1,2);
    2223             :       long n;
    2224             : 
    2225       29458 :       a = FpM_FpC_mul(M2,v, p); /* not a scalar */
    2226       29458 :       mula = FpM_red(zk_multable(nf,a), p);
    2227       29458 :       mul2 = FpM_mul(Mi2, FpM_mul(mula,M2, p), p);
    2228       29458 :       R = FpX_roots(pol_min(mul2,p), p); /* totally split mod p */
    2229       29458 :       n = lg(R)-1;
    2230       93728 :       for (i=1; i<=n; i++)
    2231             :       {
    2232       64270 :         GEN I = RgM_Rg_sub_shallow(mula, gel(R,i));
    2233       64270 :         gel(h,c++) = FpM_image(shallowconcat(H, I), p);
    2234             :       }
    2235       29458 :       if (n == dim)
    2236       22103 :         for (i=1; i<=n; i++) gel(L,iL++) = gel(h,--c);
    2237             :     }
    2238             :     else /* A2 field ==> H maximal, f = N-r = dim(A2) */
    2239        8030 :       gel(L,iL++) = H;
    2240             :   }
    2241       21042 :   setlg(L, iL);
    2242       21042 :   return primedec_end(nf, L, p, flim);
    2243             : }
    2244             : 
    2245             : GEN
    2246      691964 : idealprimedec_limit_f(GEN nf, GEN p, long f)
    2247             : {
    2248      691964 :   pari_sp av = avma;
    2249             :   GEN v;
    2250      691964 :   if (typ(p) != t_INT) pari_err_TYPE("idealprimedec",p);
    2251      691965 :   if (f < 0) pari_err_DOMAIN("idealprimedec", "f", "<", gen_0, stoi(f));
    2252      691965 :   v = primedec_aux(checknf(nf), p, f);
    2253      691928 :   v = gen_sort(v, (void*)&cmp_prime_over_p, &cmp_nodata);
    2254      691986 :   return gerepileupto(av,v);
    2255             : }
    2256             : GEN
    2257        6069 : idealprimedec_galois(GEN nf, GEN p)
    2258             : {
    2259        6069 :   pari_sp av = avma;
    2260        6069 :   GEN v = primedec_aux(checknf(nf), p, -1);
    2261        6069 :   return gerepilecopy(av,v);
    2262             : }
    2263             : GEN
    2264         203 : idealprimedec_degrees(GEN nf, GEN p)
    2265             : {
    2266         203 :   pari_sp av = avma;
    2267         203 :   GEN v = primedec_aux(checknf(nf), p, -2);
    2268         203 :   vecsmall_sort(v); return gerepileuptoleaf(av, v);
    2269             : }
    2270             : GEN
    2271      201166 : idealprimedec_limit_norm(GEN nf, GEN p, GEN B)
    2272      201166 : { return idealprimedec_limit_f(nf, p, logint(B,p)); }
    2273             : GEN
    2274      165654 : idealprimedec(GEN nf, GEN p)
    2275      165654 : { return idealprimedec_limit_f(nf, p, 0); }
    2276             : GEN
    2277        1120 : nf_pV_to_prV(GEN nf, GEN P)
    2278             : {
    2279             :   long i, l;
    2280        1120 :   GEN Q = cgetg_copy(P,&l);
    2281        1120 :   if (l == 1) return Q;
    2282        1078 :   for (i = 1; i < l; i++) gel(Q,i) = idealprimedec(nf, gel(P,i));
    2283        1078 :   return shallowconcat1(Q);
    2284             : }
    2285             : 
    2286             : /* return [Fp[x]: Fp] */
    2287             : static long
    2288         448 : ffdegree(GEN x, GEN frob, GEN p)
    2289             : {
    2290         448 :   pari_sp av = avma;
    2291         448 :   long d, f = lg(frob)-1;
    2292         448 :   GEN y = x;
    2293             : 
    2294        1827 :   for (d=1; d < f; d++)
    2295             :   {
    2296        1533 :     y = FpM_FpC_mul(frob, y, p);
    2297        1533 :     if (ZV_equal(y, x)) break;
    2298             :   }
    2299         448 :   return gc_long(av,d);
    2300             : }
    2301             : 
    2302             : static GEN
    2303        8043 : lift_to_zk(GEN v, GEN c, long N)
    2304             : {
    2305        8043 :   GEN w = zerocol(N);
    2306        8043 :   long i, l = lg(c);
    2307        8043 :   for (i=1; i<l; i++) gel(w,c[i]) = gel(v,i);
    2308        8043 :   return w;
    2309             : }
    2310             : 
    2311             : /* return t = 1 mod pr, t = 0 mod p / pr^e(pr/p) */
    2312             : static GEN
    2313      375428 : anti_uniformizer(GEN nf, GEN pr)
    2314             : {
    2315      375428 :   long N = nf_get_degree(nf), e = pr_get_e(pr);
    2316             :   GEN p, b, z;
    2317             : 
    2318      375414 :   if (e * pr_get_f(pr) == N) return gen_1;
    2319      104172 :   p = pr_get_p(pr);
    2320      104168 :   b = pr_get_tau(pr); /* ZM */
    2321      104165 :   if (e != 1)
    2322             :   {
    2323        2744 :     GEN q = powiu(pr_get_p(pr), e-1);
    2324        2744 :     b = ZM_Z_divexact(ZM_powu(b,e), q);
    2325             :   }
    2326             :   /* b = tau^e / p^(e-1), v_pr(b) = 0, v_Q(b) >= e(Q/p) for other Q | p */
    2327      104165 :   z = ZM_hnfmodid(FpM_red(b,p), p); /* ideal (p) / pr^e, coprime to pr */
    2328      104184 :   z = idealaddtoone_raw(nf, pr, z);
    2329      104185 :   return Z_ZC_sub(gen_1, FpC_center(FpC_red(z,p), p, shifti(p,-1)));
    2330             : }
    2331             : 
    2332             : #define mpr_TAU 1
    2333             : #define mpr_FFP 2
    2334             : #define mpr_NFP 5
    2335             : #define SMALLMODPR 4
    2336             : #define LARGEMODPR 6
    2337             : static GEN
    2338      858396 : modpr_TAU(GEN modpr)
    2339             : {
    2340      858396 :   GEN tau = gel(modpr,mpr_TAU);
    2341      858396 :   return isintzero(tau)? NULL: tau;
    2342             : }
    2343             : 
    2344             : /* prh = HNF matrix, which is identity but for the first line. Return a
    2345             :  * projector to ZK / prh ~ Z/prh[1,1] */
    2346             : GEN
    2347      377269 : dim1proj(GEN prh)
    2348             : {
    2349      377269 :   long i, N = lg(prh)-1;
    2350      377269 :   GEN ffproj = cgetg(N+1, t_VEC);
    2351      377324 :   GEN x, q = gcoeff(prh,1,1);
    2352      377324 :   gel(ffproj,1) = gen_1;
    2353      572441 :   for (i=2; i<=N; i++)
    2354             :   {
    2355      195225 :     x = gcoeff(prh,1,i);
    2356      195225 :     if (signe(x)) x = subii(q,x);
    2357      195117 :     gel(ffproj,i) = x;
    2358             :   }
    2359      377216 :   return ffproj;
    2360             : }
    2361             : 
    2362             : /* p not necessarily prime, but coprime to denom(basis) */
    2363             : GEN
    2364         140 : QXQV_to_FpM(GEN basis, GEN T, GEN p)
    2365             : {
    2366         140 :   long i, l = lg(basis), f = degpol(T);
    2367         140 :   GEN z = cgetg(l, t_MAT);
    2368        3444 :   for (i = 1; i < l; i++)
    2369             :   {
    2370        3304 :     GEN w = gel(basis,i);
    2371        3304 :     if (typ(w) == t_INT)
    2372           0 :       w = scalarcol_shallow(w, f);
    2373             :     else
    2374             :     {
    2375             :       GEN dx;
    2376        3304 :       w = Q_remove_denom(w, &dx);
    2377        3304 :       w = FpXQ_red(w, T, p);
    2378        3304 :       if (dx)
    2379             :       {
    2380           0 :         dx = Fp_inv(dx, p);
    2381           0 :         if (!equali1(dx)) w = FpX_Fp_mul(w, dx, p);
    2382             :       }
    2383        3304 :       w = RgX_to_RgC(w, f);
    2384             :     }
    2385        3304 :     gel(z,i) = w; /* w_i mod (T,p) */
    2386             :   }
    2387         140 :   return z;
    2388             : }
    2389             : 
    2390             : /* initialize reduction mod pr; if zk = 1, will only init data required to
    2391             :  * reduce *integral* element.  Realize (O_K/pr) as Fp[X] / (T), for a
    2392             :  * *monic* T; use variable vT for varn(T) */
    2393             : static GEN
    2394      393240 : modprinit(GEN nf, GEN pr, int zk, long vT)
    2395             : {
    2396      393240 :   pari_sp av = avma;
    2397             :   GEN res, tau, mul, x, p, T, pow, ffproj, nfproj, prh, c;
    2398             :   long N, i, k, f;
    2399             : 
    2400      393240 :   nf = checknf(nf); checkprid(pr);
    2401      393232 :   if (vT < 0) vT = nf_get_varn(nf);
    2402      393224 :   f = pr_get_f(pr);
    2403      393219 :   N = nf_get_degree(nf);
    2404      393219 :   prh = pr_hnf(nf, pr);
    2405      393241 :   tau = zk? gen_0: anti_uniformizer(nf, pr);
    2406      393223 :   p = pr_get_p(pr);
    2407             : 
    2408      393225 :   if (f == 1)
    2409             :   {
    2410      373261 :     res = cgetg(SMALLMODPR, t_COL);
    2411      373262 :     gel(res,mpr_TAU) = tau;
    2412      373262 :     gel(res,mpr_FFP) = dim1proj(prh);
    2413      373202 :     gel(res,3) = pr; return gerepilecopy(av, res);
    2414             :   }
    2415             : 
    2416       19964 :   c = cgetg(f+1, t_VECSMALL);
    2417       19964 :   ffproj = cgetg(N+1, t_MAT);
    2418      130319 :   for (k=i=1; i<=N; i++)
    2419             :   {
    2420      110355 :     x = gcoeff(prh, i,i);
    2421      110355 :     if (!is_pm1(x)) { c[k] = i; gel(ffproj,i) = col_ei(N, i); k++; }
    2422             :     else
    2423       53886 :       gel(ffproj,i) = ZC_neg(gel(prh,i));
    2424             :   }
    2425       19964 :   ffproj = rowpermute(ffproj, c);
    2426       19964 :   if (! dvdii(nf_get_index(nf), p))
    2427             :   {
    2428       17017 :     GEN basis = nf_get_zkprimpart(nf), D = nf_get_zkden(nf);
    2429       17017 :     if (N == f)
    2430             :     { /* pr inert */
    2431        8568 :       T = nf_get_pol(nf);
    2432        8568 :       T = FpX_red(T,p);
    2433        8568 :       ffproj = RgV_to_RgM(basis, lg(basis)-1);
    2434             :     }
    2435             :     else
    2436             :     {
    2437        8449 :       T = RgV_RgC_mul(basis, pr_get_gen(pr));
    2438        8449 :       T = FpX_normalize(T,p);
    2439        8449 :       basis = FqV_red(vecpermute(basis,c), T, p);
    2440        8449 :       basis = RgV_to_RgM(basis, lg(basis)-1);
    2441        8449 :       ffproj = ZM_mul(basis, ffproj);
    2442             :     }
    2443       17017 :     setvarn(T, vT);
    2444       17017 :     ffproj = FpM_red(ffproj, p);
    2445       17017 :     if (!equali1(D))
    2446             :     {
    2447        2373 :       D = modii(D,p);
    2448        2373 :       if (!equali1(D)) ffproj = FpM_Fp_mul(ffproj, Fp_inv(D,p), p);
    2449             :     }
    2450             : 
    2451       17017 :     res = cgetg(SMALLMODPR+1, t_COL);
    2452       17017 :     gel(res,mpr_TAU) = tau;
    2453       17017 :     gel(res,mpr_FFP) = ffproj;
    2454       17017 :     gel(res,3) = pr;
    2455       17017 :     gel(res,4) = T; return gerepilecopy(av, res);
    2456             :   }
    2457             : 
    2458        2947 :   if (uisprime(f))
    2459             :   {
    2460        2653 :     mul = ei_multable(nf, c[2]);
    2461        2653 :     mul = vecpermute(mul, c);
    2462             :   }
    2463             :   else
    2464             :   {
    2465             :     GEN v, u, u2, frob;
    2466             :     long deg,deg1,deg2;
    2467             : 
    2468             :     /* matrix of Frob: x->x^p over Z_K/pr = < w[c1], ..., w[cf] > over Fp */
    2469         294 :     frob = cgetg(f+1, t_MAT);
    2470        1792 :     for (i=1; i<=f; i++)
    2471             :     {
    2472        1498 :       x = pow_ei_mod_p(nf,c[i],p);
    2473        1498 :       gel(frob,i) = FpM_FpC_mul(ffproj, x, p);
    2474             :     }
    2475         294 :     u = col_ei(f,2); k = 2;
    2476         294 :     deg1 = ffdegree(u, frob, p);
    2477         735 :     while (deg1 < f)
    2478             :     {
    2479         147 :       k++; u2 = col_ei(f, k);
    2480         147 :       deg2 = ffdegree(u2, frob, p);
    2481         147 :       deg = ulcm(deg1,deg2);
    2482         147 :       if (deg == deg1) continue;
    2483         140 :       if (deg == deg2) { deg1 = deg2; u = u2; continue; }
    2484           7 :       u = ZC_add(u, u2);
    2485           7 :       while (ffdegree(u, frob, p) < deg) u = ZC_add(u, u2);
    2486           7 :       deg1 = deg;
    2487             :     }
    2488         294 :     v = lift_to_zk(u,c,N);
    2489             : 
    2490         294 :     mul = cgetg(f+1,t_MAT);
    2491         294 :     gel(mul,1) = v; /* assume w_1 = 1 */
    2492         294 :     for (i=2; i<=f; i++) gel(mul,i) = zk_ei_mul(nf,v,c[i]);
    2493             :   }
    2494             : 
    2495             :   /* Z_K/pr = Fp(v), mul = mul by v */
    2496        2947 :   mul = FpM_red(mul, p);
    2497        2947 :   mul = FpM_mul(ffproj, mul, p);
    2498             : 
    2499        2947 :   pow = get_powers(mul, p);
    2500        2947 :   T = RgV_to_RgX(FpM_deplin(pow, p), vT);
    2501        2947 :   nfproj = cgetg(f+1, t_MAT);
    2502        2947 :   for (i=1; i<=f; i++) gel(nfproj,i) = lift_to_zk(gel(pow,i), c, N);
    2503             : 
    2504        2947 :   setlg(pow, f+1);
    2505        2947 :   ffproj = FpM_mul(FpM_inv(pow, p), ffproj, p);
    2506             : 
    2507        2947 :   res = cgetg(LARGEMODPR, t_COL);
    2508        2947 :   gel(res,mpr_TAU) = tau;
    2509        2947 :   gel(res,mpr_FFP) = ffproj;
    2510        2947 :   gel(res,3) = pr;
    2511        2947 :   gel(res,4) = T;
    2512        2947 :   gel(res,mpr_NFP) = nfproj; return gerepilecopy(av, res);
    2513             : }
    2514             : 
    2515             : GEN
    2516           7 : nfmodprinit(GEN nf, GEN pr) { return modprinit(nf, pr, 0, -1); }
    2517             : GEN
    2518        7023 : zkmodprinit(GEN nf, GEN pr) { return modprinit(nf, pr, 1, -1); }
    2519             : GEN
    2520          56 : nfmodprinit0(GEN nf, GEN pr, long v) { return modprinit(nf, pr, 0, v); }
    2521             : 
    2522             : /* x may be a modpr */
    2523             : static int
    2524     1418063 : ok_modpr(GEN x)
    2525     1418063 : { return typ(x) == t_COL && lg(x) >= SMALLMODPR && lg(x) <= LARGEMODPR; }
    2526             : void
    2527         210 : checkmodpr(GEN x)
    2528             : {
    2529         210 :   if (!ok_modpr(x)) pari_err_TYPE("checkmodpr [use nfmodprinit]", x);
    2530         210 :   checkprid(modpr_get_pr(x));
    2531         210 : }
    2532             : GEN
    2533        3493 : get_modpr(GEN x)
    2534        3493 : { return ok_modpr(x)? x: NULL; }
    2535             : 
    2536             : int
    2537     3148105 : checkprid_i(GEN x)
    2538             : {
    2539     5743040 :   return (typ(x) == t_VEC && lg(x) == 6
    2540     2558303 :           && typ(gel(x,2)) == t_COL && typ(gel(x,3)) == t_INT
    2541     5706366 :           && typ(gel(x,5)) != t_COL); /* tau changed to t_MAT/t_INT in 2.6 */
    2542             : }
    2543             : void
    2544     2349970 : checkprid(GEN x)
    2545     2349970 : { if (!checkprid_i(x)) pari_err_TYPE("checkprid",x); }
    2546             : GEN
    2547      754208 : get_prid(GEN x)
    2548             : {
    2549      754208 :   long lx = lg(x);
    2550      754208 :   if (lx == 3 && typ(x) == t_VEC) x = gel(x,1);
    2551      754208 :   if (checkprid_i(x)) return x;
    2552      550298 :   if (ok_modpr(x)) {
    2553        2989 :     x = modpr_get_pr(x);
    2554        2989 :     if (checkprid_i(x)) return x;
    2555             :   }
    2556      547309 :   return NULL;
    2557             : }
    2558             : 
    2559             : static GEN
    2560      864061 : to_ff_init(GEN nf, GEN *pr, GEN *T, GEN *p, int zk)
    2561             : {
    2562      864061 :   GEN modpr = ok_modpr(*pr)? *pr: modprinit(nf, *pr, zk, -1);
    2563      864072 :   *T = modpr_get_T(modpr);
    2564      864068 :   *pr = modpr_get_pr(modpr);
    2565      864070 :   *p = pr_get_p(*pr); return modpr;
    2566             : }
    2567             : 
    2568             : /* Return an element of O_K which is set to x Mod T */
    2569             : GEN
    2570        4200 : modpr_genFq(GEN modpr)
    2571             : {
    2572        4200 :   switch(lg(modpr))
    2573             :   {
    2574             :     case SMALLMODPR: /* Fp */
    2575         917 :       return gen_1;
    2576             :     case LARGEMODPR:  /* painful case, p \mid index */
    2577        1428 :       return gmael(modpr,mpr_NFP, 2);
    2578             :     default: /* trivial case : p \nmid index */
    2579             :     {
    2580        1855 :       long v = varn( modpr_get_T(modpr) );
    2581        1855 :       return pol_x(v);
    2582             :     }
    2583             :   }
    2584             : }
    2585             : 
    2586             : GEN
    2587      853283 : nf_to_Fq_init(GEN nf, GEN *pr, GEN *T, GEN *p) {
    2588      853283 :   GEN modpr = to_ff_init(nf,pr,T,p,0);
    2589      853288 :   GEN tau = modpr_TAU(modpr);
    2590      853280 :   if (!tau) gel(modpr,mpr_TAU) = anti_uniformizer(nf, *pr);
    2591      853280 :   return modpr;
    2592             : }
    2593             : GEN
    2594       10780 : zk_to_Fq_init(GEN nf, GEN *pr, GEN *T, GEN *p) {
    2595       10780 :   return to_ff_init(nf,pr,T,p,1);
    2596             : }
    2597             : 
    2598             : /* assume x in 'basis' form (t_COL) */
    2599             : GEN
    2600     1283848 : zk_to_Fq(GEN x, GEN modpr)
    2601             : {
    2602     1283848 :   GEN pr = modpr_get_pr(modpr), p = pr_get_p(pr);
    2603     1283885 :   GEN ffproj = gel(modpr,mpr_FFP);
    2604     1283885 :   GEN T = modpr_get_T(modpr);
    2605     1283896 :   return T? FpM_FpC_mul_FpX(ffproj,x, p, varn(T)): FpV_dotproduct(ffproj,x, p);
    2606             : }
    2607             : 
    2608             : /* REDUCTION Modulo a prime ideal */
    2609             : 
    2610             : /* nf a true nf */
    2611             : static GEN
    2612     5699621 : Rg_to_ff(GEN nf, GEN x0, GEN modpr)
    2613             : {
    2614     5699621 :   GEN x = x0, den, pr = modpr_get_pr(modpr), p = pr_get_p(pr);
    2615     5699631 :   long tx = typ(x);
    2616             : 
    2617     5699631 :   if (tx == t_POLMOD) { x = gel(x,2); tx = typ(x); }
    2618     5699631 :   switch(tx)
    2619             :   {
    2620     4387490 :     case t_INT: return modii(x, p);
    2621        6923 :     case t_FRAC: return Rg_to_Fp(x, p);
    2622             :     case t_POL:
    2623      167990 :       switch(lg(x))
    2624             :       {
    2625         231 :         case 2: return gen_0;
    2626       24841 :         case 3: return Rg_to_Fp(gel(x,2), p);
    2627             :       }
    2628      142918 :       x = Q_remove_denom(x, &den);
    2629      142891 :       x = poltobasis(nf, x);
    2630             :       /* content(x) and den may not be coprime */
    2631      142706 :       break;
    2632             :     case t_COL:
    2633     1137246 :       x = Q_remove_denom(x, &den);
    2634             :       /* content(x) and den are coprime */
    2635     1137246 :       if (lg(x)-1 == nf_get_degree(nf)) break;
    2636          48 :     default: pari_err_TYPE("Rg_to_ff",x);
    2637             :       return NULL;/*LCOV_EXCL_LINE*/
    2638             :   }
    2639     1279896 :   if (den)
    2640             :   {
    2641      113943 :     long v = Z_pvalrem(den, p, &den);
    2642      113943 :     if (v)
    2643             :     {
    2644        5628 :       if (tx == t_POL) v -= ZV_pvalrem(x, p, &x);
    2645             :       /* now v = valuation(true denominator of x) */
    2646        5628 :       if (v > 0)
    2647             :       {
    2648        5110 :         GEN tau = modpr_TAU(modpr);
    2649        5110 :         if (!tau) pari_err_TYPE("zk_to_ff", x0);
    2650        5110 :         x = nfmuli(nf,x, nfpow_u(nf, tau, v));
    2651        5110 :         v -= ZV_pvalrem(x, p, &x);
    2652             :       }
    2653        5628 :       if (v > 0) pari_err_INV("Rg_to_ff", mkintmod(gen_0,p));
    2654        5600 :       if (v) return gen_0;
    2655        5292 :       if (is_pm1(den)) den = NULL;
    2656             :     }
    2657      113607 :     x = FpC_red(x, p);
    2658             :   }
    2659     1279560 :   x = zk_to_Fq(x, modpr);
    2660     1279639 :   if (den)
    2661             :   {
    2662      110268 :     GEN c = Fp_inv(den, p);
    2663      110268 :     x = typ(x) == t_INT? Fp_mul(x,c,p): FpX_Fp_mul(x,c,p);
    2664             :   }
    2665     1279639 :   return x;
    2666             : }
    2667             : 
    2668             : GEN
    2669         210 : nfreducemodpr(GEN nf, GEN x, GEN modpr)
    2670             : {
    2671         210 :   pari_sp av = avma;
    2672         210 :   nf = checknf(nf); checkmodpr(modpr);
    2673         210 :   return gerepileupto(av, algtobasis(nf, Fq_to_nf(Rg_to_ff(nf,x,modpr),modpr)));
    2674             : }
    2675             : 
    2676             : GEN
    2677         280 : nfmodpr(GEN nf, GEN x, GEN pr)
    2678             : {
    2679         280 :   pari_sp av = avma;
    2680             :   GEN T, p, modpr;
    2681         280 :   nf = checknf(nf);
    2682         280 :   modpr = nf_to_Fq_init(nf, &pr, &T, &p);
    2683         273 :   x = Rg_to_ff(nf, x, modpr);
    2684         161 :   x = Fq_to_FF(x, Tp_to_FF(T,p));
    2685         161 :   return gerepilecopy(av, x);
    2686             : }
    2687             : GEN
    2688          70 : nfmodprlift(GEN nf, GEN x, GEN pr)
    2689             : {
    2690          70 :   pari_sp av = avma;
    2691             :   GEN y, T, p, modpr;
    2692             :   long i, l, d;
    2693          70 :   nf = checknf(nf);
    2694          70 :   switch(typ(x))
    2695             :   {
    2696           7 :     case t_INT: return icopy(x);
    2697          35 :     case t_FFELT: break;
    2698             :     case t_VEC: case t_COL: case t_MAT:
    2699          28 :       y = cgetg_copy(x,&l);
    2700          28 :       for (i = 1; i < l; i++) gel(y,i) = nfmodprlift(nf,gel(x,i),pr);
    2701          28 :       return y;
    2702           0 :     default: pari_err_TYPE("nfmodprlit",x);
    2703             :   }
    2704          35 :   x = FF_to_FpXQ_i(x);
    2705          35 :   d = degpol(x);
    2706          35 :   if (d <= 0) { set_avma(av); return d? gen_0: icopy(gel(x,2)); }
    2707           7 :   modpr = nf_to_Fq_init(nf, &pr, &T, &p);
    2708           7 :   return gerepilecopy(av, Fq_to_nf(x, modpr));
    2709             : }
    2710             : 
    2711             : /* lift A from residue field to nf */
    2712             : GEN
    2713     1346023 : Fq_to_nf(GEN A, GEN modpr)
    2714             : {
    2715             :   long dA;
    2716     1346023 :   if (typ(A) == t_INT || lg(modpr) < LARGEMODPR) return A;
    2717        5551 :   dA = degpol(A);
    2718        5551 :   if (dA <= 0) return dA ? gen_0: gel(A,2);
    2719        2296 :   return ZM_ZX_mul(gel(modpr,mpr_NFP), A);
    2720             : }
    2721             : GEN
    2722           0 : FqV_to_nfV(GEN x, GEN modpr)
    2723           0 : { pari_APPLY_same(Fq_to_nf(gel(x,i), modpr)) }
    2724             : GEN
    2725        8218 : FqM_to_nfM(GEN A, GEN modpr)
    2726             : {
    2727        8218 :   long i,j,h,l = lg(A);
    2728        8218 :   GEN B = cgetg(l, t_MAT);
    2729             : 
    2730        8218 :   if (l == 1) return B;
    2731        7630 :   h = lgcols(A);
    2732       34972 :   for (j=1; j<l; j++)
    2733             :   {
    2734       27342 :     GEN Aj = gel(A,j), Bj = cgetg(h,t_COL); gel(B,j) = Bj;
    2735       27342 :     for (i=1; i<h; i++) gel(Bj,i) = Fq_to_nf(gel(Aj,i), modpr);
    2736             :   }
    2737        7630 :   return B;
    2738             : }
    2739             : GEN
    2740        9023 : FqX_to_nfX(GEN A, GEN modpr)
    2741             : {
    2742             :   long i, l;
    2743             :   GEN B;
    2744             : 
    2745        9023 :   if (typ(A)!=t_POL) return icopy(A); /* scalar */
    2746        9023 :   B = cgetg_copy(A, &l); B[1] = A[1];
    2747        9023 :   for (i=2; i<l; i++) gel(B,i) = Fq_to_nf(gel(A,i), modpr);
    2748        9023 :   return B;
    2749             : }
    2750             : 
    2751             : /* reduce A to residue field */
    2752             : GEN
    2753     5699130 : nf_to_Fq(GEN nf, GEN A, GEN modpr)
    2754             : {
    2755     5699130 :   pari_sp av = avma;
    2756     5699130 :   return gerepileupto(av, Rg_to_ff(checknf(nf), A, modpr));
    2757             : }
    2758             : /* A t_VEC/t_COL */
    2759             : GEN
    2760        4337 : nfV_to_FqV(GEN A, GEN nf,GEN modpr)
    2761             : {
    2762        4337 :   long i,l = lg(A);
    2763        4337 :   GEN B = cgetg(l,typ(A));
    2764        4337 :   for (i=1; i<l; i++) gel(B,i) = nf_to_Fq(nf,gel(A,i), modpr);
    2765        4337 :   return B;
    2766             : }
    2767             : /* A  t_MAT */
    2768             : GEN
    2769        4361 : nfM_to_FqM(GEN A, GEN nf,GEN modpr)
    2770             : {
    2771        4361 :   long i,j,h,l = lg(A);
    2772        4361 :   GEN B = cgetg(l,t_MAT);
    2773             : 
    2774        4361 :   if (l == 1) return B;
    2775        4361 :   h = lgcols(A);
    2776      130431 :   for (j=1; j<l; j++)
    2777             :   {
    2778      126070 :     GEN Aj = gel(A,j), Bj = cgetg(h,t_COL); gel(B,j) = Bj;
    2779      126070 :     for (i=1; i<h; i++) gel(Bj,i) = nf_to_Fq(nf, gel(Aj,i), modpr);
    2780             :   }
    2781        4361 :   return B;
    2782             : }
    2783             : /* A t_POL */
    2784             : GEN
    2785       12306 : nfX_to_FqX(GEN A, GEN nf,GEN modpr)
    2786             : {
    2787       12306 :   long i,l = lg(A);
    2788       12306 :   GEN B = cgetg(l,t_POL); B[1] = A[1];
    2789       12306 :   for (i=2; i<l; i++) gel(B,i) = nf_to_Fq(nf,gel(A,i),modpr);
    2790       12299 :   return normalizepol_lg(B, l);
    2791             : }
    2792             : 
    2793             : /*******************************************************************/
    2794             : /*                                                                 */
    2795             : /*                       RELATIVE ROUND 2                          */
    2796             : /*                                                                 */
    2797             : /*******************************************************************/
    2798             : /* Shallow functions */
    2799             : /* FIXME: use a bb_field and export the nfX_* routines */
    2800             : static GEN
    2801        3437 : nfX_sub(GEN nf, GEN x, GEN y)
    2802             : {
    2803        3437 :   long i, lx = lg(x), ly = lg(y);
    2804             :   GEN z;
    2805        3437 :   if (ly <= lx) {
    2806        3437 :     z = cgetg(lx,t_POL); z[1] = x[1];
    2807        3437 :     for (i=2; i < ly; i++) gel(z,i) = nfsub(nf,gel(x,i),gel(y,i));
    2808        3437 :     for (   ; i < lx; i++) gel(z,i) = gel(x,i);
    2809        3437 :     z = normalizepol_lg(z, lx);
    2810             :   } else {
    2811           0 :     z = cgetg(ly,t_POL); z[1] = y[1];
    2812           0 :     for (i=2; i < lx; i++) gel(z,i) = nfsub(nf,gel(x,i),gel(y,i));
    2813           0 :     for (   ; i < ly; i++) gel(z,i) = gneg(gel(y,i));
    2814           0 :     z = normalizepol_lg(z, ly);
    2815             :   }
    2816        3437 :   return z;
    2817             : }
    2818             : /* FIXME: quadratic multiplication */
    2819             : static GEN
    2820       55363 : nfX_mul(GEN nf, GEN a, GEN b)
    2821             : {
    2822       55363 :   long da = degpol(a), db = degpol(b), dc, lc, k;
    2823             :   GEN c;
    2824       55363 :   if (da < 0 || db < 0) return gen_0;
    2825       55363 :   dc = da + db;
    2826       55363 :   if (dc == 0) return nfmul(nf, gel(a,2),gel(b,2));
    2827       55363 :   lc = dc+3;
    2828       55363 :   c = cgetg(lc, t_POL); c[1] = a[1];
    2829      443737 :   for (k = 0; k <= dc; k++)
    2830             :   {
    2831      388374 :     long i, I = minss(k, da);
    2832      388374 :     GEN d = NULL;
    2833     1315125 :     for (i = maxss(k-db, 0); i <= I; i++)
    2834             :     {
    2835      926751 :       GEN e = nfmul(nf, gel(a, i+2), gel(b, k-i+2));
    2836      926751 :       d = d? nfadd(nf, d, e): e;
    2837             :     }
    2838      388374 :     gel(c, k+2) = d;
    2839             :   }
    2840       55363 :   return normalizepol_lg(c, lc);
    2841             : }
    2842             : /* assume b monic */
    2843             : static GEN
    2844       51926 : nfX_rem(GEN nf, GEN a, GEN b)
    2845             : {
    2846       51926 :   long da = degpol(a), db = degpol(b);
    2847       51926 :   if (da < 0) return gen_0;
    2848       51926 :   a = leafcopy(a);
    2849      178451 :   while (da >= db)
    2850             :   {
    2851       74599 :     long i, k = da;
    2852       74599 :     GEN A = gel(a, k+2);
    2853      541240 :     for (i = db-1, k--; i >= 0; i--, k--)
    2854      466641 :       gel(a,k+2) = nfsub(nf, gel(a,k+2), nfmul(nf, A, gel(b,i+2)));
    2855       74599 :     a = normalizepol_lg(a, lg(a)-1);
    2856       74599 :     da = degpol(a);
    2857             :   }
    2858       51926 :   return a;
    2859             : }
    2860             : static GEN
    2861       51926 : nfXQ_mul(GEN nf, GEN a, GEN b, GEN T)
    2862             : {
    2863       51926 :   GEN c = nfX_mul(nf, a, b);
    2864       51926 :   if (typ(c) != t_POL) return c;
    2865       51926 :   return nfX_rem(nf, c, T);
    2866             : }
    2867             : 
    2868             : static void
    2869       10444 : fill(long l, GEN H, GEN Hx, GEN I, GEN Ix)
    2870             : {
    2871             :   long i;
    2872       10444 :   if (typ(Ix) == t_VEC) /* standard */
    2873        6335 :     for (i=1; i<l; i++) { gel(H,i) = gel(Hx,i); gel(I,i) = gel(Ix,i); }
    2874             :   else /* constant ideal */
    2875        4109 :     for (i=1; i<l; i++) { gel(H,i) = gel(Hx,i); gel(I,i) = Ix; }
    2876       10444 : }
    2877             : 
    2878             : /* given MODULES x and y by their pseudo-bases, returns a pseudo-basis of the
    2879             :  * module generated by x and y. */
    2880             : static GEN
    2881        5222 : rnfjoinmodules_i(GEN nf, GEN Hx, GEN Ix, GEN Hy, GEN Iy)
    2882             : {
    2883        5222 :   long lx = lg(Hx), ly = lg(Hy), l = lx+ly-1;
    2884        5222 :   GEN H = cgetg(l, t_MAT), I = cgetg(l, t_VEC);
    2885        5222 :   fill(lx, H     , Hx, I     , Ix);
    2886        5222 :   fill(ly, H+lx-1, Hy, I+lx-1, Iy); return nfhnf(nf, mkvec2(H, I));
    2887             : }
    2888             : static GEN
    2889        1806 : rnfjoinmodules(GEN nf, GEN x, GEN y)
    2890             : {
    2891        1806 :   if (!x) return y;
    2892        1113 :   if (!y) return x;
    2893        1113 :   return rnfjoinmodules_i(nf, gel(x,1), gel(x,2), gel(y,1), gel(y,2));
    2894             : }
    2895             : 
    2896             : typedef struct {
    2897             :   GEN multab, T,p;
    2898             :   long h;
    2899             : } rnfeltmod_muldata;
    2900             : 
    2901             : static GEN
    2902       59941 : _sqr(void *data, GEN x)
    2903             : {
    2904       59941 :   rnfeltmod_muldata *D = (rnfeltmod_muldata *) data;
    2905      102025 :   GEN z = x? tablesqr(D->multab,x)
    2906      102025 :            : tablemul_ei_ej(D->multab,D->h,D->h);
    2907       59941 :   return FqV_red(z,D->T,D->p);
    2908             : }
    2909             : static GEN
    2910       10731 : _msqr(void *data, GEN x)
    2911             : {
    2912       10731 :   GEN x2 = _sqr(data, x), z;
    2913       10731 :   rnfeltmod_muldata *D = (rnfeltmod_muldata *) data;
    2914       10731 :   z = tablemul_ei(D->multab, x2, D->h);
    2915       10731 :   return FqV_red(z,D->T,D->p);
    2916             : }
    2917             : 
    2918             : /* Compute W[h]^n mod (T,p) in the extension, assume n >= 0. T a ZX */
    2919             : static GEN
    2920       17857 : rnfeltid_powmod(GEN multab, long h, GEN n, GEN T, GEN p)
    2921             : {
    2922       17857 :   pari_sp av = avma;
    2923             :   GEN y;
    2924             :   rnfeltmod_muldata D;
    2925             : 
    2926       17857 :   if (!signe(n)) return gen_1;
    2927             : 
    2928       17857 :   D.multab = multab;
    2929       17857 :   D.h = h;
    2930       17857 :   D.T = T;
    2931       17857 :   D.p = p;
    2932       17857 :   y = gen_pow_fold(NULL, n, (void*)&D, &_sqr, &_msqr);
    2933       17857 :   return gerepilecopy(av, y);
    2934             : }
    2935             : 
    2936             : /* P != 0 has at most degpol(P) roots. Look for an element in Fq which is not
    2937             :  * a root, cf repres() */
    2938             : static GEN
    2939          21 : FqX_non_root(GEN P, GEN T, GEN p)
    2940             : {
    2941          21 :   long dP = degpol(P), f, vT;
    2942             :   long i, j, k, pi, pp;
    2943             :   GEN v;
    2944             : 
    2945          21 :   if (dP == 0) return gen_1;
    2946          21 :   pp = is_bigint(p) ? dP+1: itos(p);
    2947          21 :   v = cgetg(dP + 2, t_VEC);
    2948          21 :   gel(v,1) = gen_0;
    2949          21 :   if (T)
    2950           0 :   { f = degpol(T); vT = varn(T); }
    2951             :   else
    2952          21 :   { f = 1; vT = 0; }
    2953          42 :   for (i=pi=1; i<=f; i++,pi*=pp)
    2954             :   {
    2955          21 :     GEN gi = i == 1? gen_1: pol_xn(i-1, vT), jgi = gi;
    2956          42 :     for (j=1; j<pp; j++)
    2957             :     {
    2958          42 :       for (k=1; k<=pi; k++)
    2959             :       {
    2960          21 :         GEN z = Fq_add(gel(v,k), jgi, T,p);
    2961          21 :         if (!gequal0(FqX_eval(P, z, T,p))) return z;
    2962          21 :         gel(v, j*pi+k) = z;
    2963             :       }
    2964          21 :       if (j < pp-1) jgi = Fq_add(jgi, gi, T,p); /* j*g[i] */
    2965             :     }
    2966             :   }
    2967          21 :   return NULL;
    2968             : }
    2969             : 
    2970             : /* Relative Dedekind criterion over (true) nf, applied to the order defined by a
    2971             :  * root of monic irreducible polynomial P, modulo the prime ideal pr. Assume
    2972             :  * vdisc = v_pr( disc(P) ).
    2973             :  * Return NULL if nf[X]/P is pr-maximal. Otherwise, return [flag, O, v]:
    2974             :  *   O = enlarged order, given by a pseudo-basis
    2975             :  *   flag = 1 if O is proven pr-maximal (may be 0 and O nevertheless pr-maximal)
    2976             :  *   v = v_pr(disc(O)). */
    2977             : static GEN
    2978        3472 : rnfdedekind_i(GEN nf, GEN P, GEN pr, long vdisc, long only_maximal)
    2979             : {
    2980             :   GEN Ppr, A, I, p, tau, g, h, k, base, T, gzk, hzk, prinvp, pal, nfT, modpr;
    2981             :   long m, vt, r, d, i, j, mpr;
    2982             : 
    2983        3472 :   if (vdisc < 0) pari_err_TYPE("rnfdedekind [non integral pol]", P);
    2984        3465 :   if (vdisc == 1) return NULL; /* pr-maximal */
    2985        3465 :   if (!only_maximal && !gequal1(leading_coeff(P)))
    2986           0 :     pari_err_IMPL( "the full Dedekind criterion in the non-monic case");
    2987             :   /* either monic OR only_maximal = 1 */
    2988        3465 :   m = degpol(P);
    2989        3465 :   nfT = nf_get_pol(nf);
    2990        3465 :   modpr = nf_to_Fq_init(nf,&pr, &T, &p);
    2991        3465 :   Ppr = nfX_to_FqX(P, nf, modpr);
    2992        3458 :   mpr = degpol(Ppr);
    2993        3458 :   if (mpr < m) /* non-monic => only_maximal = 1 */
    2994             :   {
    2995          21 :     if (mpr < 0) return NULL;
    2996          21 :     if (! RgX_valrem(Ppr, &Ppr))
    2997             :     { /* non-zero constant coefficient */
    2998           0 :       Ppr = RgX_shift_shallow(RgX_recip_shallow(Ppr), m - mpr);
    2999           0 :       P = RgX_recip_shallow(P);
    3000             :     }
    3001             :     else
    3002             :     {
    3003          21 :       GEN z = FqX_non_root(Ppr, T, p);
    3004          21 :       if (!z) pari_err_IMPL( "Dedekind in the difficult case");
    3005           0 :       z = Fq_to_nf(z, modpr);
    3006           0 :       if (typ(z) == t_INT)
    3007           0 :         P = RgX_translate(P, z);
    3008             :       else
    3009           0 :         P = RgXQX_translate(P, z, T);
    3010           0 :       P = RgX_recip_shallow(P);
    3011           0 :       Ppr = nfX_to_FqX(P, nf, modpr); /* degpol(P) = degpol(Ppr) = m */
    3012             :     }
    3013             :   }
    3014        3437 :   A = gel(FqX_factor(Ppr,T,p),1);
    3015        3437 :   r = lg(A); /* > 1 */
    3016        3437 :   g = gel(A,1);
    3017        3437 :   for (i=2; i<r; i++) g = FqX_mul(g, gel(A,i), T, p);
    3018        3437 :   h = FqX_div(Ppr,g, T, p);
    3019        3437 :   gzk = FqX_to_nfX(g, modpr);
    3020        3437 :   hzk = FqX_to_nfX(h, modpr);
    3021        3437 :   k = nfX_sub(nf, P, nfX_mul(nf, gzk,hzk));
    3022        3437 :   tau = pr_get_tau(pr);
    3023        3437 :   switch(typ(tau))
    3024             :   {
    3025        1442 :     case t_INT: k = gdiv(k, p); break;
    3026        1995 :     case t_MAT: k = RgX_Rg_div(tablemulvec(NULL,tau, k), p); break;
    3027             :   }
    3028        3437 :   k = nfX_to_FqX(k, nf, modpr);
    3029        3437 :   k = FqX_normalize(FqX_gcd(FqX_gcd(g,h,  T,p), k, T,p), T,p);
    3030        3437 :   d = degpol(k);  /* <= m */
    3031        3437 :   if (!d) return NULL; /* pr-maximal */
    3032        2163 :   if (only_maximal) return gen_0; /* not maximal */
    3033             : 
    3034        2142 :   A = cgetg(m+d+1,t_MAT);
    3035        2142 :   I = cgetg(m+d+1,t_VEC); base = mkvec2(A, I);
    3036             :  /* base[2] temporarily multiplied by p, for the final nfhnfmod,
    3037             :   * which requires integral ideals */
    3038        2142 :   prinvp = pr_inv_p(pr); /* again multiplied by p */
    3039       12635 :   for (j=1; j<=m; j++)
    3040             :   {
    3041       10493 :     gel(A,j) = col_ei(m, j);
    3042       10493 :     gel(I,j) = p;
    3043             :   }
    3044        2142 :   pal = FqX_to_nfX(FqX_div(Ppr,k, T,p), modpr);
    3045        4641 :   for (   ; j<=m+d; j++)
    3046             :   {
    3047        2499 :     gel(A,j) = RgX_to_RgC(pal,m);
    3048        2499 :     gel(I,j) = prinvp;
    3049        2499 :     if (j < m+d) pal = RgXQX_rem(RgX_shift_shallow(pal,1),P,nfT);
    3050             :   }
    3051             :   /* the modulus is integral */
    3052        2142 :   base = nfhnfmod(nf,base, idealmulpowprime(nf, powiu(p,m), pr, utoineg(d)));
    3053        2142 :   gel(base,2) = gdiv(gel(base,2), p); /* cancel the factor p */
    3054        2142 :   vt = vdisc - 2*d;
    3055        2142 :   return mkvec3(vt < 2? gen_1: gen_0, base, stoi(vt));
    3056             : }
    3057             : 
    3058             : /* [L:K] = n */
    3059             : static GEN
    3060         889 : triv_order(long n)
    3061             : {
    3062         889 :   GEN z = cgetg(3, t_VEC);
    3063         889 :   gel(z,1) = matid(n);
    3064         889 :   gel(z,2) = const_vec(n, gen_1); return z;
    3065             : }
    3066             : 
    3067             : /* if flag is set, return gen_1 (resp. gen_0) if the order K[X]/(P)
    3068             :  * is pr-maximal (resp. not pr-maximal). */
    3069             : GEN
    3070          84 : rnfdedekind(GEN nf, GEN P, GEN pr, long flag)
    3071             : {
    3072          84 :   pari_sp av = avma;
    3073             :   GEN z, dP;
    3074             :   long v;
    3075             : 
    3076          84 :   nf = checknf(nf);
    3077          84 :   P = RgX_nffix("rnfdedekind", nf_get_pol(nf), P, 0);
    3078          84 :   dP = RgX_disc(P); P = lift_shallow(P);
    3079          84 :   if (!pr)
    3080             :   {
    3081          21 :     GEN fa = idealfactor(nf, dP);
    3082          21 :     GEN Q = gel(fa,1), E = gel(fa,2);
    3083          21 :     pari_sp av2 = avma;
    3084          21 :     long i, l = lg(Q);
    3085          21 :     for (i = 1; i < l; i++, set_avma(av2))
    3086             :     {
    3087          21 :       v = itos(gel(E,i));
    3088          21 :       if (rnfdedekind_i(nf,P,gel(Q,i),v,1)) { set_avma(av); return gen_0; }
    3089           0 :       set_avma(av2);
    3090             :     }
    3091           0 :     set_avma(av); return gen_1;
    3092             :   }
    3093          63 :   else if (typ(pr) == t_VEC)
    3094             :   { /* flag = 1 is implicit */
    3095          63 :     if (lg(pr) == 1) { set_avma(av); return gen_1; }
    3096          63 :     if (typ(gel(pr,1)) == t_VEC)
    3097             :     { /* list of primes */
    3098          14 :       GEN Q = pr;
    3099          14 :       pari_sp av2 = avma;
    3100          14 :       long i, l = lg(Q);
    3101          14 :       for (i = 1; i < l; i++, set_avma(av2))
    3102             :       {
    3103          14 :         v = nfval(nf, dP, gel(Q,i));
    3104          14 :         if (rnfdedekind_i(nf,P,gel(Q,i),v,1)) { set_avma(av); return gen_0; }
    3105             :       }
    3106           0 :       set_avma(av); return gen_1;
    3107             :     }
    3108             :   }
    3109             :   /* single prime */
    3110          49 :   v = nfval(nf, dP, pr);
    3111          49 :   z = rnfdedekind_i(nf, P, pr, v, flag);
    3112          42 :   if (z)
    3113             :   {
    3114          21 :     if (flag) { set_avma(av); return gen_0; }
    3115          14 :     z = gerepilecopy(av, z);
    3116             :   }
    3117             :   else
    3118             :   {
    3119          21 :     set_avma(av); if (flag) return gen_1;
    3120           7 :     z = cgetg(4, t_VEC);
    3121           7 :     gel(z,1) = gen_1;
    3122           7 :     gel(z,2) = triv_order(degpol(P));
    3123           7 :     gel(z,3) = stoi(v);
    3124             :   }
    3125          21 :   return z;
    3126             : }
    3127             : 
    3128             : static int
    3129       22554 : ideal_is1(GEN x) {
    3130       22554 :   switch(typ(x))
    3131             :   {
    3132        9758 :     case t_INT: return is_pm1(x);
    3133       12033 :     case t_MAT: return RgM_isidentity(x);
    3134             :   }
    3135         763 :   return 0;
    3136             : }
    3137             : 
    3138             : /* return a in ideal A such that v_pr(a) = v_pr(A) */
    3139             : static GEN
    3140       12656 : minval(GEN nf, GEN A, GEN pr)
    3141             : {
    3142       12656 :   GEN ab = idealtwoelt(nf,A), a = gel(ab,1), b = gel(ab,2);
    3143       12656 :   if (nfval(nf,a,pr) > nfval(nf,b,pr)) a = b;
    3144       12656 :   return a;
    3145             : }
    3146             : 
    3147             : /* nf a true nf. Return NULL if power order if pr-maximal */
    3148             : static GEN
    3149        3388 : rnfmaxord(GEN nf, GEN pol, GEN pr, long vdisc)
    3150             : {
    3151        3388 :   pari_sp av = avma, av1;
    3152             :   long i, j, k, n, nn, vpol, cnt, sep;
    3153             :   GEN q, q1, p, T, modpr, W, I, p1;
    3154             :   GEN prhinv, mpi, Id;
    3155             : 
    3156        3388 :   if (DEBUGLEVEL>1) err_printf(" treating %Ps^%ld\n", pr, vdisc);
    3157        3388 :   modpr = nf_to_Fq_init(nf,&pr,&T,&p);
    3158        3388 :   av1 = avma;
    3159        3388 :   p1 = rnfdedekind_i(nf, pol, modpr, vdisc, 0);
    3160        3381 :   if (!p1) return gc_NULL(av);
    3161        2128 :   if (is_pm1(gel(p1,1))) return gerepilecopy(av,gel(p1,2));
    3162         917 :   sep = itos(gel(p1,3));
    3163         917 :   W = gmael(p1,2,1);
    3164         917 :   I = gmael(p1,2,2);
    3165         917 :   gerepileall(av1, 2, &W, &I);
    3166             : 
    3167         917 :   mpi = zk_multable(nf, pr_get_gen(pr));
    3168         917 :   n = degpol(pol); nn = n*n;
    3169         917 :   vpol = varn(pol);
    3170         917 :   q1 = q = pr_norm(pr);
    3171         917 :   while (abscmpiu(q1,n) < 0) q1 = mulii(q1,q);
    3172         917 :   Id = matid(n);
    3173         917 :   prhinv = pr_inv(pr);
    3174         917 :   av1 = avma;
    3175        4361 :   for(cnt=1;; cnt++)
    3176        3444 :   {
    3177        4361 :     GEN I0 = leafcopy(I), W0 = leafcopy(W);
    3178             :     GEN Wa, Winv, Ip, A, MW, MWmod, F, pseudo, C, G;
    3179        4361 :     GEN Tauinv = cgetg(n+1, t_VEC), Tau = cgetg(n+1, t_VEC);
    3180             : 
    3181        4361 :     if (DEBUGLEVEL>1) err_printf("    pass no %ld\n",cnt);
    3182       26579 :     for (j=1; j<=n; j++)
    3183             :     {
    3184             :       GEN tau, tauinv;
    3185       22218 :       if (ideal_is1(gel(I,j)))
    3186             :       {
    3187        9562 :         gel(I,j) = gel(Tau,j) = gel(Tauinv,j) = gen_1;
    3188        9562 :         continue;
    3189             :       }
    3190       12656 :       gel(Tau,j) = tau = minval(nf, gel(I,j), pr);
    3191       12656 :       gel(Tauinv,j) = tauinv = nfinv(nf, tau);
    3192       12656 :       gel(W,j) = nfC_nf_mul(nf, gel(W,j), tau);
    3193       12656 :       gel(I,j) = idealmul(nf, tauinv, gel(I,j)); /* v_pr(I[j]) = 0 */
    3194             :     }
    3195             :     /* W = (Z_K/pr)-basis of O/pr. O = (W0,I0) ~ (W, I) */
    3196             : 
    3197             :    /* compute MW: W_i*W_j = sum MW_k,(i,j) W_k */
    3198        4361 :     Wa = RgM_to_RgXV(W,vpol);
    3199        4361 :     Winv = nfM_inv(nf, W);
    3200        4361 :     MW = cgetg(nn+1, t_MAT);
    3201             :     /* W_1 = 1 */
    3202        4361 :     for (j=1; j<=n; j++) gel(MW, j) = gel(MW, (j-1)*n+1) = gel(Id,j);
    3203       22218 :     for (i=2; i<=n; i++)
    3204       69783 :       for (j=i; j<=n; j++)
    3205             :       {
    3206       51926 :         GEN z = nfXQ_mul(nf, gel(Wa,i), gel(Wa,j), pol);
    3207       51926 :         if (typ(z) != t_POL)
    3208           0 :           z = nfC_nf_mul(nf, gel(Winv,1), z);
    3209             :         else
    3210             :         {
    3211       51926 :           z = RgX_to_RgC(z, lg(Winv)-1);
    3212       51926 :           z = nfM_nfC_mul(nf, Winv, z);
    3213             :         }
    3214       51926 :         gel(MW, (i-1)*n+j) = gel(MW, (j-1)*n+i) = z;
    3215             :       }
    3216             : 
    3217             :     /* compute Ip =  pr-radical [ could use Ker(trace) if q large ] */
    3218        4361 :     MWmod = nfM_to_FqM(MW,nf,modpr);
    3219        4361 :     F = cgetg(n+1, t_MAT); gel(F,1) = gel(Id,1);
    3220        4361 :     for (j=2; j<=n; j++) gel(F,j) = rnfeltid_powmod(MWmod, j, q1, T,p);
    3221        4361 :     Ip = FqM_ker(F,T,p);
    3222        4361 :     if (lg(Ip) == 1) { W = W0; I = I0; break; }
    3223             : 
    3224             :     /* Fill C: W_k A_j = sum_i C_(i,j),k A_i */
    3225        4109 :     A = FqM_to_nfM(FqM_suppl(Ip,T,p), modpr);
    3226        4109 :     for (j = lg(Ip); j<=n; j++) gel(A,j) = nfC_multable_mul(gel(A,j), mpi);
    3227        4109 :     MW = nfM_mul(nf, nfM_inv(nf,A), MW);
    3228        4109 :     C = cgetg(n+1, t_MAT);
    3229       25074 :     for (k=1; k<=n; k++)
    3230             :     {
    3231       20965 :       GEN mek = vecslice(MW, (k-1)*n+1, k*n), Ck;
    3232       20965 :       gel(C,k) = Ck = cgetg(nn+1, t_COL);
    3233      139258 :       for (j=1; j<=n; j++)
    3234             :       {
    3235      118293 :         GEN z = nfM_nfC_mul(nf, mek, gel(A,j));
    3236      118293 :         for (i=1; i<=n; i++) gel(Ck, (j-1)*n+i) = nf_to_Fq(nf,gel(z,i),modpr);
    3237             :       }
    3238             :     }
    3239        4109 :     G = FqM_to_nfM(FqM_ker(C,T,p), modpr);
    3240             : 
    3241        4109 :     pseudo = rnfjoinmodules_i(nf, G,prhinv, Id,I);
    3242             :     /* express W in terms of the power basis */
    3243        4109 :     W = nfM_mul(nf, W, gel(pseudo,1));
    3244        4109 :     I = gel(pseudo,2);
    3245             :     /* restore the HNF property W[i,i] = 1. NB: W upper triangular, with
    3246             :      * W[i,i] = Tau[i] */
    3247       25074 :     for (j=1; j<=n; j++)
    3248       20965 :       if (gel(Tau,j) != gen_1)
    3249             :       {
    3250       11837 :         gel(W,j) = nfC_nf_mul(nf, gel(W,j), gel(Tauinv,j));
    3251       11837 :         gel(I,j) = idealmul(nf, gel(Tau,j), gel(I,j));
    3252             :       }
    3253        4109 :     if (DEBUGLEVEL>3) err_printf(" new order:\n%Ps\n%Ps\n", W, I);
    3254        4109 :     if (sep <= 3 || gequal(I,I0)) break;
    3255             : 
    3256        3444 :     if (gc_needed(av1,2))
    3257             :     {
    3258           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"rnfmaxord");
    3259           0 :       gerepileall(av1,2, &W,&I);
    3260             :     }
    3261             :   }
    3262         917 :   return gerepilecopy(av, mkvec2(W, I));
    3263             : }
    3264             : 
    3265             : GEN
    3266      480500 : Rg_nffix(const char *f, GEN T, GEN c, int lift)
    3267             : {
    3268      480500 :   switch(typ(c))
    3269             :   {
    3270      197181 :     case t_INT: case t_FRAC: return c;
    3271             :     case t_POL:
    3272        4648 :       if (lg(c) >= lg(T)) c = RgX_rem(c,T);
    3273        4648 :       break;
    3274             :     case t_POLMOD:
    3275      278664 :       if (!RgX_equal_var(gel(c,1), T)) pari_err_MODULUS(f, gel(c,1),T);
    3276      278363 :       c = gel(c,2);
    3277      278363 :       switch(typ(c))
    3278             :       {
    3279      243118 :         case t_POL: break;
    3280       35245 :         case t_INT: case t_FRAC: return c;
    3281           0 :         default: pari_err_TYPE(f, c);
    3282             :       }
    3283      243118 :       break;
    3284           7 :     default: pari_err_TYPE(f,c);
    3285             :   }
    3286             :   /* typ(c) = t_POL */
    3287      247766 :   if (varn(c) != varn(T)) pari_err_VAR(f, c,T);
    3288      247759 :   switch(lg(c))
    3289             :   {
    3290        7406 :     case 2: return gen_0;
    3291             :     case 3:
    3292       15666 :       c = gel(c,2); if (is_rational_t(typ(c))) return c;
    3293           0 :       pari_err_TYPE(f,c);
    3294             :   }
    3295      224687 :   RgX_check_QX(c, f);
    3296      224673 :   return lift? c: mkpolmod(c, T);
    3297             : }
    3298             : /* check whether P is a polynomials with coeffs in number field Q[y]/(T) */
    3299             : GEN
    3300      191752 : RgX_nffix(const char *f, GEN T, GEN P, int lift)
    3301             : {
    3302      191752 :   long i, l, vT = varn(T);
    3303      191752 :   GEN Q = cgetg_copy(P, &l);
    3304      191752 :   if (typ(P) != t_POL) pari_err_TYPE(stack_strcat(f," [t_POL expected]"), P);
    3305      191752 :   if (varncmp(varn(P), vT) >= 0) pari_err_PRIORITY(f, P, ">=", vT);
    3306      191738 :   Q[1] = P[1];
    3307      191738 :   for (i=2; i<l; i++) gel(Q,i) = Rg_nffix(f, T, gel(P,i), lift);
    3308      191731 :   return normalizepol_lg(Q, l);
    3309             : }
    3310             : GEN
    3311          28 : RgV_nffix(const char *f, GEN T, GEN P, int lift)
    3312             : {
    3313             :   long i, l;
    3314          28 :   GEN Q = cgetg_copy(P, &l);
    3315          28 :   for (i=1; i<l; i++) gel(Q,i) = Rg_nffix(f, T, gel(P,i), lift);
    3316          21 :   return Q;
    3317             : }
    3318             : 
    3319             : static GEN
    3320        1897 : get_d(GEN nf, GEN d)
    3321             : {
    3322        1897 :   GEN b = idealredmodpower(nf, d, 2, 100000);
    3323        1897 :   return nfmul(nf, d, nfsqr(nf,b));
    3324             : }
    3325             : 
    3326             : static GEN
    3327        3129 : pr_factorback(GEN nf, GEN fa)
    3328             : {
    3329        3129 :   GEN P = gel(fa,1), E = gel(fa,2), z = gen_1;
    3330        3129 :   long i, l = lg(P);
    3331        3129 :   for (i = 1; i < l; i++) z = idealmulpowprime(nf, z, gel(P,i), gel(E,i));
    3332        3129 :   return z;
    3333             : }
    3334             : static GEN
    3335        3129 : pr_factorback_scal(GEN nf, GEN fa)
    3336             : {
    3337        3129 :   GEN D = pr_factorback(nf,fa);
    3338        3129 :   if (typ(D) == t_MAT && RgM_isscalar(D,NULL)) D = gcoeff(D,1,1);
    3339        3129 :   return D;
    3340             : }
    3341             : 
    3342             : /* nf = base field K
    3343             :  * pol= monic polynomial in Z_K[X] defining a relative extension L = K[X]/(pol).
    3344             :  * Returns a pseudo-basis [A,I] of Z_L, set *pD to [D,d] and *pf to the
    3345             :  * index-ideal; rnf is used when lim != 0 and may be NULL */
    3346             : GEN
    3347        1848 : rnfallbase(GEN nf, GEN pol, GEN lim, GEN rnf, GEN *pD, GEN *pf, GEN *pDKP)
    3348             : {
    3349             :   long i, j, jf, l;
    3350             :   GEN fa, E, P, Ef, Pf, z, disc;
    3351             : 
    3352        1848 :   nf = checknf(nf);
    3353        1848 :   if (!gequal1(leading_coeff(pol)))
    3354           7 :     pari_err_IMPL("non-monic relative polynomials in rnfallbase");
    3355        1841 :   disc = nf_to_scalar_or_basis(nf, RgX_disc(pol));
    3356        1841 :   if (lim)
    3357             :   {
    3358             :     GEN rnfeq, zknf, dzknf, U, vU, dA, A, MB, dB, BdB, vj, B, Tabs;
    3359         273 :     GEN D = idealhnf(nf, disc);
    3360         273 :     long rU, m = nf_get_degree(nf), n = degpol(pol), N = n*m;
    3361             :     nfmaxord_t S;
    3362             : 
    3363         273 :     if (typ(lim) == t_INT)
    3364          35 :       P = ZV_union_shallow(nf_get_ramified_primes(nf),
    3365          35 :                            gel(Z_factor_limit(gcoeff(D,1,1), itou(lim)), 1));
    3366             :     else
    3367             :     {
    3368         238 :       P = cgetg_copy(lim, &l);
    3369         714 :       for (i = 1; i < l; i++)
    3370             :       {
    3371         476 :         GEN p = gel(lim,i);
    3372         476 :         if (typ(p) != t_INT) p = pr_get_p(p);
    3373         476 :         gel(P,i) = p;
    3374             :       }
    3375         238 :       P = ZV_sort_uniq(P);
    3376             :     }
    3377         273 :     if (rnf)
    3378             :     {
    3379         224 :       rnfeq = rnf_get_map(rnf);
    3380         224 :       zknf = rnf_get_nfzk(rnf);
    3381             :     }
    3382             :     else
    3383             :     {
    3384          49 :       rnfeq = nf_rnfeq(nf, pol);
    3385          49 :       zknf = nf_nfzk(nf, rnfeq);
    3386             :     }
    3387         273 :     dzknf = gel(zknf,1);
    3388         273 :     if (gequal1(dzknf)) dzknf = NULL;
    3389         273 :     Tabs = gel(rnfeq,1);
    3390         273 :     nfmaxord(&S, mkvec2(Tabs,P), 0);
    3391         273 :     B = RgXV_unscale(S.basis, S.unscale);
    3392         273 :     BdB = Q_remove_denom(B, &dB);
    3393         273 :     MB = RgXV_to_RgM(BdB, N); /* HNF */
    3394             : 
    3395         273 :     vU = cgetg(N+1, t_VEC);
    3396         273 :     vj = cgetg(N+1, t_VECSMALL);
    3397         273 :     gel(vU,1) = U = cgetg(m+1, t_MAT);
    3398         273 :     gel(U,1) = col_ei(N, 1);
    3399         273 :     A = dB? (dzknf? gdiv(dB,dzknf): dB): NULL;
    3400         273 :     if (A && gequal1(A)) A = NULL;
    3401         574 :     for (j = 2; j <= m; j++)
    3402             :     {
    3403         301 :       GEN t = gel(zknf,j);
    3404         301 :       if (A) t = ZX_Z_mul(t, A);
    3405         301 :       gel(U,j) = hnf_solve(MB, RgX_to_RgC(t, N));
    3406             :     }
    3407        1855 :     for (i = 2; i <= N; i++)
    3408             :     {
    3409        1582 :       GEN b = gel(BdB,i);
    3410        1582 :       gel(vU,i) = U = cgetg(m+1, t_MAT);
    3411        1582 :       gel(U,1) = hnf_solve(MB, RgX_to_RgC(b, N));
    3412        3507 :       for (j = 2; j <= m; j++)
    3413             :       {
    3414        1925 :         GEN t = ZX_rem(ZX_mul(b, gel(zknf,j)), Tabs);
    3415        1925 :         if (dzknf) t = gdiv(t, dzknf);
    3416        1925 :         gel(U,j) = hnf_solve(MB, RgX_to_RgC(t, N));
    3417             :       }
    3418             :     }
    3419         273 :     vj[1] = 1; U = gel(vU,1); rU = m;
    3420         707 :     for (i = j = 2; i <= N; i++)
    3421             :     {
    3422         707 :       GEN V = shallowconcat(U, gel(vU,i));
    3423         707 :       if (ZM_rank(V) != rU)
    3424             :       {
    3425         707 :         U = V; rU += m; vj[j++] = i;
    3426         707 :         if (rU == N) break;
    3427             :       }
    3428             :     }
    3429         273 :     if (dB) for(;;)
    3430         322 :     {
    3431         588 :       GEN c = gen_1, H = ZM_hnfmodid(U, dB);
    3432         588 :       long ic = 0;
    3433        5257 :       for (i = 1; i <= N; i++)
    3434        4669 :         if (cmpii(gcoeff(H,i,i), c) > 0) { c = gcoeff(H,i,i); ic = i; }
    3435         588 :       if (!ic) break;
    3436         322 :       vj[j++] = ic;
    3437         322 :       U = shallowconcat(H, gel(vU, ic));
    3438             :     }
    3439         273 :     setlg(vj, j);
    3440         273 :     B = vecpermute(B, vj);
    3441             : 
    3442         273 :     l = lg(B);
    3443         273 :     A = cgetg(l,t_MAT);
    3444        1575 :     for (j = 1; j < l; j++)
    3445             :     {
    3446        1302 :       GEN t = eltabstorel_lift(rnfeq, gel(B,j));
    3447        1302 :       gel(A,j) = Rg_to_RgC(t, n);
    3448             :     }
    3449         273 :     A = RgM_to_nfM(nf, A);
    3450         273 :     A = Q_remove_denom(A, &dA);
    3451         273 :     if (!dA)
    3452             :     { /* order is maximal */
    3453          14 :       z = triv_order(n);
    3454          14 :       if (pf) *pf = gen_1;
    3455             :     }
    3456             :     else
    3457             :     {
    3458             :       GEN fi;
    3459             :       /* the first n columns of A are probably in HNF already */
    3460         259 :       A = shallowconcat(vecslice(A,n+1,lg(A)-1), vecslice(A,1,n));
    3461         259 :       A = mkvec2(A, const_vec(l-1,gen_1));
    3462         259 :       if (DEBUGLEVEL > 2) err_printf("rnfallbase: nfhnf in dim %ld\n", l-1);
    3463         259 :       z = nfhnfmod(nf, A, nfdetint(nf,A));
    3464         259 :       gel(z,2) = gdiv(gel(z,2), dA);
    3465         259 :       fi = idealprod(nf,gel(z,2));
    3466         259 :       D = idealmul(nf, D, idealsqr(nf, fi));
    3467         259 :       if (pf) *pf = idealinv(nf, fi);
    3468             :     }
    3469         273 :     if (RgM_isscalar(D,NULL)) D = gcoeff(D,1,1);
    3470         273 :     if (pDKP) { settyp(S.dKP, t_VEC); *pDKP = S.dKP; }
    3471         273 :     *pD = mkvec2(D, get_d(nf, disc)); return z;
    3472             :   }
    3473        1568 :   pol = lift_shallow(pol);
    3474        1568 :   fa = idealfactor(nf, disc);
    3475        1568 :   P = gel(fa,1); l = lg(P); z = NULL;
    3476        1568 :   E = gel(fa,2);
    3477        1568 :   Pf = cgetg(l, t_COL);
    3478        1568 :   Ef = cgetg(l, t_COL);
    3479        5047 :   for (i = j = jf = 1; i < l; i++)
    3480             :   {
    3481        3486 :     GEN pr = gel(P,i);
    3482        3486 :     long e = itos(gel(E,i));
    3483        3486 :     if (e > 1)
    3484             :     {
    3485        2751 :       GEN vD = rnfmaxord(nf, pol, pr, e);
    3486        2744 :       if (vD)
    3487             :       {
    3488        1806 :         long ef = idealprodval(nf, gel(vD,2), pr);
    3489        1806 :         z = rnfjoinmodules(nf, z, vD);
    3490        1806 :         if (ef) { gel(Pf, jf) = pr; gel(Ef, jf++) = stoi(-ef); }
    3491        1806 :         e += 2 * ef;
    3492             :       }
    3493             :     }
    3494        3479 :     if (e) { gel(P, j) = pr; gel(E, j++) = stoi(e); }
    3495             :   }
    3496        1561 :   setlg(P,j);
    3497        1561 :   setlg(E,j);
    3498        1561 :   if (pDKP)
    3499             :   {
    3500        1505 :     GEN v = cgetg(j, t_VEC);
    3501        1505 :     for (i = 1; i < j; i++) gel(v,i) = pr_get_p(gel(P,i));
    3502        1505 :     *pDKP = ZV_sort_uniq(v);
    3503             :   }
    3504        1561 :   if (pf)
    3505             :   {
    3506        1505 :     setlg(Pf, jf);
    3507        1505 :     setlg(Ef, jf); *pf = pr_factorback_scal(nf, mkmat2(Pf,Ef));
    3508             :   }
    3509        1561 :   *pD = mkvec2(pr_factorback_scal(nf,fa), get_d(nf, disc));
    3510        1561 :   return z? z: triv_order(degpol(pol));
    3511             : }
    3512             : 
    3513             : static GEN
    3514         819 : RgX_to_algX(GEN nf, GEN x)
    3515             : {
    3516             :   long i, l;
    3517         819 :   GEN y = cgetg_copy(x, &l); y[1] = x[1];
    3518         819 :   for (i=2; i<l; i++) gel(y,i) = nf_to_scalar_or_alg(nf, gel(x,i));
    3519         819 :   return y;
    3520             : }
    3521             : 
    3522             : GEN
    3523         819 : nfX_to_monic(GEN nf, GEN T, GEN *pL)
    3524             : {
    3525             :   GEN lT, g, a;
    3526         819 :   long i, l = lg(T);
    3527         819 :   if (l == 2) return pol_0(varn(T));
    3528         819 :   if (l == 3) return pol_1(varn(T));
    3529         819 :   nf = checknf(nf);
    3530         819 :   T = Q_primpart(RgX_to_nfX(nf, T));
    3531         819 :   lT = leading_coeff(T); if (pL) *pL = lT;
    3532         819 :   if (isint1(T)) return T;
    3533         819 :   g = cgetg_copy(T, &l); g[1] = T[1]; a = lT;
    3534         819 :   gel(g, l-1) = gen_1;
    3535         819 :   gel(g, l-2) = gel(T,l-2);
    3536         819 :   if (l == 4) return g;
    3537         819 :   if (typ(lT) == t_INT)
    3538             :   {
    3539         805 :     gel(g, l-3) = gmul(a, gel(T,l-3));
    3540         805 :     for (i = l-4; i > 1; i--) { a = mulii(a,lT); gel(g,i) = gmul(a, gel(T,i)); }
    3541             :   }
    3542             :   else
    3543             :   {
    3544          14 :     gel(g, l-3) = nfmul(nf, a, gel(T,l-3));
    3545          35 :     for (i = l-3; i > 1; i--)
    3546             :     {
    3547          21 :       a = nfmul(nf,a,lT);
    3548          21 :       gel(g,i) = nfmul(nf, a, gel(T,i));
    3549             :     }
    3550             :   }
    3551         819 :   return RgX_to_algX(nf, g);
    3552             : }
    3553             : 
    3554             : GEN
    3555         434 : rnfdisc_factored(GEN nf, GEN pol, GEN *pd)
    3556             : {
    3557             :   long i, j, l;
    3558             :   GEN fa, E, P, disc, lim;
    3559             : 
    3560         434 :   nf = checknf(nf);
    3561         434 :   pol = rnfdisc_get_T(nf, pol, &lim);
    3562         434 :   pol = nfX_to_monic(nf, pol, NULL);
    3563             : 
    3564         434 :   disc = nf_to_scalar_or_basis(nf, RgX_disc(pol));
    3565         434 :   pol = lift_shallow(pol);
    3566         434 :   fa = idealfactor_partial(nf, disc, lim);
    3567         434 :   P = gel(fa,1); l = lg(P);
    3568         434 :   E = gel(fa,2);
    3569        1155 :   for (i = j = 1; i < l; i++)
    3570             :   {
    3571         721 :     long e = itos(gel(E,i));
    3572         721 :     GEN pr = gel(P,i);
    3573         721 :     if (e > 1)
    3574             :     {
    3575         637 :       GEN vD = rnfmaxord(nf, pol, pr, e);
    3576         637 :       if (vD) e += 2*idealprodval(nf, gel(vD,2), pr);
    3577             :     }
    3578         721 :     if (e) { gel(P, j) = pr; gel(E, j++) = stoi(e); }
    3579             :   }
    3580         434 :   if (pd) *pd = get_d(nf, disc);
    3581         434 :   setlg(P, j);
    3582         434 :   setlg(E, j); return fa;
    3583             : }
    3584             : GEN
    3585          63 : rnfdiscf(GEN nf, GEN pol)
    3586             : {
    3587          63 :   pari_sp av = avma;
    3588          63 :   GEN d, fa = rnfdisc_factored(nf, pol, &d);
    3589          63 :   return gerepilecopy(av, mkvec2(pr_factorback_scal(nf,fa), d));
    3590             : }
    3591             : 
    3592             : GEN
    3593          35 : gen_if_principal(GEN bnf, GEN x)
    3594             : {
    3595          35 :   pari_sp av = avma;
    3596          35 :   GEN z = bnfisprincipal0(bnf,x, nf_GEN_IF_PRINCIPAL | nf_FORCE);
    3597          35 :   return isintzero(z)? gc_NULL(av): z;
    3598             : }
    3599             : 
    3600             : static int
    3601          63 : is_pseudo_matrix(GEN O)
    3602             : {
    3603         126 :   return (typ(O) ==t_VEC && lg(O) >= 3
    3604          63 :           && typ(gel(O,1)) == t_MAT
    3605          63 :           && typ(gel(O,2)) == t_VEC
    3606         126 :           && lgcols(O) == lg(gel(O,2)));
    3607             : }
    3608             : 
    3609             : /* given bnf and a pseudo-basis of an order in HNF [A,I], tries to simplify
    3610             :  * the HNF as much as possible. The resulting matrix will be upper triangular
    3611             :  * but the diagonal coefficients will not be equal to 1. The ideals are
    3612             :  * guaranteed to be integral and primitive. */
    3613             : GEN
    3614           0 : rnfsimplifybasis(GEN bnf, GEN x)
    3615             : {
    3616           0 :   pari_sp av = avma;
    3617             :   long i, l;
    3618             :   GEN y, Az, Iz, nf, A, I;
    3619             : 
    3620           0 :   bnf = checkbnf(bnf); nf = bnf_get_nf(bnf);
    3621           0 :   if (!is_pseudo_matrix(x)) pari_err_TYPE("rnfsimplifybasis",x);
    3622           0 :   A = gel(x,1);
    3623           0 :   I = gel(x,2); l = lg(I);
    3624           0 :   y = cgetg(3, t_VEC);
    3625           0 :   Az = cgetg(l, t_MAT); gel(y,1) = Az;
    3626           0 :   Iz = cgetg(l, t_VEC); gel(y,2) = Iz;
    3627           0 :   for (i = 1; i < l; i++)
    3628             :   {
    3629             :     GEN c, d;
    3630           0 :     if (ideal_is1(gel(I,i))) {
    3631           0 :       gel(Iz,i) = gen_1;
    3632           0 :       gel(Az,i) = gel(A,i);
    3633           0 :       continue;
    3634             :     }
    3635             : 
    3636           0 :     gel(Iz,i) = Q_primitive_part(gel(I,i), &c);
    3637           0 :     gel(Az,i) = c? RgC_Rg_mul(gel(A,i),c): gel(A,i);
    3638           0 :     if (c && ideal_is1(gel(Iz,i))) continue;
    3639             : 
    3640           0 :     d = gen_if_principal(bnf, gel(Iz,i));
    3641           0 :     if (d)
    3642             :     {
    3643           0 :       gel(Iz,i) = gen_1;
    3644           0 :       gel(Az,i) = nfC_nf_mul(nf, gel(Az,i), d);
    3645             :     }
    3646             :   }
    3647           0 :   return gerepilecopy(av, y);
    3648             : }
    3649             : 
    3650             : static GEN
    3651          70 : get_order(GEN nf, GEN O, const char *s)
    3652             : {
    3653          70 :   if (typ(O) == t_POL)
    3654           7 :     return rnfpseudobasis(nf, O);
    3655          63 :   if (!is_pseudo_matrix(O)) pari_err_TYPE(s, O);
    3656          63 :   return O;
    3657             : }
    3658             : 
    3659             : GEN
    3660          14 : rnfdet(GEN nf, GEN order)
    3661             : {
    3662          14 :   pari_sp av = avma;
    3663             :   GEN A, I, D;
    3664          14 :   nf = checknf(nf);
    3665          14 :   order = get_order(nf, order, "rnfdet");
    3666          14 :   A = gel(order,1);
    3667          14 :   I = gel(order,2);
    3668          14 :   D = idealmul(nf, nfM_det(nf,A), idealprod(nf,I));
    3669          14 :   return gerepileupto(av, D);
    3670             : }
    3671             : 
    3672             : /* Given two fractional ideals a and b, gives x in a, y in b, z in b^-1,
    3673             :    t in a^-1 such that xt-yz=1. In the present version, z is in Z. */
    3674             : static void
    3675          56 : nfidealdet1(GEN nf, GEN a, GEN b, GEN *px, GEN *py, GEN *pz, GEN *pt)
    3676             : {
    3677             :   GEN x, uv, y, da, db;
    3678             : 
    3679          56 :   a = idealinv(nf,a);
    3680          56 :   a = Q_remove_denom(a, &da);
    3681          56 :   b = Q_remove_denom(b, &db);
    3682          56 :   x = idealcoprime(nf,a,b);
    3683          56 :   uv = idealaddtoone(nf, idealmul(nf,x,a), b);
    3684          56 :   y = gel(uv,2);
    3685          56 :   if (da) x = gmul(x,da);
    3686          56 :   if (db) y = gdiv(y,db);
    3687          56 :   *px = x;
    3688          56 :   *py = y;
    3689          56 :   *pz = db ? negi(db): gen_m1;
    3690          56 :   *pt = nfdiv(nf, gel(uv,1), x);
    3691          56 : }
    3692             : 
    3693             : /* given a pseudo-basis of an order in HNF [A,I] (or [A,I,D,d]), gives an
    3694             :  * n x n matrix (not in HNF) of a pseudo-basis and an ideal vector
    3695             :  * [1,1,...,1,I] such that order = Z_K^(n-1) x I.
    3696             :  * Uses the approximation theorem ==> slow. */
    3697             : GEN
    3698          28 : rnfsteinitz(GEN nf, GEN order)
    3699             : {
    3700          28 :   pari_sp av = avma;
    3701             :   long i, n, l;
    3702             :   GEN A, I, p1;
    3703             : 
    3704          28 :   nf = checknf(nf);
    3705          28 :   order = get_order(nf, order, "rnfsteinitz");
    3706          28 :   A = RgM_to_nfM(nf, gel(order,1));
    3707          28 :   I = leafcopy(gel(order,2)); n=lg(A)-1;
    3708         189 :   for (i=1; i<n; i++)
    3709             :   {
    3710         161 :     GEN c1, c2, b, a = gel(I,i);
    3711         161 :     gel(I,i) = gen_1;
    3712         161 :     if (ideal_is1(a)) continue;
    3713             : 
    3714          56 :     c1 = gel(A,i);
    3715          56 :     c2 = gel(A,i+1);
    3716          56 :     b = gel(I,i+1);
    3717          56 :     if (ideal_is1(b))
    3718             :     {
    3719           0 :       gel(A,i) = c2;
    3720           0 :       gel(A,i+1) = gneg(c1);
    3721           0 :       gel(I,i+1) = a;
    3722             :     }
    3723             :     else
    3724             :     {
    3725          56 :       pari_sp av2 = avma;
    3726             :       GEN x, y, z, t;
    3727          56 :       nfidealdet1(nf,a,b, &x,&y,&z,&t);
    3728          56 :       x = RgC_add(nfC_nf_mul(nf, c1, x), nfC_nf_mul(nf, c2, y));
    3729          56 :       y = RgC_add(nfC_nf_mul(nf, c1, z), nfC_nf_mul(nf, c2, t));
    3730          56 :       gerepileall(av2, 2, &x,&y);
    3731          56 :       gel(A,i) = x;
    3732          56 :       gel(A,i+1) = y;
    3733          56 :       gel(I,i+1) = Q_primitive_part(idealmul(nf,a,b), &p1);
    3734          56 :       if (p1) gel(A,i+1) = nfC_nf_mul(nf, gel(A,i+1), p1);
    3735             :     }
    3736             :   }
    3737          28 :   l = lg(order);
    3738          28 :   p1 = cgetg(l,t_VEC);
    3739          28 :   gel(p1,1) = A;
    3740          28 :   gel(p1,2) = I; for (i=3; i<l; i++) gel(p1,i) = gel(order,i);
    3741          28 :   return gerepilecopy(av, p1);
    3742             : }
    3743             : 
    3744             : /* Given bnf and either an order as output by rnfpseudobasis or a polynomial,
    3745             :  * and outputs a basis if it is free, an n+1-generating set if it is not */
    3746             : GEN
    3747          21 : rnfbasis(GEN bnf, GEN order)
    3748             : {
    3749          21 :   pari_sp av = avma;
    3750             :   long j, n;
    3751             :   GEN nf, A, I, cl, col, a;
    3752             : 
    3753          21 :   bnf = checkbnf(bnf); nf = bnf_get_nf(bnf);
    3754          21 :   order = get_order(nf, order, "rnfbasis");
    3755          21 :   I = gel(order,2); n = lg(I)-1;
    3756          21 :   j=1; while (j<n && ideal_is1(gel(I,j))) j++;
    3757          21 :   if (j<n)
    3758             :   {
    3759           7 :     order = rnfsteinitz(nf,order);
    3760           7 :     I = gel(order,2);
    3761             :   }
    3762          21 :   A = gel(order,1);
    3763          21 :   col= gel(A,n); A = vecslice(A, 1, n-1);
    3764          21 :   cl = gel(I,n);
    3765          21 :   a = gen_if_principal(bnf, cl);
    3766          21 :   if (!a)
    3767             :   {
    3768           7 :     GEN v = idealtwoelt(nf, cl);
    3769           7 :     A = shallowconcat(A, gmul(gel(v,1), col));
    3770           7 :     a = gel(v,2);
    3771             :   }
    3772          21 :   A = shallowconcat(A, nfC_nf_mul(nf, col, a));
    3773          21 :   return gerepilecopy(av, A);
    3774             : }
    3775             : 
    3776             : /* Given bnf and either an order as output by rnfpseudobasis or a polynomial,
    3777             :  * and outputs a basis (not pseudo) in Hermite Normal Form if it exists, zero
    3778             :  * if not
    3779             :  */
    3780             : GEN
    3781           7 : rnfhnfbasis(GEN bnf, GEN order)
    3782             : {
    3783           7 :   pari_sp av = avma;
    3784             :   long j, n;
    3785             :   GEN nf, A, I, a;
    3786             : 
    3787           7 :   bnf = checkbnf(bnf); nf = bnf_get_nf(bnf);
    3788           7 :   order = get_order(nf, order, "rnfbasis");
    3789           7 :   A = gel(order,1); A = RgM_shallowcopy(A);
    3790           7 :   I = gel(order,2); n = lg(A)-1;
    3791          42 :   for (j=1; j<=n; j++)
    3792             :   {
    3793          35 :     if (ideal_is1(gel(I,j))) continue;
    3794          14 :     a = gen_if_principal(bnf, gel(I,j));
    3795          14 :     if (!a) { set_avma(av); return gen_0; }
    3796          14 :     gel(A,j) = nfC_nf_mul(nf, gel(A,j), a);
    3797             :   }
    3798           7 :   return gerepilecopy(av,A);
    3799             : }
    3800             : 
    3801             : static long
    3802           7 : rnfisfree_aux(GEN bnf, GEN order)
    3803             : {
    3804             :   long n, j;
    3805             :   GEN nf, P, I;
    3806             : 
    3807           7 :   bnf = checkbnf(bnf);
    3808           7 :   if (is_pm1( bnf_get_no(bnf) )) return 1;
    3809           0 :   nf = bnf_get_nf(bnf);
    3810           0 :   order = get_order(nf, order, "rnfisfree");
    3811           0 :   I = gel(order,2); n = lg(I)-1;
    3812           0 :   j=1; while (j<=n && ideal_is1(gel(I,j))) j++;
    3813           0 :   if (j>n) return 1;
    3814             : 
    3815           0 :   P = gel(I,j);
    3816           0 :   for (j++; j<=n; j++)
    3817           0 :     if (!ideal_is1(gel(I,j))) P = idealmul(nf,P,gel(I,j));
    3818           0 :   return gequal0( isprincipal(bnf,P) );
    3819             : }
    3820             : 
    3821             : long
    3822           7 : rnfisfree(GEN bnf, GEN order)
    3823           7 : { pari_sp av = avma; return gc_long(av, rnfisfree_aux(bnf,order)); }
    3824             : 
    3825             : /**********************************************************************/
    3826             : /**                                                                  **/
    3827             : /**                   COMPOSITUM OF TWO NUMBER FIELDS                **/
    3828             : /**                                                                  **/
    3829             : /**********************************************************************/
    3830             : static GEN
    3831        2233 : compositum_fix(GEN nf, GEN A)
    3832             : {
    3833             :   int ok;
    3834        2233 :   if (nf)
    3835             :   {
    3836         777 :     long i, l = lg(A);
    3837         777 :     A = shallowcopy(A);
    3838         777 :     for (i=2; i<l; i++) gel(A,i) = basistoalg(nf, gel(A,i));
    3839         777 :     ok = nfissquarefree(nf,A);
    3840             :   }
    3841             :   else
    3842             :   {
    3843        1456 :     A = Q_primpart(A); RgX_check_ZX(A,"polcompositum");
    3844        1456 :     ok = ZX_is_squarefree(A);
    3845             :   }
    3846        2233 :   if (!ok) pari_err_DOMAIN("polcompositum","issquarefree(arg)","=",gen_0,A);
    3847        2226 :   return A;
    3848             : }
    3849             : INLINE long
    3850           7 : nextk(long k) { return k>0 ? -k : 1-k; }
    3851             : 
    3852             : /* modular version */
    3853             : GEN
    3854        1183 : nfcompositum(GEN nf, GEN A, GEN B, long flag)
    3855             : {
    3856        1183 :   pari_sp av = avma;
    3857             :   int same;
    3858             :   long v, k;
    3859             :   GEN C, D, LPRS;
    3860             : 
    3861        1183 :   if (typ(A)!=t_POL) pari_err_TYPE("polcompositum",A);
    3862        1183 :   if (typ(B)!=t_POL) pari_err_TYPE("polcompositum",B);
    3863        1183 :   if (degpol(A)<=0 || degpol(B)<=0) pari_err_CONSTPOL("polcompositum");
    3864        1183 :   v = varn(A);
    3865        1183 :   if (varn(B) != v) pari_err_VAR("polcompositum", A,B);
    3866        1183 :   if (nf)
    3867             :   {
    3868         434 :     nf = checknf(nf);
    3869         427 :     if (varncmp(v,nf_get_varn(nf))>=0) pari_err_PRIORITY("polcompositum", nf, ">=",  v);
    3870             :   }
    3871        1141 :   same = (A == B || RgX_equal(A,B));
    3872        1141 :   A = compositum_fix(nf,A);
    3873        1134 :   if (!same) B = compositum_fix(nf,B);
    3874             : 
    3875        1134 :   D = LPRS = NULL; /* -Wall */
    3876        1134 :   k = same? -1: 1;
    3877        1134 :   if (nf)
    3878             :   {
    3879         392 :     long v0 = fetch_var();
    3880             :     GEN q;
    3881           7 :     for(;; k = nextk(k))
    3882           7 :     {
    3883         399 :       GEN chgvar = deg1pol_shallow(stoi(k),pol_x(v0),v);
    3884         399 :       GEN B1 = poleval(B,chgvar);
    3885         399 :       C = RgX_resultant_all(A,B1,&q);
    3886         399 :       C = gsubst(C,v0,pol_x(v));
    3887         399 :       if (nfissquarefree(nf,C)) break;
    3888             :     }
    3889         392 :     C = lift_if_rational(C);
    3890         392 :     if (flag&1)
    3891             :     {
    3892             :       GEN H0, H1;
    3893         182 :       H0 = gsubst(gel(q,2),v0,pol_x(v));
    3894         182 :       H1 = gsubst(gel(q,3),v0,pol_x(v));
    3895         182 :       if (typ(H0) != t_POL) H0 = scalarpol_shallow(H0,v);
    3896         182 :       if (typ(H1) != t_POL) H1 = scalarpol_shallow(H1,v);
    3897         182 :       H0 = lift_if_rational(H0);
    3898         182 :       H1 = lift_if_rational(H1);
    3899         182 :       LPRS = mkvec2(H0,H1);
    3900             :     }
    3901             :   }
    3902             :   else
    3903             :   {
    3904         742 :     B = leafcopy(B); setvarn(B,fetch_var_higher());
    3905         742 :     C = ZX_ZXY_resultant_all(A, B, &k, (flag&1)? &LPRS: NULL);
    3906         742 :     setvarn(C, v);
    3907             :   }
    3908             :   /* C = Res_Y (A(Y), B(X + kY)) guaranteed squarefree */
    3909        1134 :   if (same)
    3910             :   {
    3911          42 :     D = RgX_rescale(A, stoi(1 - k));
    3912          42 :     C = RgX_div(C, D);
    3913          42 :     if (degpol(C) <= 0)
    3914           0 :       C = mkvec(D);
    3915             :     else
    3916          42 :       C = shallowconcat(nf? gel(nffactor(nf,C),1): ZX_DDF(C), D);
    3917             :   }
    3918        1092 :   else if (flag & 2)
    3919         532 :     C = mkvec(C);
    3920             :   else
    3921         560 :     C = nf? gel(nffactor(nf,C),1): ZX_DDF(C);
    3922        1134 :   gen_sort_inplace(C, (void*)(nf?&cmp_RgX: &cmpii), &gen_cmp_RgX, NULL);
    3923        1134 :   if (flag&1)
    3924             :   { /* a,b,c root of A,B,C = compositum, c = b - k a */
    3925         770 :     long i, l = lg(C);
    3926         770 :     GEN a, b, mH0 = RgX_neg(gel(LPRS,1)), H1 = gel(LPRS,2);
    3927         770 :     setvarn(mH0,v);
    3928         770 :     setvarn(H1,v);
    3929        1561 :     for (i=1; i<l; i++)
    3930             :     {
    3931         791 :       GEN D = gel(C,i);
    3932         791 :       a = RgXQ_mul(mH0, nf? RgXQ_inv(H1,D): QXQ_inv(H1,D), D);
    3933         791 :       b = gadd(pol_x(v), gmulsg(k,a));
    3934         791 :       if (degpol(D) == 1) b = RgX_rem(b,D);
    3935         791 :       gel(C,i) = mkvec4(D, mkpolmod(a,D), mkpolmod(b,D), stoi(-k));
    3936             :     }
    3937             :   }
    3938        1134 :   (void)delete_var();
    3939        1134 :   settyp(C, t_VEC);
    3940        1134 :   if (flag&2) C = gel(C,1);
    3941        1134 :   return gerepilecopy(av, C);
    3942             : }
    3943             : GEN
    3944         749 : polcompositum0(GEN A, GEN B, long flag)
    3945         749 : { return nfcompositum(NULL,A,B,flag); }
    3946             : 
    3947             : GEN
    3948          42 : compositum(GEN pol1,GEN pol2) { return polcompositum0(pol1,pol2,0); }
    3949             : GEN
    3950         490 : compositum2(GEN pol1,GEN pol2) { return polcompositum0(pol1,pol2,1); }
    3951             : 
    3952             : /* Assume A,B irreducible (in particular squarefree) and define linearly
    3953             :  * disjoint extensions: no factorisation needed */
    3954             : GEN
    3955         385 : ZX_compositum_disjoint(GEN A, GEN B)
    3956             : {
    3957         385 :   long k = 1;
    3958         385 :   return ZX_ZXY_rnfequation(A, B, &k);
    3959             : }
    3960             : 
    3961             : static GEN
    3962           7 : lastel(GEN x) { return gel(x, lg(x)-1); }
    3963             : 
    3964             : static GEN
    3965         140 : nfsplitting_composite(GEN P)
    3966             : {
    3967         140 :   GEN F = gel(ZX_factor(P), 1), Q = NULL;
    3968         140 :   long i, n = lg(F)-1;
    3969         280 :   for (i = 1; i <= n; i++)
    3970             :   {
    3971         140 :     GEN Fi = gel(F, i);
    3972         140 :     if (degpol(Fi) == 1) continue;
    3973         126 :     Q = Q ? lastel(compositum(Q, Fi)): Fi;
    3974             :   }
    3975         140 :   return Q ? Q: pol_x(varn(P));
    3976             : }
    3977             : GEN
    3978         147 : nfsplitting(GEN T, GEN D)
    3979             : {
    3980         147 :   pari_sp av = avma;
    3981             :   long d, v;
    3982             :   GEN F, K;
    3983         147 :   T = get_nfpol(T,&K);
    3984         140 :   if (!K)
    3985             :   {
    3986         133 :     if (typ(T) != t_POL) pari_err_TYPE("nfsplitting",T);
    3987         133 :     T = Q_primpart(T);
    3988         133 :     RgX_check_ZX(T,"nfsplitting");
    3989             :   }
    3990         140 :   T = nfsplitting_composite(T);
    3991         140 :   d = degpol(T);
    3992         140 :   if (d<=1) return pol_x(varn(T));
    3993         112 :   if (!K) {
    3994         105 :     if (!isint1(leading_coeff(T))) K = T = polredbest(T,0);
    3995         105 :     K = T;
    3996             :   }
    3997         112 :   if (D)
    3998             :   {
    3999          21 :     if (typ(D) != t_INT || signe(D) < 1) pari_err_TYPE("nfsplitting",D);
    4000             :   }
    4001             :   else
    4002             :   {
    4003          91 :     char *data = stack_strcat(pari_datadir, "/galdata");
    4004          91 :     long dmax = pari_is_dir(data)? 11: 7;
    4005          91 :     D = (d <= dmax)? gel(polgalois(T,DEFAULTPREC), 1): mpfact(d);
    4006             :   }
    4007         112 :   d = itos_or_0(D);
    4008         112 :   v = varn(T);
    4009         112 :   T = leafcopy(T); setvarn(T, fetch_var_higher());
    4010         112 :   for(F = T;;)
    4011          35 :   {
    4012         147 :     GEN P = gel(nffactor(K, F), 1), Q = gel(P,lg(P)-1);
    4013         147 :     if (degpol(gel(P,1)) == degpol(Q)) break;
    4014         119 :     F = rnfequation(K,Q);
    4015         119 :     if (degpol(F) == d) break;
    4016             :   }
    4017         112 :   if (umodiu(D,degpol(F)))
    4018             :   {
    4019           7 :     char *sD = itostr(D);
    4020           7 :     pari_warn(warner,stack_strcat("ignoring incorrect degree bound ",sD));
    4021             :   }
    4022         112 :   (void)delete_var();
    4023         112 :   setvarn(F,v);
    4024         112 :   return gerepilecopy(av, F);
    4025             : }

Generated by: LCOV version 1.13