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.16.2 lcov report (development 29115-f22e516b23) Lines: 2213 2337 94.7 %
Date: 2024-03-28 08:06:56 Functions: 168 172 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; either version 2 of the License, or (at your option) any later
       8             : version. It is distributed in the hope that it will be useful, but WITHOUT
       9             : ANY WARRANTY WHATSOEVER.
      10             : Check the License for details. You should have received a copy of it, along
      11             : with the package; see the file 'COPYING'. If not, write to the Free Software
      12             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      13             : 
      14             : /*******************************************************************/
      15             : /*                                                                 */
      16             : /*                       MAXIMAL ORDERS                            */
      17             : /*                                                                 */
      18             : /*******************************************************************/
      19             : #include "pari.h"
      20             : #include "paripriv.h"
      21             : 
      22             : #define DEBUGLEVEL DEBUGLEVEL_nf
      23             : 
      24             : /* allow p = -1 from factorizations, avoid oo loop on p = 1 */
      25             : static long
      26       12285 : safe_Z_pvalrem(GEN x, GEN p, GEN *z)
      27             : {
      28       12285 :   if (is_pm1(p))
      29             :   {
      30          28 :     if (signe(p) > 0) return gvaluation(x,p); /*error*/
      31          21 :     *z = absi(x); return 1;
      32             :   }
      33       12257 :   return Z_pvalrem(x, p, z);
      34             : }
      35             : /* D an integer, P a ZV, return a factorization matrix for D over P, removing
      36             :  * entries with 0 exponent. */
      37             : static GEN
      38        3584 : fact_from_factors(GEN D, GEN P, long flag)
      39             : {
      40        3584 :   long i, l = lg(P), iq = 1;
      41        3584 :   GEN Q = cgetg(l+1,t_COL);
      42        3584 :   GEN E = cgetg(l+1,t_COL);
      43       15862 :   for (i=1; i<l; i++)
      44             :   {
      45       12285 :     GEN p = gel(P,i);
      46             :     long k;
      47       12285 :     if (flag && !equalim1(p))
      48             :     {
      49          14 :       p = gcdii(p, D);
      50          14 :       if (is_pm1(p)) continue;
      51             :     }
      52       12285 :     k = safe_Z_pvalrem(D, p, &D);
      53       12278 :     if (k) { gel(Q,iq) = p; gel(E,iq) = utoipos(k); iq++; }
      54             :   }
      55        3577 :   D = absi_shallow(D);
      56        3577 :   if (!equali1(D))
      57             :   {
      58         616 :     long k = Z_isanypower(D, &D);
      59         616 :     if (!k) k = 1;
      60         616 :     gel(Q,iq) = D; gel(E,iq) = utoipos(k); iq++;
      61             :   }
      62        3577 :   setlg(Q,iq);
      63        3577 :   setlg(E,iq); return mkmat2(Q,E);
      64             : }
      65             : 
      66             : /* d a t_INT; f a t_MAT factorisation of some t_INT sharing some divisors
      67             :  * with d, or a prime (t_INT). Return a factorization F of d: "primes"
      68             :  * entries in f _may_ be composite, and are included as is in d. */
      69             : static GEN
      70        2090 : update_fact(GEN d, GEN f)
      71             : {
      72             :   GEN P;
      73        2090 :   switch (typ(f))
      74             :   {
      75        2076 :     case t_INT: case t_VEC: case t_COL: return f;
      76          14 :     case t_MAT:
      77          14 :       if (lg(f) == 3) { P = gel(f,1); break; }
      78             :     /*fall through*/
      79             :     default:
      80           7 :       pari_err_TYPE("nfbasis [factorization expected]",f);
      81             :       return NULL;/*LCOV_EXCL_LINE*/
      82             :   }
      83           7 :   return fact_from_factors(d, P, 1);
      84             : }
      85             : 
      86             : /* T = C T0(X/L); C = L^d / lt(T0), d = deg(T)
      87             :  * disc T = C^2(d - 1) L^-(d(d-1)) disc T0 = (L^d / lt(T0)^2)^(d-1) disc T0 */
      88             : static GEN
      89      811608 : set_disc(nfmaxord_t *S)
      90             : {
      91             :   GEN L, dT;
      92             :   long d;
      93      811608 :   if (S->T0 == S->T) return ZX_disc(S->T);
      94      249014 :   d = degpol(S->T0);
      95      249016 :   L = S->unscale;
      96      249016 :   if (typ(L) == t_FRAC && abscmpii(gel(L,1), gel(L,2)) < 0)
      97       11905 :     dT = ZX_disc(S->T); /* more efficient */
      98             :   else
      99             :   {
     100      237111 :     GEN l0 = leading_coeff(S->T0);
     101      237109 :     GEN a = gpowgs(gdiv(gpowgs(L, d), sqri(l0)), d-1);
     102      237109 :     dT = gmul(a, ZX_disc(S->T0)); /* more efficient */
     103             :   }
     104      249002 :   return S->dT = dT;
     105             : }
     106             : 
     107             : /* dT != 0 */
     108             : static GEN
     109      787893 : poldiscfactors_i(GEN T, GEN dT, long flag)
     110             : {
     111             :   GEN U, fa, Z, E, P, Tp;
     112             :   long i, l;
     113             : 
     114      787893 :   fa = absZ_factor_limit_strict(dT, minuu(tridiv_bound(dT), maxprime()), &U);
     115      787879 :   if (!U) return fa;
     116         741 :   Z = mkcol(gel(U,1)); P = gel(fa,1); Tp = NULL;
     117        1594 :   while (lg(Z) != 1)
     118             :   { /* pop and handle last element of Z */
     119         853 :     GEN p = veclast(Z), r;
     120         853 :     setlg(Z, lg(Z)-1);
     121         853 :     if (!Tp) /* first time: p is composite and not a power */
     122         741 :       Tp = ZX_deriv(T);
     123             :     else
     124             :     {
     125         112 :       (void)Z_isanypower(p, &p);
     126         112 :       if ((flag || lgefint(p)==3) && BPSW_psp(p))
     127          89 :       { P = vec_append(P, p); continue; }
     128             :     }
     129         764 :     r = FpX_gcd_check(T, Tp, p);
     130         764 :     if (r)
     131          56 :       Z = shallowconcat(Z, Z_cba(r, diviiexact(p,r)));
     132         708 :     else if (flag)
     133           7 :       P = shallowconcat(P, gel(Z_factor(p),1));
     134             :     else
     135         701 :       P = vec_append(P, p);
     136             :   }
     137         741 :   ZV_sort_inplace(P); l = lg(P); E = cgetg(l, t_COL);
     138        6686 :   for (i = 1; i < l; i++) gel(E,i) = utoipos(Z_pvalrem(dT, gel(P,i), &dT));
     139         741 :   return mkmat2(P,E);
     140             : }
     141             : 
     142             : GEN
     143          42 : poldiscfactors(GEN T, long flag)
     144             : {
     145          42 :   pari_sp av = avma;
     146             :   GEN dT;
     147          42 :   if (typ(T) != t_POL || !RgX_is_ZX(T)) pari_err_TYPE("poldiscfactors",T);
     148          42 :   if (flag < 0 || flag > 1) pari_err_FLAG("poldiscfactors");
     149          42 :   dT = ZX_disc(T);
     150          42 :   if (!signe(dT)) retmkvec2(gen_0, Z_factor(gen_0));
     151          35 :   return gerepilecopy(av, mkvec2(dT, poldiscfactors_i(T, dT, flag)));
     152             : }
     153             : 
     154             : static void
     155      811632 : nfmaxord_check_args(nfmaxord_t *S, GEN T, long flag)
     156             : {
     157      811632 :   GEN dT, L, E, P, fa = NULL;
     158             :   pari_timer t;
     159      811632 :   long l, ty = typ(T);
     160             : 
     161      811632 :   if (DEBUGLEVEL) timer_start(&t);
     162      811632 :   if (ty == t_VEC) {
     163       23730 :     if (lg(T) != 3) pari_err_TYPE("nfmaxord",T);
     164       23730 :     fa = gel(T,2); T = gel(T,1); ty = typ(T);
     165             :   }
     166      811632 :   if (ty != t_POL) pari_err_TYPE("nfmaxord",T);
     167      811632 :   T = Q_primpart(T);
     168      811604 :   if (degpol(T) <= 0) pari_err_CONSTPOL("nfmaxord");
     169      811600 :   RgX_check_ZX(T, "nfmaxord");
     170      811597 :   S->T0 = T;
     171      811597 :   S->T = T = ZX_Q_normalize(T, &L);
     172      811609 :   S->unscale = L;
     173      811609 :   S->dT = dT = set_disc(S);
     174      811585 :   S->certify = 1;
     175      811585 :   if (!signe(dT)) pari_err_IRREDPOL("nfmaxord",T);
     176      811585 :   if (fa)
     177             :   {
     178       23730 :     const long MIN = 100; /* include at least all p < 101 */
     179       23730 :     GEN P0 = NULL, U;
     180       23730 :     S->certify = 0;
     181       23730 :     if (!isint1(L)) fa = update_fact(dT, fa);
     182       23723 :     switch(typ(fa))
     183             :     {
     184         238 :       case t_MAT:
     185         238 :         if (!is_Z_factornon0(fa)) pari_err_TYPE("nfmaxord",fa);
     186         231 :         P0 = gel(fa,1); /* fall through */
     187        3577 :       case t_VEC: case t_COL:
     188        3577 :         if (!P0)
     189             :         {
     190        3346 :           if (!RgV_is_ZV(fa)) pari_err_TYPE("nfmaxord",fa);
     191        3346 :           P0 = fa;
     192             :         }
     193        3577 :         P = gel(absZ_factor_limit_strict(dT, MIN, &U), 1);
     194        3577 :         if (lg(P) != 0) { settyp(P, typ(P0)); P0 = shallowconcat(P0,P); }
     195        3577 :         P0 = ZV_sort_uniq_shallow(P0);
     196        3577 :         fa = fact_from_factors(dT, P0, 0);
     197        3570 :         break;
     198       20132 :       case t_INT:
     199       20132 :         fa = absZ_factor_limit(dT, (signe(fa) <= 0)? 1: maxuu(itou(fa), MIN));
     200       20132 :         break;
     201           7 :       default:
     202           7 :         pari_err_TYPE("nfmaxord",fa);
     203             :     }
     204             :   }
     205             :   else
     206             :   {
     207      787855 :     S->certify = !(flag & nf_PARTIALFACT);
     208      787855 :     fa = poldiscfactors_i(T, dT, 0);
     209             :   }
     210      811550 :   P = gel(fa,1); l = lg(P);
     211      811550 :   E = gel(fa,2);
     212      811550 :   if (l > 1 && is_pm1(gel(P,1)))
     213             :   {
     214          21 :     l--;
     215          21 :     P = vecslice(P, 2, l);
     216          21 :     E = vecslice(E, 2, l);
     217             :   }
     218      811551 :   S->dTP = P;
     219      811551 :   S->dTE = vec_to_vecsmall(E);
     220      811552 :   if (DEBUGLEVEL>2) timer_printf(&t, "disc. factorisation");
     221      811552 : }
     222             : 
     223             : static int
     224      215942 : fnz(GEN x,long j)
     225             : {
     226             :   long i;
     227      694609 :   for (i=1; i<j; i++)
     228      527681 :     if (signe(gel(x,i))) return 0;
     229      166928 :   return 1;
     230             : }
     231             : /* return list u[i], 2 by 2 coprime with the same prime divisors as ab */
     232             : static GEN
     233         266 : get_coprimes(GEN a, GEN b)
     234             : {
     235         266 :   long i, k = 1;
     236         266 :   GEN u = cgetg(3, t_COL);
     237         266 :   gel(u,1) = a;
     238         266 :   gel(u,2) = b;
     239             :   /* u1,..., uk 2 by 2 coprime */
     240         973 :   while (k+1 < lg(u))
     241             :   {
     242         707 :     GEN d, c = gel(u,k+1);
     243         707 :     if (is_pm1(c)) { k++; continue; }
     244        1204 :     for (i=1; i<=k; i++)
     245             :     {
     246         777 :       GEN ui = gel(u,i);
     247         777 :       if (is_pm1(ui)) continue;
     248         441 :       d = gcdii(c, ui);
     249         441 :       if (d == gen_1) continue;
     250         441 :       c = diviiexact(c, d);
     251         441 :       gel(u,i) = diviiexact(ui, d);
     252         441 :       u = vec_append(u, d);
     253             :     }
     254         427 :     gel(u,++k) = c;
     255             :   }
     256        1239 :   for (i = k = 1; i < lg(u); i++)
     257         973 :     if (!is_pm1(gel(u,i))) gel(u,k++) = gel(u,i);
     258         266 :   setlg(u, k); return u;
     259             : }
     260             : 
     261             : /*******************************************************************/
     262             : /*                                                                 */
     263             : /*                            ROUND 4                              */
     264             : /*                                                                 */
     265             : /*******************************************************************/
     266             : typedef struct {
     267             :   /* constants */
     268             :   long pisprime; /* -1: unknown, 1: prime,  0: composite */
     269             :   GEN p, f; /* goal: factor f p-adically */
     270             :   long df;
     271             :   GEN pdf; /* p^df = reduced discriminant of f */
     272             :   long mf; /* */
     273             :   GEN psf, pmf; /* stability precision for f, wanted precision for f */
     274             :   long vpsf; /* v_p(p_f) */
     275             :   /* these are updated along the way */
     276             :   GEN phi; /* a p-integer, in Q[X] */
     277             :   GEN phi0; /* a p-integer, in Q[X] from testb2 / testc2, to be composed with
     278             :              * phi when correct precision is known */
     279             :   GEN chi; /* characteristic polynomial of phi (mod psc) in Z[X] */
     280             :   GEN nu; /* irreducible divisor of chi mod p, in Z[X] */
     281             :   GEN invnu; /* numerator ( 1/ Mod(nu, chi) mod pmr ) */
     282             :   GEN Dinvnu;/* denominator ( ... ) */
     283             :   long vDinvnu; /* v_p(Dinvnu) */
     284             :   GEN prc, psc; /* reduced discriminant of chi, stability precision for chi */
     285             :   long vpsc; /* v_p(p_c) */
     286             :   GEN ns, nsf, precns; /* cached Newton sums for nsf and their precision */
     287             : } decomp_t;
     288             : static GEN maxord_i(decomp_t *S, GEN p, GEN f, long mf, GEN w, long flag);
     289             : static GEN dbasis(GEN p, GEN f, long mf, GEN alpha, GEN U);
     290             : static GEN maxord(GEN p,GEN f,long mf);
     291             : static GEN ZX_Dedekind(GEN F, GEN *pg, GEN p);
     292             : 
     293             : static void
     294         476 : fix_PE(GEN *pP, GEN *pE, long i, GEN u, GEN N)
     295             : {
     296             :   GEN P, E;
     297         476 :   long k, l = lg(u), lP = lg(*pP);
     298             :   pari_sp av;
     299             : 
     300         476 :   *pP = P = shallowconcat(*pP, vecslice(u, 2, l-1));
     301         476 :   *pE = E = vecsmall_lengthen(*pE, lP + l-2);
     302         476 :   gel(P,i) = gel(u,1); av = avma;
     303         476 :   E[i] = Z_pvalrem(N, gel(P,i), &N);
     304         959 :   for (k=lP, lP=lg(P); k < lP; k++) E[k] = Z_pvalrem(N, gel(P,k), &N);
     305         476 :   set_avma(av);
     306         476 : }
     307             : static long
     308      684941 : diag_denomval(GEN M, GEN p)
     309             : {
     310             :   long j, v, l;
     311      684941 :   if (typ(M) != t_MAT) return 0;
     312      462348 :   v = 0; l = lg(M);
     313     2066343 :   for (j=1; j<l; j++)
     314             :   {
     315     1603997 :     GEN t = gcoeff(M,j,j);
     316     1603997 :     if (typ(t) == t_FRAC) v += Z_pval(gel(t,2), p);
     317             :   }
     318      462346 :   return v;
     319             : }
     320             : 
     321             : /* n > 1 is composite, not a pure power, and has no prime divisor < 2^14;
     322             :  * return a BPSW divisor of n and smallest k-th root of largest coprime cofactor */
     323             : static GEN
     324         189 : Z_fac(GEN n)
     325             : {
     326         189 :   GEN p = icopy(n), part = ifac_start(p, 0);
     327             :   long e;
     328         189 :   ifac_next(&part , &p, &e); n = diviiexact(n, powiu(p, e));
     329         189 :   (void)Z_isanypower(n, &n); return mkvec2(p, n);
     330             : }
     331             : 
     332             : /* Warning: data computed for T = ZX_Q_normalize(T0). If S.unscale !=
     333             :  * gen_1, caller must take steps to correct the components if it wishes
     334             :  * to stick to the original T0. Return a vector of p-maximal orders, for
     335             :  * those p s.t p^2 | disc(T) [ = S->dTP ]*/
     336             : static GEN
     337      811626 : get_maxord(nfmaxord_t *S, GEN T0, long flag)
     338             : {
     339             :   GEN P, E;
     340             :   VOLATILE GEN O;
     341             :   VOLATILE long lP, i, k;
     342             : 
     343      811626 :   nfmaxord_check_args(S, T0, flag);
     344      811552 :   P = S->dTP; lP = lg(P);
     345      811552 :   E = S->dTE;
     346      811552 :   O = cgetg(1, t_VEC);
     347     3139423 :   for (i=1; i<lP; i++)
     348             :   {
     349             :     VOLATILE pari_sp av;
     350             :     /* includes the silly case where P[i] = -1 */
     351     2327849 :     if (E[i] <= 1)
     352             :     {
     353     1302802 :       if (S->certify)
     354             :       {
     355     1294973 :         GEN p = gel(P,i);
     356     1294973 :         if (signe(p) > 0 && !BPSW_psp(p))
     357             :         {
     358         189 :           fix_PE(&P, &E, i, Z_fac(p), S->dT);
     359         189 :           lP = lg(P); i--; continue;
     360             :         }
     361             :       }
     362     1302630 :       O = vec_append(O, gen_1); continue;
     363             :     }
     364     1025047 :     av = avma;
     365     1025047 :     pari_CATCH(CATCH_ALL) {
     366         266 :       GEN u, err = pari_err_last();
     367             :       long l;
     368         266 :       switch(err_get_num(err))
     369             :       {
     370         266 :         case e_INV:
     371             :         {
     372         266 :           GEN p, x = err_get_compo(err, 2);
     373         266 :           if (typ(x) == t_INTMOD)
     374             :           { /* caught false prime, update factorization */
     375         266 :             p = gcdii(gel(x,1), gel(x,2));
     376         266 :             u = diviiexact(gel(x,1),p);
     377         266 :             if (DEBUGLEVEL) pari_warn(warner,"impossible inverse: %Ps", x);
     378         266 :             gerepileall(av, 2, &p, &u);
     379             : 
     380         266 :             u = get_coprimes(p, u); l = lg(u);
     381             :             /* no small factors, but often a prime power */
     382         798 :             for (k = 1; k < l; k++) (void)Z_isanypower(gel(u,k), &gel(u,k));
     383         266 :             break;
     384             :           }
     385             :           /* fall through */
     386             :         }
     387             :         case e_PRIME: case e_IRREDPOL:
     388             :         { /* we're here because we failed BPSW_isprime(), no point in
     389             :            * reporting a possible counter-example to the BPSW test */
     390           0 :           GEN p = gel(P,i);
     391           0 :           set_avma(av);
     392           0 :           if (DEBUGLEVEL)
     393           0 :             pari_warn(warner,"large composite in nfmaxord:loop(), %Ps", p);
     394           0 :           if (expi(p) < 100)
     395           0 :             u = gel(Z_factor(p), 1); /* p < 2^100 should take ~20ms */
     396           0 :           else if (S->certify)
     397           0 :             u = Z_fac(p);
     398             :           else
     399             :           { /* give up, probably not maximal */
     400           0 :             GEN B, g, k = ZX_Dedekind(S->T, &g, p);
     401           0 :             k = FpX_normalize(k, p);
     402           0 :             B = dbasis(p, S->T, E[i], NULL, FpX_div(S->T,k,p));
     403           0 :             O = vec_append(O, B);
     404           0 :             pari_CATCH_reset(); continue;
     405             :           }
     406           0 :           break;
     407             :         }
     408           0 :         default: pari_err(0, err);
     409             :           return NULL;/*LCOV_EXCL_LINE*/
     410             :       }
     411         266 :       fix_PE(&P, &E, i, u, S->dT);
     412         266 :       lP = lg(P); av = avma;
     413     1025318 :     } pari_RETRY {
     414     1025318 :       GEN p = gel(P,i), O2;
     415     1025318 :       if (DEBUGLEVEL>2) err_printf("Treating p^k = %Ps^%ld\n",p,E[i]);
     416     1025318 :       O2 = maxord(p,S->T,E[i]);
     417     1025054 :       if (S->certify && (odd(E[i]) || E[i] != 2*diag_denomval(O2, p))
     418      609353 :                      && !BPSW_psp(p))
     419             :       {
     420          21 :         fix_PE(&P, &E, i, gel(Z_factor(p), 1), S->dT);
     421          21 :         lP = lg(P); i--;
     422             :       }
     423             :       else
     424     1025039 :         O = vec_append(O, O2);
     425     1025061 :     } pari_ENDCATCH;
     426             :   }
     427      811574 :   S->dTP = P; S->dTE = E; return O;
     428             : }
     429             : 
     430             : /* M a QM, return denominator of diagonal. All denominators are powers of
     431             :  * a given integer */
     432             : static GEN
     433       97325 : diag_denom(GEN M)
     434             : {
     435       97325 :   GEN d = gen_1;
     436       97325 :   long j, l = lg(M);
     437      675844 :   for (j=1; j<l; j++)
     438             :   {
     439      578519 :     GEN t = gcoeff(M,j,j);
     440      578519 :     if (typ(t) == t_INT) continue;
     441      206672 :     t = gel(t,2);
     442      206672 :     if (abscmpii(t,d) > 0) d = t;
     443             :   }
     444       97325 :   return d;
     445             : }
     446             : static void
     447      745148 : setPE(GEN D, GEN P, GEN *pP, GEN *pE)
     448             : {
     449      745148 :   long k, j, l = lg(P);
     450             :   GEN P2, E2;
     451      745148 :   *pP = P2 = cgetg(l, t_VEC);
     452      745151 :   *pE = E2 = cgetg(l, t_VECSMALL);
     453     2865284 :   for (k = j = 1; j < l; j++)
     454             :   {
     455     2120093 :     long v = Z_pvalrem(D, gel(P,j), &D);
     456     2120105 :     if (v) { gel(P2,k) = gel(P,j); E2[k] = v; k++; }
     457             :   }
     458      745191 :   setlg(P2, k);
     459      745187 :   setlg(E2, k);
     460      745183 : }
     461             : void
     462      100430 : nfmaxord(nfmaxord_t *S, GEN T0, long flag)
     463             : {
     464      100430 :   GEN O = get_maxord(S, T0, flag);
     465      100431 :   GEN f = S->T, P = S->dTP, a = NULL, da = NULL;
     466      100431 :   long n = degpol(f), lP = lg(P), i, j, k;
     467      100432 :   int centered = 0;
     468      100432 :   pari_sp av = avma;
     469             :   /* r1 & basden not initialized here */
     470      100432 :   S->r1 = -1;
     471      100432 :   S->basden = NULL;
     472      351489 :   for (i=1; i<lP; i++)
     473             :   {
     474      251058 :     GEN M, db, b = gel(O,i);
     475      251058 :     if (b == gen_1) continue;
     476       97325 :     db = diag_denom(b);
     477       97325 :     if (db == gen_1) continue;
     478             : 
     479             :     /* db = denom(b), (da,db) = 1. Compute da Im(b) + db Im(a) */
     480       97325 :     b = Q_muli_to_int(b,db);
     481       97325 :     if (!da) { da = db; a = b; }
     482             :     else
     483             :     { /* optimization: easy as long as both matrix are diagonal */
     484      131027 :       j=2; while (j<=n && fnz(gel(a,j),j) && fnz(gel(b,j),j)) j++;
     485       49028 :       k = j-1; M = cgetg(2*n-k+1,t_MAT);
     486      180055 :       for (j=1; j<=k; j++)
     487             :       {
     488      131026 :         gel(M,j) = gel(a,j);
     489      131026 :         gcoeff(M,j,j) = mulii(gcoeff(a,j,j),gcoeff(b,j,j));
     490             :       }
     491             :       /* could reduce mod M(j,j) but not worth it: usually close to da*db */
     492      269103 :       for (  ; j<=n;     j++) gel(M,j) = ZC_Z_mul(gel(a,j), db);
     493      269102 :       for (  ; j<=2*n-k; j++) gel(M,j) = ZC_Z_mul(gel(b,j+k-n), da);
     494       49026 :       da = mulii(da,db);
     495       49028 :       a = ZM_hnfmodall_i(M, da, hnf_MODID|hnf_CENTER);
     496       49027 :       gerepileall(av, 2, &a, &da);
     497       49027 :       centered = 1;
     498             :     }
     499             :   }
     500      100431 :   if (da)
     501             :   {
     502       48298 :     GEN index = diviiexact(da, gcoeff(a,1,1));
     503      227393 :     for (j=2; j<=n; j++) index = mulii(index, diviiexact(da, gcoeff(a,j,j)));
     504       48297 :     if (!centered) a = ZM_hnfcenter(a);
     505       48297 :     a = RgM_Rg_div(a, da);
     506       48298 :     S->index = index;
     507       48298 :     S->dK = diviiexact(S->dT, sqri(index));
     508             :   }
     509             :   else
     510             :   {
     511       52133 :     S->index = gen_1;
     512       52133 :     S->dK = S->dT;
     513       52133 :     a = matid(n);
     514             :   }
     515      100428 :   setPE(S->dK, P, &S->dKP, &S->dKE);
     516      100432 :   S->basis = RgM_to_RgXV(a, varn(f));
     517      100430 : }
     518             : GEN
     519         938 : nfbasis(GEN x, GEN *pdK)
     520             : {
     521         938 :   pari_sp av = avma;
     522             :   nfmaxord_t S;
     523             :   GEN B;
     524         938 :   nfmaxord(&S, x, 0);
     525         938 :   B = RgXV_unscale(S.basis, S.unscale);
     526         938 :   if (pdK) *pdK = S.dK;
     527         938 :   return gc_all(av, pdK? 2: 1, &B, pdK);
     528             : }
     529             : /* field discriminant: faster than nfmaxord, use local data only */
     530             : static GEN
     531      711195 : maxord_disc(nfmaxord_t *S, GEN x)
     532             : {
     533      711195 :   GEN O = get_maxord(S, x, 0), I = gen_1;
     534      711167 :   long n = degpol(S->T), lP = lg(O), i, j;
     535     2787758 :   for (i = 1; i < lP; i++)
     536             :   {
     537     2076603 :     GEN b = gel(O,i);
     538     2076603 :     if (b == gen_1) continue;
     539     2706477 :     for (j = 1; j <= n; j++)
     540             :     {
     541     2115123 :       GEN c = gcoeff(b,j,j);
     542     2115123 :       if (typ(c) == t_FRAC) I = mulii(I, gel(c,2)) ;
     543             :     }
     544             :   }
     545      711155 :   return diviiexact(S->dT, sqri(I));
     546             : }
     547             : GEN
     548       66429 : nfdisc(GEN x)
     549             : {
     550       66429 :   pari_sp av = avma;
     551             :   nfmaxord_t S;
     552       66429 :   return gerepileuptoint(av, maxord_disc(&S, x));
     553             : }
     554             : GEN
     555      644773 : nfdiscfactors(GEN x)
     556             : {
     557      644773 :   pari_sp av = avma;
     558      644773 :   GEN E, P, D, nf = checknf_i(x);
     559      644772 :   if (nf)
     560             :   {
     561           7 :     D = nf_get_disc(nf);
     562           7 :     P = nf_get_ramified_primes(nf);
     563             :   }
     564             :   else
     565             :   {
     566             :     nfmaxord_t S;
     567      644765 :     D = maxord_disc(&S, x);
     568      644703 :     P = S.dTP;
     569             :   }
     570      644710 :   setPE(D, P, &P, &E); settyp(P, t_COL);
     571      644750 :   return gerepilecopy(av, mkvec2(D, mkmat2(P, zc_to_ZC(E))));
     572             : }
     573             : 
     574             : static ulong
     575     1593662 : Flx_checkdeflate(GEN x)
     576             : {
     577     1593662 :   ulong d = 0, i, lx = (ulong)lg(x);
     578     2541975 :   for (i=3; i<lx; i++)
     579     1697439 :     if (x[i]) { d = ugcd(d,i-2); if (d == 1) break; }
     580     1593661 :   return d;
     581             : }
     582             : 
     583             : /* product of (monic) irreducible factors of f over Fp[X]
     584             :  * Assume f reduced mod p, otherwise valuation at x may be wrong */
     585             : static GEN
     586     1593666 : Flx_radical(GEN f, ulong p)
     587             : {
     588     1593666 :   long v0 = Flx_valrem(f, &f);
     589             :   ulong du, d, e;
     590             :   GEN u;
     591             : 
     592     1593660 :   d = Flx_checkdeflate(f);
     593     1593715 :   if (!d) return v0? polx_Flx(f[1]): pol1_Flx(f[1]);
     594      997709 :   if (u_lvalrem(d,p, &e)) f = Flx_deflate(f, d/e); /* f(x^p^i) -> f(x) */
     595      997716 :   u = Flx_gcd(f, Flx_deriv(f, p), p); /* (f,f') */
     596      997696 :   du = degpol(u);
     597      997695 :   if (du)
     598             :   {
     599      313338 :     if (du == (ulong)degpol(f))
     600           0 :       f = Flx_radical(Flx_deflate(f,p), p);
     601             :     else
     602             :     {
     603      313338 :       u = Flx_normalize(u, p);
     604      313344 :       f = Flx_div(f, u, p);
     605      313342 :       if (p <= du)
     606             :       {
     607       66250 :         GEN w = (degpol(f) >= degpol(u))? Flx_rem(f, u, p): f;
     608       66250 :         w = Flxq_powu(w, du, u, p);
     609       66250 :         w = Flx_div(u, Flx_gcd(w,u,p), p); /* u / gcd(u, v^(deg u-1)) */
     610       66251 :         f = Flx_mul(f, Flx_radical(Flx_deflate(w,p), p), p);
     611             :       }
     612             :     }
     613             :   }
     614      997712 :   if (v0) f = Flx_shift(f, 1);
     615      997708 :   return f;
     616             : }
     617             : /* Assume f reduced mod p, otherwise valuation at x may be wrong */
     618             : static GEN
     619        5486 : FpX_radical(GEN f, GEN p)
     620             : {
     621             :   GEN u;
     622             :   long v0;
     623        5486 :   if (lgefint(p) == 3)
     624             :   {
     625        1733 :     ulong q = p[2];
     626        1733 :     return Flx_to_ZX( Flx_radical(ZX_to_Flx(f, q), q) );
     627             :   }
     628        3753 :   v0 = ZX_valrem(f, &f);
     629        3753 :   u = FpX_gcd(f,FpX_deriv(f, p), p);
     630        3493 :   if (degpol(u)) f = FpX_div(f, u, p);
     631        3493 :   if (v0) f = RgX_shift(f, 1);
     632        3493 :   return f;
     633             : }
     634             : /* f / a */
     635             : static GEN
     636     1525725 : zx_z_div(GEN f, ulong a)
     637             : {
     638     1525725 :   long i, l = lg(f);
     639     1525725 :   GEN g = cgetg(l, t_VECSMALL);
     640     1525713 :   g[1] = f[1];
     641     5154527 :   for (i = 2; i < l; i++) g[i] = f[i] / a;
     642     1525713 :   return g;
     643             : }
     644             : /* Dedekind criterion; return k = gcd(g,h, (f-gh)/p), where
     645             :  *   f = \prod f_i^e_i, g = \prod f_i, h = \prod f_i^{e_i-1}
     646             :  * k = 1 iff Z[X]/(f) is p-maximal */
     647             : static GEN
     648     1531259 : ZX_Dedekind(GEN F, GEN *pg, GEN p)
     649             : {
     650             :   GEN k, h, g, f, f2;
     651     1531259 :   ulong q = p[2];
     652     1531259 :   if (lgefint(p) == 3 && q < (1UL << BITS_IN_HALFULONG))
     653     1525649 :   {
     654     1525776 :     ulong q2 = q*q;
     655     1525776 :     f2 = ZX_to_Flx(F, q2);
     656     1525750 :     f = Flx_red(f2, q);
     657     1525688 :     g = Flx_radical(f, q);
     658     1525722 :     h = Flx_div(f, g, q);
     659     1525702 :     k = zx_z_div(Flx_sub(f2, Flx_mul(g,h,q2), q2), q);
     660     1525718 :     k = Flx_gcd(k, Flx_gcd(g,h,q), q);
     661     1525707 :     k = Flx_to_ZX(k);
     662     1525642 :     g = Flx_to_ZX(g);
     663             :   }
     664             :   else
     665             :   {
     666        5483 :     f2 = FpX_red(F, sqri(p));
     667        5486 :     f = FpX_red(f2, p);
     668        5486 :     g = FpX_radical(f, p);
     669        5220 :     h = FpX_div(f, g, p);
     670        5220 :     k = ZX_Z_divexact(ZX_sub(f2, ZX_mul(g,h)), p);
     671        5220 :     k = FpX_gcd(FpX_red(k, p), FpX_gcd(g,h,p), p);
     672             :   }
     673     1530873 :   *pg = g; return k;
     674             : }
     675             : 
     676             : /* p-maximal order of Z[x]/f; mf = v_p(Disc(f)) or < 0 [unknown].
     677             :  * Return gen_1 if p-maximal */
     678             : static GEN
     679     1531258 : maxord(GEN p, GEN f, long mf)
     680             : {
     681     1531258 :   const pari_sp av = avma;
     682     1531258 :   GEN res, g, k = ZX_Dedekind(f, &g, p);
     683     1530874 :   long dk = degpol(k);
     684     1530869 :   if (DEBUGLEVEL>2) err_printf("  ZX_Dedekind: gcd has degree %ld\n", dk);
     685     1530933 :   if (!dk) { set_avma(av); return gen_1; }
     686      870120 :   if (mf < 0) mf = ZpX_disc_val(f, p);
     687      870123 :   k = FpX_normalize(k, p);
     688      870128 :   if (2*dk >= mf-1)
     689      418830 :     res = dbasis(p, f, mf, NULL, FpX_div(f,k,p));
     690             :   else
     691             :   {
     692             :     GEN w, F1, F2;
     693             :     decomp_t S;
     694      451298 :     F1 = FpX_factor(k,p);
     695      451331 :     F2 = FpX_factor(FpX_div(g,k,p),p);
     696      451329 :     w = merge_sort_uniq(gel(F1,1),gel(F2,1),(void*)cmpii,&gen_cmp_RgX);
     697      451325 :     res = maxord_i(&S, p, f, mf, w, 0);
     698             :   }
     699      870178 :   return gerepilecopy(av,res);
     700             : }
     701             : /* T monic separable ZX, p prime */
     702             : GEN
     703           0 : ZpX_primedec(GEN T, GEN p)
     704             : {
     705           0 :   const pari_sp av = avma;
     706           0 :   GEN w, F1, F2, res, g, k = ZX_Dedekind(T, &g, p);
     707             :   decomp_t S;
     708           0 :   if (!degpol(k)) return zm_to_ZM(FpX_degfact(T, p));
     709           0 :   k = FpX_normalize(k, p);
     710           0 :   F1 = FpX_factor(k,p);
     711           0 :   F2 = FpX_factor(FpX_div(g,k,p),p);
     712           0 :   w = merge_sort_uniq(gel(F1,1),gel(F2,1),(void*)cmpii,&gen_cmp_RgX);
     713           0 :   res = maxord_i(&S, p, T, ZpX_disc_val(T, p), w, -1);
     714           0 :   if (!res)
     715             :   {
     716           0 :     long f = degpol(S.nu), e = degpol(T) / f;
     717           0 :     set_avma(av); retmkmat2(mkcols(f), mkcols(e));
     718             :   }
     719           0 :   return gerepilecopy(av,res);
     720             : }
     721             : 
     722             : static GEN
     723     4650486 : Zlx_sylvester_echelon(GEN f1, GEN f2, long early_abort, ulong p, ulong pm)
     724             : {
     725     4650486 :   long j, n = degpol(f1);
     726     4650486 :   GEN h, a = cgetg(n+1,t_MAT);
     727     4650466 :   f1 = Flx_get_red(f1, pm);
     728     4650481 :   h = Flx_rem(f2,f1,pm);
     729    16311136 :   for (j=1;; j++)
     730             :   {
     731    16311136 :     gel(a,j) = Flx_to_Flv(h, n);
     732    16310332 :     if (j == n) break;
     733    11660029 :     h = Flx_rem(Flx_shift(h, 1), f1, pm);
     734             :   }
     735     4650303 :   return zlm_echelon(a, early_abort, p, pm);
     736             : }
     737             : /* Sylvester's matrix, mod p^m (assumes f1 monic). If early_abort
     738             :  * is set, return NULL if one pivot is 0 mod p^m */
     739             : static GEN
     740       72926 : ZpX_sylvester_echelon(GEN f1, GEN f2, long early_abort, GEN p, GEN pm)
     741             : {
     742       72926 :   long j, n = degpol(f1);
     743       72926 :   GEN h, a = cgetg(n+1,t_MAT);
     744       72926 :   h = FpXQ_red(f2,f1,pm);
     745      416085 :   for (j=1;; j++)
     746             :   {
     747      416085 :     gel(a,j) = RgX_to_RgC(h, n);
     748      416085 :     if (j == n) break;
     749      343159 :     h = FpX_rem(RgX_shift_shallow(h, 1), f1, pm);
     750             :   }
     751       72926 :   return ZpM_echelon(a, early_abort, p, pm);
     752             : }
     753             : 
     754             : /* polynomial gcd mod p^m (assumes f1 monic). Return a QpX ! */
     755             : static GEN
     756      245079 : Zlx_gcd(GEN f1, GEN f2, ulong p, ulong pm)
     757             : {
     758      245079 :   pari_sp av = avma;
     759      245079 :   GEN a = Zlx_sylvester_echelon(f1,f2,0,p,pm);
     760      245079 :   long c, l = lg(a), sv = f1[1];
     761      748845 :   for (c = 1; c < l; c++)
     762             :   {
     763      748845 :     ulong t = ucoeff(a,c,c);
     764      748845 :     if (t)
     765             :     {
     766      245079 :       a = Flx_to_ZX(Flv_to_Flx(gel(a,c), sv));
     767      245078 :       if (t == 1) return gerepilecopy(av, a);
     768       74787 :       return gerepileupto(av, RgX_Rg_div(a, utoipos(t)));
     769             :     }
     770             :   }
     771           0 :   set_avma(av);
     772           0 :   a = cgetg(2,t_POL); a[1] = sv; return a;
     773             : }
     774             : GEN
     775      253261 : ZpX_gcd(GEN f1, GEN f2, GEN p, GEN pm)
     776             : {
     777      253261 :   pari_sp av = avma;
     778             :   GEN a;
     779             :   long c, l, v;
     780      253261 :   if (lgefint(pm) == 3)
     781             :   {
     782      245081 :     ulong q = pm[2];
     783      245081 :     return Zlx_gcd(ZX_to_Flx(f1, q), ZX_to_Flx(f2,q), p[2], q);
     784             :   }
     785        8180 :   a = ZpX_sylvester_echelon(f1,f2,0,p,pm);
     786        8180 :   l = lg(a); v = varn(f1);
     787       51330 :   for (c = 1; c < l; c++)
     788             :   {
     789       51330 :     GEN t = gcoeff(a,c,c);
     790       51330 :     if (signe(t))
     791             :     {
     792        8180 :       a = RgV_to_RgX(gel(a,c), v);
     793        8180 :       if (equali1(t)) return gerepilecopy(av, a);
     794        2367 :       return gerepileupto(av, RgX_Rg_div(a, t));
     795             :     }
     796             :   }
     797           0 :   set_avma(av); return pol_0(v);
     798             : }
     799             : 
     800             : /* Return m > 0, such that p^m ~ 2^16 for initial value of m; assume p prime */
     801             : static long
     802     4366097 : init_m(GEN p)
     803             : {
     804             :   ulong pp;
     805     4366097 :   if (lgefint(p) > 3) return 1;
     806     4364924 :   pp = p[2]; /* m ~ 16 / log2(pp) */
     807     4364924 :   if (pp < 41) switch(pp)
     808             :   {
     809     1184026 :     case 2: return 16;
     810      350245 :     case 3: return 10;
     811      244853 :     case 5: return 6;
     812      150458 :     case 7: return 5;
     813      207145 :     case 11: case 13: return 4;
     814      298894 :     default: return 3;
     815             :   }
     816     1929303 :   return pp < 257? 2: 1;
     817             : }
     818             : 
     819             : /* reduced resultant mod p^m (assumes x monic) */
     820             : GEN
     821      985828 : ZpX_reduced_resultant(GEN x, GEN y, GEN p, GEN pm)
     822             : {
     823      985828 :   pari_sp av = avma;
     824             :   GEN z;
     825      985828 :   if (lgefint(pm) == 3)
     826             :   {
     827      974445 :     ulong q = pm[2];
     828      974445 :     z = Zlx_sylvester_echelon(ZX_to_Flx(x,q), ZX_to_Flx(y,q),0,p[2],q);
     829      974511 :     if (lg(z) > 1)
     830             :     {
     831      974511 :       ulong c = ucoeff(z,1,1);
     832      974511 :       if (c) return gc_utoipos(av, c);
     833             :     }
     834             :   }
     835             :   else
     836             :   {
     837       11383 :     z = ZpX_sylvester_echelon(x,y,0,p,pm);
     838       11383 :     if (lg(z) > 1)
     839             :     {
     840       11383 :       GEN c = gcoeff(z,1,1);
     841       11383 :       if (signe(c)) return gerepileuptoint(av, c);
     842             :     }
     843             :   }
     844      127627 :   set_avma(av); return gen_0;
     845             : }
     846             : /* Assume Res(f,g) divides p^M. Return Res(f, g), using dynamic p-adic
     847             :  * precision (until result is nonzero or p^M). */
     848             : GEN
     849      926054 : ZpX_reduced_resultant_fast(GEN f, GEN g, GEN p, long M)
     850             : {
     851      926054 :   GEN R, q = NULL;
     852             :   long m;
     853      926054 :   m = init_m(p); if (m < 1) m = 1;
     854       59771 :   for(;; m <<= 1) {
     855      985822 :     if (M < 2*m) break;
     856       91704 :     q = q? sqri(q): powiu(p, m); /* p^m */
     857       91704 :     R = ZpX_reduced_resultant(f,g, p, q); if (signe(R)) return R;
     858             :   }
     859      894118 :   q = powiu(p, M);
     860      894125 :   R = ZpX_reduced_resultant(f,g, p, q); return signe(R)? R: q;
     861             : }
     862             : 
     863             : /* v_p(Res(x,y) mod p^m), assumes (lc(x),p) = 1 */
     864             : static long
     865     3484279 : ZpX_resultant_val_i(GEN x, GEN y, GEN p, GEN pm)
     866             : {
     867     3484279 :   pari_sp av = avma;
     868             :   GEN z;
     869             :   long i, l, v;
     870     3484279 :   if (lgefint(pm) == 3)
     871             :   {
     872     3430916 :     ulong q = pm[2], pp = p[2];
     873     3430916 :     z = Zlx_sylvester_echelon(ZX_to_Flx(x,q), ZX_to_Flx(y,q), 1, pp, q);
     874     3431027 :     if (!z) return gc_long(av,-1); /* failure */
     875     3246240 :     v = 0; l = lg(z);
     876    13518055 :     for (i = 1; i < l; i++) v += u_lval(ucoeff(z,i,i), pp);
     877             :   }
     878             :   else
     879             :   {
     880       53363 :     z = ZpX_sylvester_echelon(x, y, 1, p, pm);
     881       53363 :     if (!z) return gc_long(av,-1); /* failure */
     882       52675 :     v = 0; l = lg(z);
     883      194330 :     for (i = 1; i < l; i++) v += Z_pval(gcoeff(z,i,i), p);
     884             :   }
     885     3298909 :   return v;
     886             : }
     887             : 
     888             : /* assume (lc(f),p) = 1; no assumption on g */
     889             : long
     890     3440130 : ZpX_resultant_val(GEN f, GEN g, GEN p, long M)
     891             : {
     892     3440130 :   pari_sp av = avma;
     893     3440130 :   GEN q = NULL;
     894             :   long v, m;
     895     3440130 :   m = init_m(p); if (m < 2) m = 2;
     896       44118 :   for(;; m <<= 1) {
     897     3484241 :     if (m > M) m = M;
     898     3484241 :     q = q? sqri(q): powiu(p, m); /* p^m */
     899     3484280 :     v = ZpX_resultant_val_i(f,g, p, q); if (v >= 0) return gc_long(av,v);
     900      185474 :     if (m == M) return gc_long(av,M);
     901             :   }
     902             : }
     903             : 
     904             : /* assume f separable and (lc(f),p) = 1 */
     905             : long
     906      182384 : ZpX_disc_val(GEN f, GEN p)
     907             : {
     908      182384 :   pari_sp av = avma;
     909             :   long v;
     910      182384 :   if (degpol(f) == 1) return 0;
     911      182384 :   v = ZpX_resultant_val(f, ZX_deriv(f), p, LONG_MAX);
     912      182387 :   return gc_long(av,v);
     913             : }
     914             : 
     915             : /* *e a ZX, *d, *z in Z, *d = p^(*vd). Simplify e / d by cancelling a
     916             :  * common factor p^v; if z!=NULL, update it by cancelling the same power of p */
     917             : static void
     918     3530130 : update_den(GEN p, GEN *e, GEN *d, long *vd, GEN *z)
     919             : {
     920             :   GEN newe;
     921     3530130 :   long ve = ZX_pvalrem(*e, p, &newe);
     922     3530273 :   if (ve) {
     923             :     GEN newd;
     924     1735858 :     long v = minss(*vd, ve);
     925     1735853 :     if (v) {
     926     1735943 :       if (v == *vd)
     927             :       { /* rare, denominator cancelled */
     928      378915 :         if (ve != v) newe = ZX_Z_mul(newe, powiu(p, ve - v));
     929      378915 :         newd = gen_1;
     930      378915 :         *vd = 0;
     931      378915 :         if (z) *z =diviiexact(*z, powiu(p, v));
     932             :       }
     933             :       else
     934             :       { /* v = ve < vd, generic case */
     935     1357028 :         GEN q = powiu(p, v);
     936     1357070 :         newd = diviiexact(*d, q);
     937     1356894 :         *vd -= v;
     938     1356894 :         if (z) *z = diviiexact(*z, q);
     939             :       }
     940     1735772 :       *e = newe;
     941     1735772 :       *d = newd;
     942             :     }
     943             :   }
     944     3530097 : }
     945             : 
     946             : /* return denominator, a power of p */
     947             : static GEN
     948     2724649 : QpX_denom(GEN x)
     949             : {
     950     2724649 :   long i, l = lg(x);
     951     2724649 :   GEN maxd = gen_1;
     952     9401037 :   for (i=2; i<l; i++)
     953             :   {
     954     6676391 :     GEN d = gel(x,i);
     955     6676391 :     if (typ(d) == t_FRAC && cmpii(gel(d,2), maxd) > 0) maxd = gel(d,2);
     956             :   }
     957     2724646 :   return maxd;
     958             : }
     959             : static GEN
     960      505944 : QpXV_denom(GEN x)
     961             : {
     962      505944 :   long l = lg(x), i;
     963      505944 :   GEN maxd = gen_1;
     964     1502099 :   for (i = 1; i < l; i++)
     965             :   {
     966      996156 :     GEN d = QpX_denom(gel(x,i));
     967      996155 :     if (cmpii(d, maxd) > 0) maxd = d;
     968             :   }
     969      505943 :   return maxd;
     970             : }
     971             : 
     972             : static GEN
     973     1728514 : QpX_remove_denom(GEN x, GEN p, GEN *pdx, long *pv)
     974             : {
     975     1728514 :   *pdx = QpX_denom(x);
     976     1728492 :   if (*pdx == gen_1) { *pv = 0; *pdx = NULL; }
     977             :   else {
     978     1255588 :     x = Q_muli_to_int(x,*pdx);
     979     1255512 :     *pv = Z_pval(*pdx, p);
     980             :   }
     981     1728425 :   return x;
     982             : }
     983             : 
     984             : /* p^v * f o g mod (T,q). q = p^vq  */
     985             : static GEN
     986      284801 : compmod(GEN p, GEN f, GEN g, GEN T, GEN q, long v)
     987             : {
     988      284801 :   GEN D = NULL, z, df, dg, qD;
     989      284801 :   long vD = 0, vdf, vdg;
     990             : 
     991      284801 :   f = QpX_remove_denom(f, p, &df, &vdf);
     992      284798 :   if (typ(g) == t_VEC) /* [num,den,v_p(den)] */
     993           0 :   { vdg = itos(gel(g,3)); dg = gel(g,2); g = gel(g,1); }
     994             :   else
     995      284798 :     g = QpX_remove_denom(g, p, &dg, &vdg);
     996      284801 :   if (df) { D = df; vD = vdf; }
     997      284801 :   if (dg) {
     998       55344 :     long degf = degpol(f);
     999       55344 :     D = mul_content(D, powiu(dg, degf));
    1000       55344 :     vD += degf * vdg;
    1001             :   }
    1002      284801 :   qD = D ? mulii(q, D): q;
    1003      284792 :   if (dg) f = FpX_rescale(f, dg, qD);
    1004      284792 :   z = FpX_FpXQ_eval(f, g, T, qD);
    1005      284795 :   if (!D) {
    1006           0 :     if (v) {
    1007           0 :       if (v > 0)
    1008           0 :         z = ZX_Z_mul(z, powiu(p, v));
    1009             :       else
    1010           0 :         z = RgX_Rg_div(z, powiu(p, -v));
    1011             :     }
    1012           0 :     return z;
    1013             :   }
    1014      284795 :   update_den(p, &z, &D, &vD, NULL);
    1015      284801 :   qD = mulii(D,q);
    1016      284791 :   if (v) vD -= v;
    1017      284791 :   z = FpX_center_i(z, qD, shifti(qD,-1));
    1018      284798 :   if (vD > 0)
    1019      284798 :     z = RgX_Rg_div(z, powiu(p, vD));
    1020           0 :   else if (vD < 0)
    1021           0 :     z = ZX_Z_mul(z, powiu(p, -vD));
    1022      284803 :   return z;
    1023             : }
    1024             : 
    1025             : /* fast implementation of ZM_hnfmodid(M, D) / D, D = p^k */
    1026             : static GEN
    1027      451325 : ZpM_hnfmodid(GEN M, GEN p, GEN D)
    1028             : {
    1029      451325 :   long i, l = lg(M);
    1030      451325 :   M = RgM_Rg_div(ZpM_echelon(M,0,p,D), D);
    1031     2011970 :   for (i = 1; i < l; i++)
    1032     1560640 :     if (gequal0(gcoeff(M,i,i))) gcoeff(M,i,i) = gen_1;
    1033      451330 :   return M;
    1034             : }
    1035             : 
    1036             : /* Return Z-basis for Z[a] + U(a)/p Z[a] in Z[t]/(f), mf = v_p(disc f), U
    1037             :  * a ZX. Special cases: a = t is coded as NULL, U = 0 is coded as NULL */
    1038             : static GEN
    1039      617184 : dbasis(GEN p, GEN f, long mf, GEN a, GEN U)
    1040             : {
    1041      617184 :   long n = degpol(f), i, dU;
    1042             :   GEN b, h;
    1043             : 
    1044      617184 :   if (n == 1) return matid(1);
    1045      617184 :   if (a && gequalX(a)) a = NULL;
    1046      617184 :   if (DEBUGLEVEL>5)
    1047             :   {
    1048           0 :     err_printf("  entering Dedekind Basis with parameters p=%Ps\n",p);
    1049           0 :     err_printf("  f = %Ps,\n  a = %Ps\n",f, a? a: pol_x(varn(f)));
    1050             :   }
    1051      617191 :   if (a)
    1052             :   {
    1053      198358 :     GEN pd = powiu(p, mf >> 1);
    1054      198358 :     GEN da, pdp = mulii(pd,p), D = pdp;
    1055             :     long vda;
    1056      198355 :     dU = U ? degpol(U): 0;
    1057      198355 :     b = cgetg(n+1, t_MAT);
    1058      198355 :     h = scalarpol(pd, varn(f));
    1059      198357 :     a = QpX_remove_denom(a, p, &da, &vda);
    1060      198354 :     if (da) D = mulii(D, da);
    1061      198352 :     gel(b,1) = scalarcol_shallow(pd, n);
    1062      564466 :     for (i=2; i<=n; i++)
    1063             :     {
    1064      366115 :       if (i == dU+1)
    1065           0 :         h = compmod(p, U, mkvec3(a,da,stoi(vda)), f, pdp, (mf>>1) - 1);
    1066             :       else
    1067             :       {
    1068      366115 :         h = FpXQ_mul(h, a, f, D);
    1069      366112 :         if (da) h = ZX_Z_divexact(h, da);
    1070             :       }
    1071      366093 :       gel(b,i) = RgX_to_RgC(h,n);
    1072             :     }
    1073      198351 :     return ZpM_hnfmodid(b, p, pd);
    1074             :   }
    1075             :   else
    1076             :   {
    1077      418833 :     if (!U) return matid(n);
    1078      418833 :     dU = degpol(U);
    1079      418831 :     if (dU == n) return matid(n);
    1080      418831 :     U = FpX_normalize(U, p);
    1081      418840 :     b = cgetg(n+1, t_MAT);
    1082     1618946 :     for (i = 1; i <= dU; i++) gel(b,i) = vec_ei(n, i);
    1083      418843 :     h = RgX_Rg_div(U, p);
    1084      470243 :     for ( ; i <= n; i++)
    1085             :     {
    1086      470243 :       gel(b, i) = RgX_to_RgC(h,n);
    1087      470243 :       if (i == n) break;
    1088       51394 :       h = RgX_shift_shallow(h,1);
    1089             :     }
    1090      418849 :     return b;
    1091             :   }
    1092             : }
    1093             : 
    1094             : static GEN
    1095      505945 : get_partial_order_as_pols(GEN p, GEN f)
    1096             : {
    1097      505945 :   GEN O = maxord(p, f, -1);
    1098      505935 :   long v = varn(f);
    1099      505935 :   return O == gen_1? pol_x_powers(degpol(f), v): RgM_to_RgXV(O, v);
    1100             : }
    1101             : 
    1102             : static long
    1103        2227 : p_is_prime(decomp_t *S)
    1104             : {
    1105        2227 :   if (S->pisprime < 0) S->pisprime = BPSW_psp(S->p);
    1106        2227 :   return S->pisprime;
    1107             : }
    1108             : static GEN ZpX_monic_factor_squarefree(GEN f, GEN p, long prec);
    1109             : 
    1110             : /* if flag = 0, maximal order, else factorization to precision r = flag */
    1111             : static GEN
    1112      253260 : Decomp(decomp_t *S, long flag)
    1113             : {
    1114      253260 :   pari_sp av = avma;
    1115             :   GEN fred, pr2, pr, pk, ph2, ph, b1, b2, a, e, de, f1, f2, dt, th, chip;
    1116      253260 :   GEN p = S->p;
    1117      253260 :   long vde, vdt, k, r = maxss(flag, 2*S->df + 1);
    1118             : 
    1119      253260 :   if (DEBUGLEVEL>5) err_printf("  entering Decomp: %Ps^%ld\n  f = %Ps\n",
    1120             :                                p, r, S->f);
    1121      253260 :   else if (DEBUGLEVEL>2) err_printf("  entering Decomp\n");
    1122      253260 :   chip = FpX_red(S->chi, p);
    1123      253255 :   if (!FpX_valrem(chip, S->nu, p, &b1))
    1124             :   {
    1125           0 :     if (!p_is_prime(S)) pari_err_PRIME("Decomp",p);
    1126           0 :     pari_err_BUG("Decomp (not a factor)");
    1127             :   }
    1128      253262 :   b2 = FpX_div(chip, b1, p);
    1129      253261 :   a = FpX_mul(FpXQ_inv(b2, b1, p), b2, p);
    1130             :   /* E = e / de, e in Z[X], de in Z,  E = a(phi) mod (f, p) */
    1131      253257 :   th = QpX_remove_denom(S->phi, p, &dt, &vdt);
    1132      253262 :   if (dt)
    1133             :   {
    1134      122230 :     long dega = degpol(a);
    1135      122230 :     vde = dega * vdt;
    1136      122230 :     de = powiu(dt, dega);
    1137      122229 :     pr = mulii(p, de);
    1138      122228 :     a = FpX_rescale(a, dt, pr);
    1139             :   }
    1140             :   else
    1141             :   {
    1142      131032 :     vde = 0;
    1143      131032 :     de = gen_1;
    1144      131032 :     pr = p;
    1145             :   }
    1146      253261 :   e = FpX_FpXQ_eval(a, th, S->f, pr);
    1147      253256 :   update_den(p, &e, &de, &vde, NULL);
    1148             : 
    1149      253261 :   pk = p; k = 1;
    1150             :   /* E, (1 - E) tend to orthogonal idempotents in Zp[X]/(f) */
    1151     1171190 :   while (k < r + vde)
    1152             :   { /* E <-- E^2(3-2E) mod p^2k, with E = e/de */
    1153             :     GEN D;
    1154      917933 :     pk = sqri(pk); k <<= 1;
    1155      917915 :     e = ZX_mul(ZX_sqr(e), Z_ZX_sub(mului(3,de), gmul2n(e,1)));
    1156      917955 :     de= mulii(de, sqri(de));
    1157      917915 :     vde *= 3;
    1158      917915 :     D = mulii(pk, de);
    1159      917917 :     e = FpX_rem(e, centermod(S->f, D), D); /* e/de defined mod pk */
    1160      917903 :     update_den(p, &e, &de, &vde, NULL);
    1161             :   }
    1162             :   /* required precision of the factors */
    1163      253257 :   pr = powiu(p, r); pr2 = shifti(pr, -1);
    1164      253258 :   ph = mulii(de,pr);ph2 = shifti(ph, -1);
    1165      253259 :   e = FpX_center_i(FpX_red(e, ph), ph, ph2);
    1166      253258 :   fred = FpX_red(S->f, ph);
    1167             : 
    1168      253257 :   f1 = ZpX_gcd(fred, Z_ZX_sub(de, e), p, ph); /* p-adic gcd(f, 1-e) */
    1169      253261 :   if (!is_pm1(de))
    1170             :   {
    1171      122229 :     fred = FpX_red(fred, pr);
    1172      122230 :     f1 = FpX_red(f1, pr);
    1173             :   }
    1174      253259 :   f2 = FpX_div(fred,f1, pr);
    1175      253259 :   f1 = FpX_center_i(f1, pr, pr2);
    1176      253258 :   f2 = FpX_center_i(f2, pr, pr2);
    1177             : 
    1178      253255 :   if (DEBUGLEVEL>5)
    1179           0 :     err_printf("  leaving Decomp: f1 = %Ps\nf2 = %Ps\ne = %Ps\nde= %Ps\n", f1,f2,e,de);
    1180             : 
    1181      253255 :   if (flag < 0)
    1182             :   {
    1183           0 :     GEN m = vconcat(ZpX_primedec(f1, p), ZpX_primedec(f2, p));
    1184           0 :     return sort_factor(m, (void*)&cmpii, &cmp_nodata);
    1185             :   }
    1186      253255 :   else if (flag)
    1187             :   {
    1188         287 :     gerepileall(av, 2, &f1, &f2);
    1189         287 :     return shallowconcat(ZpX_monic_factor_squarefree(f1, p, flag),
    1190             :                          ZpX_monic_factor_squarefree(f2, p, flag));
    1191             :   } else {
    1192             :     GEN D, d1, d2, B1, B2, M;
    1193             :     long n, n1, n2, i;
    1194      252968 :     gerepileall(av, 4, &f1, &f2, &e, &de);
    1195      252974 :     D = de;
    1196      252974 :     B1 = get_partial_order_as_pols(p,f1); n1 = lg(B1)-1;
    1197      252974 :     B2 = get_partial_order_as_pols(p,f2); n2 = lg(B2)-1; n = n1+n2;
    1198      252972 :     d1 = QpXV_denom(B1);
    1199      252972 :     d2 = QpXV_denom(B2); if (cmpii(d1, d2) < 0) d1 = d2;
    1200      252973 :     if (d1 != gen_1) {
    1201      155606 :       B1 = Q_muli_to_int(B1, d1);
    1202      155606 :       B2 = Q_muli_to_int(B2, d1);
    1203      155605 :       D = mulii(d1, D);
    1204             :     }
    1205      252972 :     fred = centermod_i(S->f, D, shifti(D,-1));
    1206      252973 :     M = cgetg(n+1, t_MAT);
    1207      799203 :     for (i=1; i<=n1; i++)
    1208      546230 :       gel(M,i) = RgX_to_RgC(FpX_rem(FpX_mul(gel(B1,i),e,D), fred, D), n);
    1209      252973 :     e = Z_ZX_sub(de, e); B2 -= n1;
    1210      702899 :     for (   ; i<=n; i++)
    1211      449926 :       gel(M,i) = RgX_to_RgC(FpX_rem(FpX_mul(gel(B2,i),e,D), fred, D), n);
    1212      252973 :     return ZpM_hnfmodid(M, p, D);
    1213             :   }
    1214             : }
    1215             : 
    1216             : /* minimum extension valuation: L/E */
    1217             : static void
    1218      618596 : vstar(GEN p,GEN h, long *L, long *E)
    1219             : {
    1220      618596 :   long first, j, k, v, w, m = degpol(h);
    1221             : 
    1222      618594 :   first = 1; k = 1; v = 0;
    1223     2549760 :   for (j=1; j<=m; j++)
    1224             :   {
    1225     1931165 :     GEN c = gel(h, m-j+2);
    1226     1931165 :     if (signe(c))
    1227             :     {
    1228     1856139 :       w = Z_pval(c,p);
    1229     1856140 :       if (first || w*k < v*j) { v = w; k = j; }
    1230     1856140 :       first = 0;
    1231             :     }
    1232             :   }
    1233             :   /* v/k = min_j ( v_p(h_{m-j}) / j ) */
    1234      618595 :   w = (long)ugcd(v,k);
    1235      618601 :   *L = v/w;
    1236      618601 :   *E = k/w;
    1237      618601 : }
    1238             : 
    1239             : static GEN
    1240       63742 : redelt_i(GEN a, GEN N, GEN p, GEN *pda, long *pvda)
    1241             : {
    1242             :   GEN z;
    1243       63742 :   a = Q_remove_denom(a, pda);
    1244       63742 :   *pvda = 0;
    1245       63742 :   if (*pda)
    1246             :   {
    1247       63742 :     long v = Z_pvalrem(*pda, p, &z);
    1248       63741 :     if (v) {
    1249       63741 :       *pda = powiu(p, v);
    1250       63741 :       *pvda = v;
    1251       63741 :       N  = mulii(*pda, N);
    1252             :     }
    1253             :     else
    1254           0 :       *pda = NULL;
    1255       63741 :     if (!is_pm1(z)) a = ZX_Z_mul(a, Fp_inv(z, N));
    1256             :   }
    1257       63741 :   return centermod(a, N);
    1258             : }
    1259             : /* reduce the element a modulo N [ a power of p ], taking first care of the
    1260             :  * denominators */
    1261             : static GEN
    1262       48018 : redelt(GEN a, GEN N, GEN p)
    1263             : {
    1264             :   GEN da;
    1265             :   long vda;
    1266       48018 :   a = redelt_i(a, N, p, &da, &vda);
    1267       48018 :   if (da) a = RgX_Rg_div(a, da);
    1268       48018 :   return a;
    1269             : }
    1270             : 
    1271             : /* compute the c first Newton sums modulo pp of the
    1272             :    characteristic polynomial of a/d mod chi, d > 0 power of p (NULL = gen_1),
    1273             :    a, chi in Zp[X], vda = v_p(da)
    1274             :    ns = Newton sums of chi */
    1275             : static GEN
    1276      699226 : newtonsums(GEN p, GEN a, GEN da, long vda, GEN chi, long c, GEN pp, GEN ns)
    1277             : {
    1278             :   GEN va, pa, dpa, s;
    1279      699226 :   long j, k, vdpa, lns = lg(ns);
    1280             :   pari_sp av;
    1281             : 
    1282      699226 :   a = centermod(a, pp); av = avma;
    1283      699188 :   dpa = pa = NULL; /* -Wall */
    1284      699188 :   vdpa = 0;
    1285      699188 :   va = zerovec(c);
    1286     2878830 :   for (j = 1; j <= c; j++)
    1287             :   { /* pa/dpa = (a/d)^(j-1) mod (chi, pp), dpa = p^vdpa */
    1288             :     long l;
    1289     2186121 :     pa = j == 1? a: FpXQ_mul(pa, a, chi, pp);
    1290     2186276 :     l = lg(pa); if (l == 2) break;
    1291     2186276 :     if (lns < l) l = lns;
    1292             : 
    1293     2186276 :     if (da) {
    1294     2052046 :       dpa = j == 1? da: mulii(dpa, da);
    1295     2051858 :       vdpa += vda;
    1296     2051858 :       update_den(p, &pa, &dpa, &vdpa, &pp);
    1297             :     }
    1298     2186016 :     s = mulii(gel(pa,2), gel(ns,2)); /* k = 2 */
    1299    10749638 :     for (k = 3; k < l; k++) s = addii(s, mulii(gel(pa,k), gel(ns,k)));
    1300     2185857 :     if (da) {
    1301             :       GEN r;
    1302     2051658 :       s = dvmdii(s, dpa, &r);
    1303     2051569 :       if (r != gen_0) return NULL;
    1304             :     }
    1305     2179326 :     gel(va,j) = centermodii(s, pp, shifti(pp,-1));
    1306             : 
    1307     2179333 :     if (gc_needed(av, 1))
    1308             :     {
    1309           7 :       if(DEBUGMEM>1) pari_warn(warnmem, "newtonsums");
    1310           7 :       gerepileall(av, dpa?4:3, &pa, &va, &pp, &dpa);
    1311             :     }
    1312             :   }
    1313      692709 :   for (; j <= c; j++) gel(va,j) = gen_0;
    1314      692709 :   return va;
    1315             : }
    1316             : 
    1317             : /* compute the characteristic polynomial of a/da mod chi (a in Z[X]), given
    1318             :  * by its Newton sums to a precision of pp using Newton sums */
    1319             : static GEN
    1320      692710 : newtoncharpoly(GEN pp, GEN p, GEN NS)
    1321             : {
    1322      692710 :   long n = lg(NS)-1, j, k;
    1323      692710 :   GEN c = cgetg(n + 2, t_VEC), pp2 = shifti(pp,-1);
    1324             : 
    1325      692753 :   gel(c,1) = (n & 1 ? gen_m1: gen_1);
    1326     2862626 :   for (k = 2; k <= n+1; k++)
    1327             :   {
    1328     2169956 :     pari_sp av2 = avma;
    1329     2169956 :     GEN s = gen_0;
    1330             :     ulong z;
    1331     2169956 :     long v = u_pvalrem(k - 1, p, &z);
    1332     9134260 :     for (j = 1; j < k; j++)
    1333             :     {
    1334     6965016 :       GEN t = mulii(gel(NS,j), gel(c,k-j));
    1335     6964075 :       if (!odd(j)) t = negi(t);
    1336     6964278 :       s = addii(s, t);
    1337             :     }
    1338     2169244 :     if (v) {
    1339      829392 :       s = gdiv(s, powiu(p, v));
    1340      829414 :       if (typ(s) != t_INT) return NULL;
    1341             :     }
    1342     2169182 :     s = mulii(s, Fp_inv(utoipos(z), pp));
    1343     2169532 :     gel(c,k) = gerepileuptoint(av2, Fp_center_i(s, pp, pp2));
    1344             :   }
    1345     1841510 :   for (k = odd(n)? 1: 2; k <= n+1; k += 2) gel(c,k) = negi(gel(c,k));
    1346      692681 :   return gtopoly(c, 0);
    1347             : }
    1348             : 
    1349             : static void
    1350      699156 : manage_cache(decomp_t *S, GEN f, GEN pp)
    1351             : {
    1352      699156 :   GEN t = S->precns;
    1353             : 
    1354      699156 :   if (!t) t = mulii(S->pmf, powiu(S->p, S->df));
    1355      699153 :   if (cmpii(t, pp) < 0) t = pp;
    1356             : 
    1357      699154 :   if (!S->precns || !RgX_equal(f, S->nsf) || cmpii(S->precns, t) < 0)
    1358             :   {
    1359      516679 :     if (DEBUGLEVEL>4)
    1360           0 :       err_printf("  Precision for cached Newton sums for %Ps: %Ps -> %Ps\n",
    1361           0 :                  f, S->precns? S->precns: gen_0, t);
    1362      516679 :     S->nsf = f;
    1363      516679 :     S->ns = FpX_Newton(f, degpol(f), t);
    1364      516705 :     S->precns = t;
    1365             :   }
    1366      699198 : }
    1367             : 
    1368             : /* return NULL if a mod f is not an integer
    1369             :  * The denominator of any integer in Zp[X]/(f) divides pdr */
    1370             : static GEN
    1371      699220 : mycaract(decomp_t *S, GEN f, GEN a, GEN pp, GEN pdr)
    1372             : {
    1373             :   pari_sp av;
    1374             :   GEN d, chi, prec1, prec2, prec3, ns;
    1375      699220 :   long vd, n = degpol(f);
    1376             : 
    1377      699219 :   if (gequal0(a)) return pol_0(varn(f));
    1378             : 
    1379      699219 :   a = QpX_remove_denom(a, S->p, &d, &vd);
    1380      699189 :   prec1 = pp;
    1381      699189 :   if (lgefint(S->p) == 3)
    1382      699137 :     prec1 = mulii(prec1, powiu(S->p, factorial_lval(n, itou(S->p))));
    1383      699161 :   if (d)
    1384             :   {
    1385      634761 :     GEN p1 = powiu(d, n);
    1386      634796 :     prec2 = mulii(prec1, p1);
    1387      634766 :     prec3 = mulii(prec1, gmin_shallow(mulii(p1, d), pdr));
    1388             :   }
    1389             :   else
    1390       64400 :     prec2 = prec3 = prec1;
    1391      699164 :   manage_cache(S, f, prec3);
    1392             : 
    1393      699226 :   av = avma;
    1394      699226 :   ns = newtonsums(S->p, a, d, vd, f, n, prec2, S->ns);
    1395      699150 :   if (!ns) return NULL;
    1396      692708 :   chi = newtoncharpoly(prec1, S->p, ns);
    1397      692780 :   if (!chi) return NULL;
    1398      692696 :   setvarn(chi, varn(f));
    1399      692696 :   return gerepileupto(av, centermod(chi, pp));
    1400             : }
    1401             : 
    1402             : static GEN
    1403      635433 : get_nu(GEN chi, GEN p, long *ptl)
    1404             : { /* split off powers of x first for efficiency */
    1405      635433 :   long v = ZX_valrem(FpX_red(chi,p), &chi), n;
    1406             :   GEN P;
    1407      635412 :   if (!degpol(chi)) { *ptl = 1; return pol_x(varn(chi)); }
    1408      470296 :   P = gel(FpX_factor(chi,p), 1); n = lg(P)-1;
    1409      470321 :   *ptl = v? n+1: n; return gel(P,n);
    1410             : }
    1411             : 
    1412             : /* Factor characteristic polynomial chi of phi mod p. If it splits, update
    1413             :  * S->{phi, chi, nu} and return 1. In any case, set *nu to an irreducible
    1414             :  * factor mod p of chi */
    1415             : static int
    1416      475494 : split_char(decomp_t *S, GEN chi, GEN phi, GEN phi0, GEN *nu)
    1417             : {
    1418             :   long l;
    1419      475494 :   *nu  = get_nu(chi, S->p, &l);
    1420      475500 :   if (l == 1) return 0; /* single irreducible factor: doesn't split */
    1421             :   /* phi o phi0 mod (p, f) */
    1422      122228 :   S->phi = compmod(S->p, phi, phi0, S->f, S->p, 0);
    1423      122229 :   S->chi = chi;
    1424      122229 :   S->nu = *nu; return 1;
    1425             : }
    1426             : 
    1427             : /* Return the prime element in Zp[phi], a t_INT (iff *Ep = 1) or QX;
    1428             :  * nup, chip are ZX. phi = NULL codes X
    1429             :  * If *Ep < oE or Ep divides Ediv (!=0) return NULL (uninteresting) */
    1430             : static GEN
    1431      558327 : getprime(decomp_t *S, GEN phi, GEN chip, GEN nup, long *Lp, long *Ep,
    1432             :          long oE, long Ediv)
    1433             : {
    1434             :   GEN z, chin, q, qp;
    1435             :   long r, s;
    1436             : 
    1437      558327 :   if (phi && dvdii(constant_coeff(chip), S->psc))
    1438             :   {
    1439        1642 :     chip = mycaract(S, S->chi, phi, S->pmf, S->prc);
    1440        1642 :     if (dvdii(constant_coeff(chip), S->pmf))
    1441        1253 :       chip = ZXQ_charpoly(phi, S->chi, varn(chip));
    1442             :   }
    1443      558325 :   if (degpol(nup) == 1)
    1444             :   {
    1445      519682 :     GEN c = gel(nup,2); /* nup = X + c */
    1446      519682 :     chin = signe(c)? RgX_translate(chip, negi(c)): chip;
    1447             :   }
    1448             :   else
    1449       38640 :     chin = ZXQ_charpoly(nup, chip, varn(chip));
    1450             : 
    1451      558334 :   vstar(S->p, chin, Lp, Ep);
    1452      558341 :   if (*Ep < oE || (Ediv && Ediv % *Ep == 0)) return NULL;
    1453             : 
    1454      439526 :   if (*Ep == 1) return S->p;
    1455      303287 :   (void)cbezout(*Lp, -*Ep, &r, &s); /* = 1 */
    1456      303292 :   if (r <= 0)
    1457             :   {
    1458       60025 :     long t = 1 + ((-r) / *Ep);
    1459       60025 :     r += t * *Ep;
    1460       60025 :     s += t * *Lp;
    1461             :   }
    1462             :   /* r > 0 minimal such that r L/E - s = 1/E
    1463             :    * pi = nu^r / p^s is an element of valuation 1/E,
    1464             :    * so is pi + O(p) since 1/E < 1. May compute nu^r mod p^(s+1) */
    1465      303292 :   q = powiu(S->p, s); qp = mulii(q, S->p);
    1466      303281 :   nup = FpXQ_powu(nup, r, S->chi, qp);
    1467      303290 :   if (!phi) return RgX_Rg_div(nup, q); /* phi = X : no composition */
    1468       48018 :   z = compmod(S->p, nup, phi, S->chi, qp, -s);
    1469       48018 :   return signe(z)? z: NULL;
    1470             : }
    1471             : 
    1472             : static int
    1473      274500 : update_phi(decomp_t *S)
    1474             : {
    1475      274500 :   GEN PHI = NULL, prc, psc, X = pol_x(varn(S->f));
    1476             :   long k, vpsc;
    1477      274500 :   for (k = 1;; k++)
    1478             :   {
    1479      276795 :     prc = ZpX_reduced_resultant_fast(S->chi, ZX_deriv(S->chi), S->p, S->vpsc);
    1480             :     /* if prc == S->psc then either chi is not separable or
    1481             :        the reduced discriminant of chi is too large */
    1482      276793 :     if (!equalii(prc, S->psc)) break;
    1483             : 
    1484             :     /* increase precision */
    1485        2295 :     S->vpsc = maxss(S->vpsf, S->vpsc + 1);
    1486        2295 :     S->psc = (S->vpsc == S->vpsf)? S->psf: mulii(S->psc, S->p);
    1487             : 
    1488        2295 :     PHI = S->phi;
    1489        2295 :     if (S->phi0) PHI = compmod(S->p, PHI, S->phi0, S->f, S->psc, 0);
    1490             :     /* change phi (in case not separable) */
    1491        2295 :     PHI = gadd(PHI, ZX_Z_mul(X, mului(k, S->p)));
    1492        2295 :     S->chi = mycaract(S, S->f, PHI, S->psc, S->pdf);
    1493             :   }
    1494      274497 :   psc = mulii(sqri(prc), S->p);
    1495      274493 :   vpsc = 2*Z_pval(prc, S->p) + 1;
    1496             : 
    1497      274493 :   if (!PHI) /* break out of above loop immediately (k = 1) */
    1498             :   {
    1499      272197 :     PHI = S->phi;
    1500      272197 :     if (S->phi0) PHI = compmod(S->p, PHI, S->phi0, S->f, psc, 0);
    1501      272198 :     if (S->phi0 || cmpii(psc,S->psc) > 0)
    1502             :     {
    1503             :       for(;;)
    1504             :       {
    1505      112837 :         S->chi = mycaract(S, S->f, PHI, psc, S->pdf);
    1506      112842 :         prc = ZpX_reduced_resultant_fast(S->chi, ZX_deriv(S->chi), S->p, vpsc);
    1507      112842 :         if (!equalii(prc, psc)) break;
    1508         497 :         psc = mulii(psc, S->p); vpsc++;
    1509             :         /* it can happen that S->chi is never squarefree: then change PHI */
    1510         497 :         if (vpsc > 2*S->mf) PHI = gadd(PHI, ZX_Z_mul(X, S->p));
    1511             :       }
    1512      112345 :       psc = mulii(sqri(prc), S->p);
    1513      112342 :       vpsc = 2*Z_pval(prc, S->p) + 1;
    1514             :     }
    1515             :   }
    1516      274496 :   S->phi = PHI;
    1517      274496 :   S->chi = FpX_red(S->chi, psc);
    1518             : 
    1519             :   /* may happen if p is unramified */
    1520      274497 :   if (is_pm1(prc)) return 0;
    1521      230391 :   S->prc = prc;
    1522      230391 :   S->psc = psc;
    1523      230391 :   S->vpsc = vpsc; return 1;
    1524             : }
    1525             : 
    1526             : /* return 1 if at least 2 factors mod p ==> chi splits
    1527             :  * Replace S->phi such that F increases (to D) */
    1528             : static int
    1529       66539 : testb2(decomp_t *S, long D, GEN theta)
    1530             : {
    1531       66539 :   long v = varn(S->chi), dlim = degpol(S->chi)-1;
    1532       66539 :   GEN T0 = S->phi, chi, phi, nu;
    1533       66539 :   if (DEBUGLEVEL>4) err_printf("  Increasing Fa\n");
    1534             :   for (;;)
    1535             :   {
    1536       66599 :     phi = gadd(theta, random_FpX(dlim, v, S->p));
    1537       66599 :     chi = mycaract(S, S->chi, phi, S->psf, S->prc);
    1538             :     /* phi nonprimary ? */
    1539       66598 :     if (split_char(S, chi, phi, T0, &nu)) return 1;
    1540       66599 :     if (degpol(nu) == D) break;
    1541             :   }
    1542             :   /* F_phi=lcm(F_alpha, F_theta)=D and E_phi=E_alpha */
    1543       66539 :   S->phi0 = T0;
    1544       66539 :   S->chi = chi;
    1545       66539 :   S->phi = phi;
    1546       66539 :   S->nu = nu; return 0;
    1547             : }
    1548             : 
    1549             : /* return 1 if at least 2 factors mod p ==> chi can be split.
    1550             :  * compute a new S->phi such that E = lcm(Ea, Et);
    1551             :  * A a ZX, T a t_INT (iff Et = 1, probably impossible ?) or QX */
    1552             : static int
    1553       48018 : testc2(decomp_t *S, GEN A, long Ea, GEN T, long Et)
    1554             : {
    1555       48018 :   GEN c, chi, phi, nu, T0 = S->phi;
    1556             : 
    1557       48018 :   if (DEBUGLEVEL>4) err_printf("  Increasing Ea\n");
    1558       48018 :   if (Et == 1) /* same as other branch, split for efficiency */
    1559           0 :     c = A; /* Et = 1 => s = 1, r = 0, t = 0 */
    1560             :   else
    1561             :   {
    1562             :     long r, s, t;
    1563       48018 :     (void)cbezout(Ea, Et, &r, &s); t = 0;
    1564       48115 :     while (r < 0) { r = r + Et; t++; }
    1565       48199 :     while (s < 0) { s = s + Ea; t++; }
    1566             : 
    1567             :     /* A^s T^r / p^t */
    1568       48017 :     c = RgXQ_mul(RgXQ_powu(A, s, S->chi), RgXQ_powu(T, r, S->chi), S->chi);
    1569       48018 :     c = RgX_Rg_div(c, powiu(S->p, t));
    1570       48018 :     c = redelt(c, S->psc, S->p);
    1571             :   }
    1572       48018 :   phi = RgX_add(c,  pol_x(varn(S->chi)));
    1573       48017 :   chi = mycaract(S, S->chi, phi, S->psf, S->prc);
    1574       48018 :   if (split_char(S, chi, phi, T0, &nu)) return 1;
    1575             :   /* E_phi = lcm(E_alpha,E_theta) */
    1576       48018 :   S->phi0 = T0;
    1577       48018 :   S->chi = chi;
    1578       48018 :   S->phi = phi;
    1579       48018 :   S->nu = nu; return 0;
    1580             : }
    1581             : 
    1582             : /* Return h^(-degpol(P)) P(x * h) if result is integral, NULL otherwise */
    1583             : static GEN
    1584       59523 : ZX_rescale_inv(GEN P, GEN h)
    1585             : {
    1586       59523 :   long i, l = lg(P);
    1587       59523 :   GEN Q = cgetg(l,t_POL), hi = h;
    1588       59523 :   gel(Q,l-1) = gel(P,l-1);
    1589      172444 :   for (i=l-2; i>=2; i--)
    1590             :   {
    1591             :     GEN r;
    1592      172444 :     gel(Q,i) = dvmdii(gel(P,i), hi, &r);
    1593      172440 :     if (signe(r)) return NULL;
    1594      172440 :     if (i == 2) break;
    1595      112919 :     hi = mulii(hi,h);
    1596             :   }
    1597       59521 :   Q[1] = P[1]; return Q;
    1598             : }
    1599             : 
    1600             : /* x p^-eq nu^-er mod p */
    1601             : static GEN
    1602      299777 : get_gamma(decomp_t *S, GEN x, long eq, long er)
    1603             : {
    1604      299777 :   GEN q, g = x, Dg = powiu(S->p, eq);
    1605      299773 :   long vDg = eq;
    1606      299773 :   if (er)
    1607             :   {
    1608       22635 :     if (!S->invnu)
    1609             :     {
    1610       15724 :       while (gdvd(S->chi, S->nu)) S->nu = RgX_Rg_add(S->nu, S->p);
    1611       15724 :       S->invnu = QXQ_inv(S->nu, S->chi);
    1612       15724 :       S->invnu = redelt_i(S->invnu, S->psc, S->p, &S->Dinvnu, &S->vDinvnu);
    1613             :     }
    1614       22635 :     if (S->Dinvnu) {
    1615       22635 :       Dg = mulii(Dg, powiu(S->Dinvnu, er));
    1616       22635 :       vDg += er * S->vDinvnu;
    1617             :     }
    1618       22635 :     q = mulii(S->p, Dg);
    1619       22635 :     g = ZX_mul(g, FpXQ_powu(S->invnu, er, S->chi, q));
    1620       22635 :     g = FpX_rem(g, S->chi, q);
    1621       22635 :     update_den(S->p, &g, &Dg, &vDg, NULL);
    1622       22634 :     g = centermod(g, mulii(S->p, Dg));
    1623             :   }
    1624      299773 :   if (!is_pm1(Dg)) g = RgX_Rg_div(g, Dg);
    1625      299777 :   return g;
    1626             : }
    1627             : static GEN
    1628      352772 : get_g(decomp_t *S, long Ea, long L, long E, GEN beta, GEN *pchig,
    1629             :       long *peq, long *per)
    1630             : {
    1631             :   long eq, er;
    1632      352772 :   GEN g, chig, chib = NULL;
    1633             :   for(;;) /* at most twice */
    1634             :   {
    1635      359298 :     if (L < 0)
    1636             :     {
    1637       60260 :       chib = ZXQ_charpoly(beta, S->chi, varn(S->chi));
    1638       60260 :       vstar(S->p, chib, &L, &E);
    1639             :     }
    1640      359301 :     eq = L / E; er = L*Ea / E - eq*Ea;
    1641             :     /* floor(L Ea/E) = eq Ea + er */
    1642      359301 :     if (er || !chib)
    1643             :     { /* g might not be an integer ==> chig = NULL */
    1644      299778 :       g = get_gamma(S, beta, eq, er);
    1645      299777 :       chig = mycaract(S, S->chi, g, S->psc, S->prc);
    1646             :     }
    1647             :     else
    1648             :     { /* g = beta/p^eq, special case of the above */
    1649       59523 :       GEN h = powiu(S->p, eq);
    1650       59523 :       g = RgX_Rg_div(beta, h);
    1651       59523 :       chig = ZX_rescale_inv(chib, h); /* chib(x h) / h^N */
    1652       59521 :       if (chig) chig = FpX_red(chig, S->pmf);
    1653             :     }
    1654             :     /* either success or second consecutive failure */
    1655      359299 :     if (chig || chib) break;
    1656             :     /* if g fails the v*-test, v(beta) was wrong. Retry once */
    1657        6526 :     L = -1;
    1658             :   }
    1659      352773 :   *pchig = chig; *peq = eq; *per = er; return g;
    1660             : }
    1661             : 
    1662             : /* return 1 if at least 2 factors mod p ==> chi can be split */
    1663             : static int
    1664      236782 : loop(decomp_t *S, long Ea)
    1665             : {
    1666      236782 :   pari_sp av = avma;
    1667      236782 :   GEN beta = FpXQ_powu(S->nu, Ea, S->chi, S->p);
    1668      236780 :   long N = degpol(S->f), v = varn(S->f);
    1669      236780 :   S->invnu = NULL;
    1670             :   for (;;)
    1671      115990 :   { /* beta tends to a factor of chi */
    1672             :     long L, i, Fg, eq, er;
    1673      352770 :     GEN chig = NULL, d, g, nug;
    1674             : 
    1675      352770 :     if (DEBUGLEVEL>4) err_printf("  beta = %Ps\n", beta);
    1676      352770 :     L = ZpX_resultant_val(S->chi, beta, S->p, S->mf+1);
    1677      352772 :     if (L > S->mf) L = -1; /* from scratch */
    1678      352772 :     g = get_g(S, Ea, L, N, beta, &chig, &eq, &er);
    1679      352773 :     if (DEBUGLEVEL>4) err_printf("  (eq,er) = (%ld,%ld)\n", eq,er);
    1680             :     /* g = beta p^-eq  nu^-er (a unit), chig = charpoly(g) */
    1681      469366 :     if (split_char(S, chig, g,S->phi, &nug)) return 1;
    1682             : 
    1683      232583 :     Fg = degpol(nug);
    1684      232583 :     if (Fg == 1)
    1685             :     { /* frequent special case nug = x - d */
    1686             :       long Le, Ee;
    1687             :       GEN chie, nue, e, pie;
    1688      157935 :       d = negi(gel(nug,2));
    1689      157934 :       chie = RgX_translate(chig, d);
    1690      157935 :       nue = pol_x(v);
    1691      157934 :       e = RgX_Rg_sub(g, d);
    1692      157933 :       pie = getprime(S, e, chie, nue, &Le, &Ee,  0,Ea);
    1693      157934 :       if (pie) return testc2(S, S->nu, Ea, pie, Ee);
    1694             :     }
    1695             :     else
    1696             :     {
    1697       74648 :       long Fa = degpol(S->nu), vdeng;
    1698             :       GEN deng, numg, nume;
    1699       78144 :       if (Fa % Fg) return testb2(S, ulcm(Fa,Fg), g);
    1700             :       /* nu & nug irreducible mod p, deg nug | deg nu. To improve beta, look
    1701             :        * for a root d of nug in Fp[phi] such that v_p(g - d) > 0 */
    1702        8109 :       if (ZX_equal(nug, S->nu))
    1703        5882 :         d = pol_x(v);
    1704             :       else
    1705             :       {
    1706        2227 :         if (!p_is_prime(S)) pari_err_PRIME("FpX_ffisom",S->p);
    1707        2227 :         d = FpX_ffisom(nug, S->nu, S->p);
    1708             :       }
    1709             :       /* write g = numg / deng, e = nume / deng */
    1710        8109 :       numg = QpX_remove_denom(g, S->p, &deng, &vdeng);
    1711       12798 :       for (i = 1; i <= Fg; i++)
    1712             :       {
    1713             :         GEN chie, nue, e;
    1714       12798 :         if (i != 1) d = FpXQ_pow(d, S->p, S->nu, S->p); /* next root */
    1715       12798 :         nume = ZX_sub(numg, ZX_Z_mul(d, deng));
    1716             :         /* test e = nume / deng */
    1717       12798 :         if (ZpX_resultant_val(S->chi, nume, S->p, vdeng*N+1) <= vdeng*N)
    1718        4689 :           continue;
    1719        8109 :         e = RgX_Rg_div(nume, deng);
    1720        8109 :         chie = mycaract(S, S->chi, e, S->psc, S->prc);
    1721        9569 :         if (split_char(S, chie, e,S->phi, &nue)) return 1;
    1722        6073 :         if (RgX_is_monomial(nue))
    1723             :         { /* v_p(e) = v_p(g - d) > 0 */
    1724             :           long Le, Ee;
    1725             :           GEN pie;
    1726        6073 :           pie = getprime(S, e, chie, nue, &Le, &Ee,  0,Ea);
    1727        6073 :           if (pie) return testc2(S, S->nu, Ea, pie, Ee);
    1728        4613 :           break;
    1729             :         }
    1730             :       }
    1731        4613 :       if (i > Fg)
    1732             :       {
    1733           0 :         if (!p_is_prime(S)) pari_err_PRIME("nilord",S->p);
    1734           0 :         pari_err_BUG("nilord (no root)");
    1735             :       }
    1736             :     }
    1737      115989 :     if (eq) d = gmul(d, powiu(S->p,  eq));
    1738      115988 :     if (er) d = gmul(d, gpowgs(S->nu, er));
    1739      115988 :     beta = gsub(beta, d);
    1740             : 
    1741      115990 :     if (gc_needed(av,1))
    1742             :     {
    1743           0 :       if (DEBUGMEM > 1) pari_warn(warnmem, "nilord");
    1744           0 :       gerepileall(av, S->invnu? 6: 4, &beta, &(S->precns), &(S->ns), &(S->nsf), &(S->invnu), &(S->Dinvnu));
    1745             :     }
    1746             :   }
    1747             : }
    1748             : 
    1749             : /* E and F cannot decrease; return 1 if O = Zp[phi], 2 if we can get a
    1750             :  * decomposition and 0 otherwise */
    1751             : static long
    1752      391498 : progress(decomp_t *S, GEN *ppa, long *pE)
    1753             : {
    1754      391498 :   long E = *pE, F;
    1755      391498 :   GEN pa = *ppa;
    1756      391498 :   S->phi0 = NULL; /* no delayed composition */
    1757             :   for(;;)
    1758        2826 :   {
    1759             :     long l, La, Ea; /* N.B If E = 0, getprime cannot return NULL */
    1760      394324 :     GEN pia  = getprime(S, NULL, S->chi, S->nu, &La, &Ea, E,0);
    1761      394340 :     if (pia) { /* success, we break out in THIS loop */
    1762      391514 :       pa = (typ(pia) == t_POL)? RgX_RgXQ_eval(pia, S->phi, S->f): pia;
    1763      391516 :       E = Ea;
    1764      391516 :       if (La == 1) break; /* no need to change phi so that nu = pia */
    1765             :     }
    1766             :     /* phi += prime elt */
    1767       64689 :     S->phi = typ(pa) == t_INT? RgX_Rg_add_shallow(S->phi, pa)
    1768      159946 :                              : RgX_add(S->phi, pa);
    1769             :     /* recompute char. poly. chi from scratch */
    1770      159945 :     S->chi = mycaract(S, S->f, S->phi, S->psf, S->pdf);
    1771      159944 :     S->nu = get_nu(S->chi, S->p, &l);
    1772      159943 :     if (l > 1) return 2;
    1773      159943 :     if (!update_phi(S)) return 1; /* unramified */
    1774      159945 :     if (pia) break;
    1775             :   }
    1776      391515 :   *pE = E; *ppa = pa; F = degpol(S->nu);
    1777      391512 :   if (DEBUGLEVEL>4) err_printf("  (E, F) = (%ld,%ld)\n", E, F);
    1778      391512 :   if (E * F == degpol(S->f)) return 1;
    1779      236780 :   if (loop(S, E)) return 2;
    1780      114557 :   if (!update_phi(S)) return 1;
    1781       70446 :   return 0;
    1782             : }
    1783             : 
    1784             : /* flag != 0 iff we're looking for the p-adic factorization,
    1785             :    in which case it is the p-adic precision we want */
    1786             : static GEN
    1787      452095 : maxord_i(decomp_t *S, GEN p, GEN f, long mf, GEN w, long flag)
    1788             : {
    1789      452095 :   long oE, n = lg(w)-1; /* factor of largest degree */
    1790      452095 :   GEN opa, D = ZpX_reduced_resultant_fast(f, ZX_deriv(f), p, mf);
    1791      452098 :   S->pisprime = -1;
    1792      452098 :   S->p = p;
    1793      452098 :   S->mf = mf;
    1794      452098 :   S->nu = gel(w,n);
    1795      452098 :   S->df = Z_pval(D, p);
    1796      452100 :   S->pdf = powiu(p, S->df);
    1797      452080 :   S->phi = pol_x(varn(f));
    1798      452088 :   S->chi = S->f = f;
    1799      452088 :   if (n > 1) return Decomp(S, flag); /* FIXME: use bezout_lift_fact */
    1800             : 
    1801      321057 :   if (DEBUGLEVEL>4)
    1802           0 :     err_printf("  entering Nilord: %Ps^%ld\n  f = %Ps, nu = %Ps\n",
    1803             :                p, S->df, S->f, S->nu);
    1804      321057 :   else if (DEBUGLEVEL>2) err_printf("  entering Nilord\n");
    1805      321057 :   S->psf = S->psc = mulii(sqri(D), p);
    1806      321047 :   S->vpsf = S->vpsc = 2*S->df + 1;
    1807      321047 :   S->prc = mulii(D, p);
    1808      321045 :   S->chi = FpX_red(S->f, S->psc);
    1809      321057 :   S->pmf = powiu(p, S->mf+1);
    1810      321055 :   S->precns = NULL;
    1811      321055 :   for(opa = NULL, oE = 0;;)
    1812       70446 :   {
    1813      391501 :     long n = progress(S, &opa, &oE);
    1814      391514 :     if (n == 1) return flag? NULL: dbasis(p, S->f, S->mf, S->phi, S->chi);
    1815      192674 :     if (n == 2) return Decomp(S, flag);
    1816             :   }
    1817             : }
    1818             : 
    1819             : static int
    1820         812 : expo_is_squarefree(GEN e)
    1821             : {
    1822         812 :   long i, l = lg(e);
    1823        1183 :   for (i=1; i<l; i++)
    1824         945 :     if (e[i] != 1) return 0;
    1825         238 :   return 1;
    1826             : }
    1827             : /* pure round 4 */
    1828             : static GEN
    1829         770 : ZpX_round4(GEN f, GEN p, GEN w, long prec)
    1830             : {
    1831             :   decomp_t S;
    1832         770 :   GEN L = maxord_i(&S, p, f, ZpX_disc_val(f,p), w, prec);
    1833         770 :   return L? L: mkvec(f);
    1834             : }
    1835             : /* f a squarefree ZX with leading_coeff 1, degree > 0. Return list of
    1836             :  * irreducible factors in Zp[X] (computed mod p^prec) */
    1837             : static GEN
    1838        1071 : ZpX_monic_factor_squarefree(GEN f, GEN p, long prec)
    1839             : {
    1840        1071 :   pari_sp av = avma;
    1841             :   GEN L, fa, w, e;
    1842             :   long i, l;
    1843        1071 :   if (degpol(f) == 1) return mkvec(f);
    1844         812 :   fa = FpX_factor(f,p); w = gel(fa,1); e = gel(fa,2);
    1845             :   /* no repeated factors: Hensel lift */
    1846         812 :   if (expo_is_squarefree(e)) return ZpX_liftfact(f, w, powiu(p,prec), p, prec);
    1847         574 :   l = lg(w);
    1848         574 :   if (l == 2)
    1849             :   {
    1850         371 :     L = ZpX_round4(f,p,w,prec);
    1851         371 :     if (lg(L) == 2) { set_avma(av); return mkvec(f); }
    1852             :   }
    1853             :   else
    1854             :   { /* >= 2 factors mod p: partial Hensel lift */
    1855         203 :     GEN D = ZpX_reduced_resultant_fast(f, ZX_deriv(f), p, ZpX_disc_val(f,p));
    1856         203 :     long r = maxss(2*Z_pval(D,p)+1, prec);
    1857         203 :     GEN W = cgetg(l, t_VEC);
    1858         665 :     for (i = 1; i < l; i++)
    1859         462 :       gel(W,i) = e[i] == 1? gel(w,i): FpX_powu(gel(w,i), e[i], p);
    1860         203 :     L = ZpX_liftfact(f, W, powiu(p,r), p, r);
    1861         665 :     for (i = 1; i < l; i++)
    1862         462 :       gel(L,i) = e[i] == 1? mkvec(gel(L,i))
    1863         462 :                           : ZpX_round4(gel(L,i), p, mkvec(gel(w,i)), prec);
    1864         203 :     L = shallowconcat1(L);
    1865             :   }
    1866         343 :   return gerepilecopy(av, L);
    1867             : }
    1868             : 
    1869             : /* assume T a ZX with leading_coeff 1, degree > 0 */
    1870             : GEN
    1871         490 : ZpX_monic_factor(GEN T, GEN p, long prec)
    1872             : {
    1873             :   GEN Q, P, E, F;
    1874             :   long L, l, i, v;
    1875             : 
    1876         490 :   if (degpol(T) == 1) return mkmat2(mkcol(T), mkcol(gen_1));
    1877         490 :   v = ZX_valrem(T, &T);
    1878         490 :   Q = ZX_squff(T, &F); l = lg(Q); L = v? l + 1: l;
    1879         490 :   P = cgetg(L, t_VEC);
    1880         490 :   E = cgetg(L, t_VEC);
    1881         987 :   for (i = 1; i < l; i++)
    1882             :   {
    1883         497 :     GEN w = ZpX_monic_factor_squarefree(gel(Q,i), p, prec);
    1884         497 :     gel(P,i) = w; settyp(w, t_COL);
    1885         497 :     gel(E,i) = const_col(lg(w)-1, utoipos(F[i]));
    1886             :   }
    1887         490 :   if (v) { gel(P,i) = pol_x(varn(T)); gel(E,i) = utoipos(v); }
    1888         490 :   return mkmat2(shallowconcat1(P), shallowconcat1(E));
    1889             : }
    1890             : 
    1891             : /* DT = multiple of disc(T) or NULL
    1892             :  * Return a multiple of the denominator of an algebraic integer (in Q[X]/(T))
    1893             :  * when expressed in terms of the power basis */
    1894             : GEN
    1895       44065 : indexpartial(GEN T, GEN DT)
    1896             : {
    1897       44065 :   pari_sp av = avma;
    1898             :   long i, nb;
    1899       44065 :   GEN fa, E, P, U, res = gen_1, dT = ZX_deriv(T);
    1900             : 
    1901       44062 :   if (!DT) DT = ZX_disc(T);
    1902       44062 :   fa = absZ_factor_limit_strict(DT, 0, &U);
    1903       44064 :   P = gel(fa,1);
    1904       44064 :   E = gel(fa,2); nb = lg(P)-1;
    1905      211878 :   for (i = 1; i <= nb; i++)
    1906             :   {
    1907      167822 :     long e = itou(gel(E,i)), e2 = e >> 1;
    1908      167821 :     GEN p = gel(P,i), q = p;
    1909      167821 :     if (e2 >= 2) q = ZpX_reduced_resultant_fast(T, dT, p, e2);
    1910      167821 :     res = mulii(res, q);
    1911             :   }
    1912       44056 :   if (U)
    1913             :   {
    1914        1896 :     long e = itou(gel(U,2)), e2 = e >> 1;
    1915        1896 :     GEN p = gel(U,1), q = powiu(p, odd(e)? e2+1: e2);
    1916        1896 :     res = mulii(res, q);
    1917             :   }
    1918       44056 :   return gerepileuptoint(av,res);
    1919             : }
    1920             : 
    1921             : /*******************************************************************/
    1922             : /*                                                                 */
    1923             : /*    2-ELT REPRESENTATION FOR PRIME IDEALS (dividing index)       */
    1924             : /*                                                                 */
    1925             : /*******************************************************************/
    1926             : /* to compute norm of elt in basis form */
    1927             : typedef struct {
    1928             :   long r1;
    1929             :   GEN M;  /* via embed_norm */
    1930             : 
    1931             :   GEN D, w, T; /* via resultant if M = NULL */
    1932             : } norm_S;
    1933             : 
    1934             : static GEN
    1935      501847 : get_norm(norm_S *S, GEN a)
    1936             : {
    1937      501847 :   if (S->M)
    1938             :   {
    1939             :     long e;
    1940      500441 :     GEN N = grndtoi( embed_norm(RgM_RgC_mul(S->M, a), S->r1), &e );
    1941      500479 :     if (e > -5) pari_err_PREC( "get_norm");
    1942      500479 :     return N;
    1943             :   }
    1944        1406 :   if (S->w) a = RgV_RgC_mul(S->w, a);
    1945        1406 :   return ZX_resultant_all(S->T, a, S->D, 0);
    1946             : }
    1947             : static void
    1948      213342 : init_norm(norm_S *S, GEN nf, GEN p)
    1949             : {
    1950      213342 :   GEN T = nf_get_pol(nf), M = nf_get_M(nf);
    1951      213343 :   long N = degpol(T), ex = gexpo(M) + gexpo(mului(8 * N, p));
    1952             : 
    1953      213350 :   S->r1 = nf_get_r1(nf);
    1954      213352 :   if (N * ex <= gprecision(M) - 20)
    1955             :   { /* enough prec to use embed_norm */
    1956      213173 :     S->M = M;
    1957      213173 :     S->D = NULL;
    1958      213173 :     S->w = NULL;
    1959      213173 :     S->T = NULL;
    1960             :   }
    1961             :   else
    1962             :   {
    1963         191 :     GEN w = leafcopy(nf_get_zkprimpart(nf)), D = nf_get_zkden(nf), Dp = sqri(p);
    1964             :     long i;
    1965         191 :     if (!equali1(D))
    1966             :     {
    1967         191 :       GEN w1 = D;
    1968         191 :       long v = Z_pval(D, p);
    1969         191 :       D = powiu(p, v);
    1970         191 :       Dp = mulii(D, Dp);
    1971         191 :       gel(w, 1) = remii(w1, Dp);
    1972             :     }
    1973        3969 :     for (i=2; i<=N; i++) gel(w,i) = FpX_red(gel(w,i), Dp);
    1974         191 :     S->M = NULL;
    1975         191 :     S->D = D;
    1976         191 :     S->w = w;
    1977         191 :     S->T = T;
    1978             :   }
    1979      213364 : }
    1980             : /* f = f(pr/p), q = p^(f+1), a in pr.
    1981             :  * Return 1 if v_pr(a) = 1, and 0 otherwise */
    1982             : static int
    1983      501845 : is_uniformizer(GEN a, GEN q, norm_S *S) { return !dvdii(get_norm(S,a), q); }
    1984             : 
    1985             : /* Return x * y, x, y are t_MAT (Fp-basis of in O_K/p), assume (x,y)=1.
    1986             :  * Either x or y may be NULL (= O_K), not both */
    1987             : static GEN
    1988      699238 : mul_intersect(GEN x, GEN y, GEN p)
    1989             : {
    1990      699238 :   if (!x) return y;
    1991      372683 :   if (!y) return x;
    1992      263831 :   return FpM_intersect_i(x, y, p);
    1993             : }
    1994             : /* Fp-basis of (ZK/pr): applied to the primes found in primedec_aux()
    1995             :  * true nf */
    1996             : static GEN
    1997      305643 : Fp_basis(GEN nf, GEN pr)
    1998             : {
    1999             :   long i, j, l;
    2000             :   GEN x, y;
    2001             :   /* already in basis form (from Buchman-Lenstra) ? */
    2002      305643 :   if (typ(pr) == t_MAT) return pr;
    2003             :   /* ordinary prid (from Kummer) */
    2004       72882 :   x = pr_hnf(nf, pr);
    2005       72882 :   l = lg(x);
    2006       72882 :   y = cgetg(l, t_MAT);
    2007      614394 :   for (i=j=1; i<l; i++)
    2008      541512 :     if (gequal1(gcoeff(x,i,i))) gel(y,j++) = gel(x,i);
    2009       72882 :   setlg(y, j); return y;
    2010             : }
    2011             : /* Let Ip = prod_{ P | p } P be the p-radical. The list L contains the
    2012             :  * P (mod Ip) seen as sub-Fp-vector spaces of ZK/Ip.
    2013             :  * Return the list of (Ip / P) (mod Ip).
    2014             :  * N.B: All ideal multiplications are computed as intersections of Fp-vector
    2015             :  * spaces. true nf */
    2016             : static GEN
    2017      213370 : get_LV(GEN nf, GEN L, GEN p, long N)
    2018             : {
    2019      213370 :   long i, l = lg(L)-1;
    2020             :   GEN LV, LW, A, B;
    2021             : 
    2022      213370 :   LV = cgetg(l+1, t_VEC);
    2023      213368 :   if (l == 1) { gel(LV,1) = matid(N); return LV; }
    2024      108851 :   LW = cgetg(l+1, t_VEC);
    2025      414494 :   for (i=1; i<=l; i++) gel(LW,i) = Fp_basis(nf, gel(L,i));
    2026             : 
    2027             :   /* A[i] = L[1]...L[i-1], i = 2..l */
    2028      108852 :   A = cgetg(l+1, t_VEC); gel(A,1) = NULL;
    2029      305645 :   for (i=1; i < l; i++) gel(A,i+1) = mul_intersect(gel(A,i), gel(LW,i), p);
    2030             :   /* B[i] = L[i+1]...L[l], i = 1..(l-1) */
    2031      108852 :   B = cgetg(l+1, t_VEC); gel(B,l) = NULL;
    2032      305648 :   for (i=l; i>=2; i--) gel(B,i-1) = mul_intersect(gel(B,i), gel(LW,i), p);
    2033      414500 :   for (i=1; i<=l; i++) gel(LV,i) = mul_intersect(gel(A,i), gel(B,i), p);
    2034      108852 :   return LV;
    2035             : }
    2036             : 
    2037             : static void
    2038           0 : errprime(GEN p) { pari_err_PRIME("idealprimedec",p); }
    2039             : 
    2040             : /* P = Fp-basis (over O_K/p) for pr.
    2041             :  * V = Z-basis for I_p/pr. ramif != 0 iff some pr|p is ramified.
    2042             :  * Return a p-uniformizer for pr. Assume pr not inert, i.e. m > 0 */
    2043             : static GEN
    2044      295531 : uniformizer(GEN nf, norm_S *S, GEN P, GEN V, GEN p, int ramif)
    2045             : {
    2046      295531 :   long i, l, f, m = lg(P)-1, N = nf_get_degree(nf);
    2047             :   GEN u, Mv, x, q;
    2048             : 
    2049      295531 :   f = N - m; /* we want v_p(Norm(x)) = p^f */
    2050      295531 :   q = powiu(p,f+1);
    2051             : 
    2052      295504 :   u = FpM_FpC_invimage(shallowconcat(P, V), col_ei(N,1), p);
    2053      295530 :   setlg(u, lg(P));
    2054      295530 :   u = centermod(ZM_ZC_mul(P, u), p);
    2055      295520 :   if (is_uniformizer(u, q, S)) return u;
    2056      157993 :   if (signe(gel(u,1)) <= 0) /* make sure u[1] in ]-p,p] */
    2057      134137 :     gel(u,1) = addii(gel(u,1), p); /* try u + p */
    2058             :   else
    2059       23856 :     gel(u,1) = subii(gel(u,1), p); /* try u - p */
    2060      157986 :   if (!ramif || is_uniformizer(u, q, S)) return u;
    2061             : 
    2062             :   /* P/p ramified, u in P^2, not in Q for all other Q|p */
    2063       86063 :   Mv = zk_multable(nf, Z_ZC_sub(gen_1,u));
    2064       86069 :   l = lg(P);
    2065      115552 :   for (i=1; i<l; i++)
    2066             :   {
    2067      115552 :     x = centermod(ZC_add(u, ZM_ZC_mul(Mv, gel(P,i))), p);
    2068      115552 :     if (is_uniformizer(x, q, S)) return x;
    2069             :   }
    2070           0 :   errprime(p);
    2071             :   return NULL; /* LCOV_EXCL_LINE */
    2072             : }
    2073             : 
    2074             : /*******************************************************************/
    2075             : /*                                                                 */
    2076             : /*                   BUCHMANN-LENSTRA ALGORITHM                    */
    2077             : /*                                                                 */
    2078             : /*******************************************************************/
    2079             : static GEN
    2080     4245463 : mk_pr(GEN p, GEN u, long e, long f, GEN t)
    2081     4245463 : { return mkvec5(p, u, utoipos(e), utoipos(f), t); }
    2082             : 
    2083             : /* nf a true nf, u in Z[X]/(T); pr = p Z_K + u Z_K of ramification index e */
    2084             : GEN
    2085     3790534 : idealprimedec_kummer(GEN nf,GEN u,long e,GEN p)
    2086             : {
    2087     3790534 :   GEN t, T = nf_get_pol(nf);
    2088     3790534 :   long f = degpol(u), N = degpol(T);
    2089             : 
    2090     3790509 :   if (f == N)
    2091             :   { /* inert */
    2092      622446 :     u = scalarcol_shallow(p,N);
    2093      622449 :     t = gen_1;
    2094             :   }
    2095             :   else
    2096             :   {
    2097     3168063 :     t = centermod(poltobasis(nf, FpX_div(T, u, p)), p);
    2098     3167774 :     u = centermod(poltobasis(nf, u), p);
    2099     3167793 :     if (e == 1)
    2100             :     { /* make sure v_pr(u) = 1 (automatic if e>1) */
    2101     2891996 :       GEN cw, w = Q_primitive_part(nf_to_scalar_or_alg(nf, u), &cw);
    2102     2892199 :       long v = cw? f - Q_pval(cw, p) * N: f;
    2103     2892198 :       if (ZpX_resultant_val(T, w, p, v + 1) > v)
    2104             :       {
    2105      107858 :         GEN c = gel(u,1);
    2106      107858 :         gel(u,1) = signe(c) > 0? subii(c, p): addii(c, p);
    2107             :       }
    2108             :     }
    2109     3168081 :     t = zk_multable(nf, t);
    2110             :   }
    2111     3790391 :   return mk_pr(p,u,e,f,t);
    2112             : }
    2113             : 
    2114             : typedef struct {
    2115             :   GEN nf, p;
    2116             :   long I;
    2117             : } eltmod_muldata;
    2118             : 
    2119             : static GEN
    2120      816594 : sqr_mod(void *data, GEN x)
    2121             : {
    2122      816594 :   eltmod_muldata *D = (eltmod_muldata*)data;
    2123      816594 :   return FpC_red(nfsqri(D->nf, x), D->p);
    2124             : }
    2125             : static GEN
    2126      344185 : ei_msqr_mod(void *data, GEN x)
    2127             : {
    2128      344185 :   GEN x2 = sqr_mod(data, x);
    2129      344178 :   eltmod_muldata *D = (eltmod_muldata*)data;
    2130      344178 :   return FpC_red(zk_ei_mul(D->nf, x2, D->I), D->p);
    2131             : }
    2132             : /* nf a true nf; compute lift(nf.zk[I]^p mod p) */
    2133             : static GEN
    2134      735175 : pow_ei_mod_p(GEN nf, long I, GEN p)
    2135             : {
    2136      735175 :   pari_sp av = avma;
    2137             :   eltmod_muldata D;
    2138      735175 :   long N = nf_get_degree(nf);
    2139      735177 :   GEN y = col_ei(N,I);
    2140      735182 :   if (I == 1) return y;
    2141      519486 :   D.nf = nf;
    2142      519486 :   D.p = p;
    2143      519486 :   D.I = I;
    2144      519486 :   y = gen_pow_fold(y, p, (void*)&D, &sqr_mod, &ei_msqr_mod);
    2145      519488 :   return gerepileupto(av,y);
    2146             : }
    2147             : 
    2148             : /* nf a true nf; return a Z basis of Z_K's p-radical, phi = x--> x^p-x */
    2149             : static GEN
    2150      213367 : pradical(GEN nf, GEN p, GEN *phi)
    2151             : {
    2152      213367 :   long i, N = nf_get_degree(nf);
    2153             :   GEN q,m,frob,rad;
    2154             : 
    2155             :   /* matrix of Frob: x->x^p over Z_K/p */
    2156      213365 :   frob = cgetg(N+1,t_MAT);
    2157      938906 :   for (i=1; i<=N; i++) gel(frob,i) = pow_ei_mod_p(nf,i,p);
    2158             : 
    2159      213363 :   m = frob; q = p;
    2160      302832 :   while (abscmpiu(q,N) < 0) { q = mulii(q,p); m = FpM_mul(m, frob, p); }
    2161      213360 :   rad = FpM_ker(m, p); /* m = Frob^k, s.t p^k >= N */
    2162      938859 :   for (i=1; i<=N; i++) gcoeff(frob,i,i) = subiu(gcoeff(frob,i,i), 1);
    2163      213326 :   *phi = frob; return rad;
    2164             : }
    2165             : 
    2166             : /* return powers of a: a^0, ... , a^d,  d = dim A */
    2167             : static GEN
    2168      158386 : get_powers(GEN mul, GEN p)
    2169             : {
    2170      158386 :   long i, d = lgcols(mul);
    2171      158386 :   GEN z, pow = cgetg(d+2,t_MAT), P = pow+1;
    2172             : 
    2173      158386 :   gel(P,0) = scalarcol_shallow(gen_1, d-1);
    2174      158389 :   z = gel(mul,1);
    2175      751512 :   for (i=1; i<=d; i++)
    2176             :   {
    2177      593130 :     gel(P,i) = z; /* a^i */
    2178      593130 :     if (i!=d) z = FpM_FpC_mul(mul, z, p);
    2179             :   }
    2180      158382 :   return pow;
    2181             : }
    2182             : 
    2183             : /* minimal polynomial of a in A (dim A = d).
    2184             :  * mul = multiplication table by a in A */
    2185             : static GEN
    2186      116869 : pol_min(GEN mul, GEN p)
    2187             : {
    2188      116869 :   pari_sp av = avma;
    2189      116869 :   GEN z = FpM_deplin(get_powers(mul, p), p);
    2190      116869 :   return gerepilecopy(av, RgV_to_RgX(z,0));
    2191             : }
    2192             : 
    2193             : static GEN
    2194      409837 : get_pr(GEN nf, norm_S *S, GEN p, GEN P, GEN V, int ramif, long N, long flim)
    2195             : {
    2196             :   GEN u, t;
    2197             :   long e, f;
    2198             : 
    2199      409837 :   if (typ(P) == t_VEC)
    2200             :   { /* already done (Kummer) */
    2201       72881 :     f = pr_get_f(P);
    2202       72881 :     if (flim > 0 && f > flim) return NULL;
    2203       71989 :     if (flim == -2) return (GEN)f;
    2204       71989 :     return P;
    2205             :   }
    2206      336956 :   f = N - (lg(P)-1);
    2207      336956 :   if (flim > 0 && f > flim) return NULL;
    2208      335357 :   if (flim == -2) return (GEN)f;
    2209             :   /* P = (p,u) prime. t is an anti-uniformizer: Z_K + t/p Z_K = P^(-1),
    2210             :    * so that v_P(t) = e(P/p)-1 */
    2211      334986 :   if (f == N) {
    2212       39458 :     u = scalarcol_shallow(p,N);
    2213       39458 :     t = gen_1;
    2214       39458 :     e = 1;
    2215             :   } else {
    2216             :     GEN mt;
    2217      295528 :     u = uniformizer(nf, S, P, V, p, ramif);
    2218      295491 :     t = FpM_deplin(zk_multable(nf,u), p);
    2219      295532 :     mt = zk_multable(nf, t);
    2220      295525 :     e = ramif? 1 + ZC_nfval(t,mk_pr(p,u,0,0,mt)): 1;
    2221      295510 :     t = mt;
    2222             :   }
    2223      334968 :   return mk_pr(p,u,e,f,t);
    2224             : }
    2225             : 
    2226             : /* true nf */
    2227             : static GEN
    2228      213368 : primedec_end(GEN nf, GEN L, GEN p, long flim)
    2229             : {
    2230      213368 :   long i, j, l = lg(L), N = nf_get_degree(nf);
    2231      213370 :   GEN LV = get_LV(nf, L,p,N);
    2232      213370 :   int ramif = dvdii(nf_get_disc(nf), p);
    2233      213341 :   norm_S S; init_norm(&S, nf, p);
    2234      622795 :   for (i = j = 1; i < l; i++)
    2235             :   {
    2236      409836 :     GEN P = get_pr(nf, &S, p, gel(L,i), gel(LV,i), ramif, N, flim);
    2237      409823 :     if (!P) continue;
    2238      407332 :     gel(L,j++) = P;
    2239      407332 :     if (flim == -1) return P;
    2240             :   }
    2241      212959 :   setlg(L, j); return L;
    2242             : }
    2243             : 
    2244             : /* prime ideal decomposition of p; if flim>0, restrict to f(P,p) <= flim
    2245             :  * if flim = -1 return only the first P
    2246             :  * if flim = -2 return only the f(P/p) in a t_VECSMALL; true nf */
    2247             : static GEN
    2248     2699584 : primedec_aux(GEN nf, GEN p, long flim)
    2249             : {
    2250     2699584 :   const long TYP = (flim == -2)? t_VECSMALL: t_VEC;
    2251     2699584 :   GEN E, F, L, Ip, phi, f, g, h, UN, T = nf_get_pol(nf);
    2252             :   long i, k, c, iL, N;
    2253             :   int kummer;
    2254             : 
    2255     2699589 :   F = FpX_factor(T, p);
    2256     2699734 :   E = gel(F,2);
    2257     2699734 :   F = gel(F,1);
    2258             : 
    2259     2699734 :   k = lg(F); if (k == 1) errprime(p);
    2260     2699734 :   if ( !dvdii(nf_get_index(nf),p) ) /* p doesn't divide index */
    2261             :   {
    2262     2484780 :     L = cgetg(k, TYP);
    2263     6058928 :     for (i=1; i<k; i++)
    2264             :     {
    2265     4129612 :       GEN t = gel(F,i);
    2266     4129612 :       long f = degpol(t);
    2267     4129627 :       if (flim > 0 && f > flim) { setlg(L, i); break; }
    2268     3578830 :       if (flim == -2)
    2269           0 :         L[i] = f;
    2270             :       else
    2271     3578830 :         gel(L,i) = idealprimedec_kummer(nf, t, E[i],p);
    2272     3578915 :       if (flim == -1) return gel(L,1);
    2273             :     }
    2274     2480112 :     return L;
    2275             :   }
    2276             : 
    2277      214759 :   kummer = 0;
    2278      214759 :   g = FpXV_prod(F, p);
    2279      214752 :   h = FpX_div(T,g,p);
    2280      214756 :   f = FpX_red(ZX_Z_divexact(ZX_sub(ZX_mul(g,h), T), p), p);
    2281             : 
    2282      214729 :   N = degpol(T);
    2283      214737 :   L = cgetg(N+1,TYP);
    2284      214752 :   iL = 1;
    2285      522822 :   for (i=1; i<k; i++)
    2286      309456 :     if (E[i] == 1 || signe(FpX_rem(f,gel(F,i),p)))
    2287       72882 :     {
    2288       74270 :       GEN t = gel(F,i);
    2289       74270 :       kummer = 1;
    2290       74270 :       gel(L,iL++) = idealprimedec_kummer(nf, t, E[i],p);
    2291       74275 :       if (flim == -1) return gel(L,1);
    2292             :     }
    2293             :     else /* F[i] | (f,g,h), happens at least once by Dedekind criterion */
    2294      235188 :       E[i] = 0;
    2295             : 
    2296             :   /* phi matrix of x -> x^p - x in algebra Z_K/p */
    2297      213366 :   Ip = pradical(nf,p,&phi);
    2298             : 
    2299             :   /* split etale algebra Z_K / (p,Ip) */
    2300      213347 :   h = cgetg(N+1,t_VEC);
    2301      213354 :   if (kummer)
    2302             :   { /* split off Kummer factors */
    2303       46128 :     GEN mb, b = NULL;
    2304      172592 :     for (i=1; i<k; i++)
    2305      126464 :       if (!E[i]) b = b? FpX_mul(b, gel(F,i), p): gel(F,i);
    2306       46128 :     if (!b) errprime(p);
    2307       46128 :     b = FpC_red(poltobasis(nf,b), p);
    2308       46127 :     mb = FpM_red(zk_multable(nf,b), p);
    2309             :     /* Fp-base of ideal (Ip, b) in ZK/p */
    2310       46126 :     gel(h,1) = FpM_image(shallowconcat(mb,Ip), p);
    2311             :   }
    2312             :   else
    2313      167226 :     gel(h,1) = Ip;
    2314             : 
    2315      213354 :   UN = col_ei(N, 1);
    2316      465808 :   for (c=1; c; c--)
    2317             :   { /* Let A:= (Z_K/p) / Ip etale; split A2 := A / Im H ~ Im M2
    2318             :        H * ? + M2 * Mi2 = Id_N ==> M2 * Mi2 projector A --> A2 */
    2319      252439 :     GEN M, Mi, M2, Mi2, phi2, mat1, H = gel(h,c); /* maximal rank */
    2320      252439 :     long dim, r = lg(H)-1;
    2321             : 
    2322      252439 :     M   = FpM_suppl(shallowconcat(H,UN), p);
    2323      252438 :     Mi  = FpM_inv(M, p);
    2324      252436 :     M2  = vecslice(M, r+1,N); /* M = (H|M2) invertible */
    2325      252437 :     Mi2 = rowslice(Mi,r+1,N);
    2326             :     /* FIXME: FpM_mul(,M2) could be done with vecpermute */
    2327      252440 :     phi2 = FpM_mul(Mi2, FpM_mul(phi,M2, p), p);
    2328      252434 :     mat1 = FpM_ker(phi2, p);
    2329      252442 :     dim = lg(mat1)-1; /* A2 product of 'dim' fields */
    2330      252442 :     if (dim > 1)
    2331             :     { /* phi2 v = 0 => a = M2 v in Ker phi, a not in Fp.1 + H */
    2332      116869 :       GEN R, a, mula, mul2, v = gel(mat1,2);
    2333             :       long n;
    2334             : 
    2335      116869 :       a = FpM_FpC_mul(M2,v, p); /* not a scalar */
    2336      116866 :       mula = FpM_red(zk_multable(nf,a), p);
    2337      116862 :       mul2 = FpM_mul(Mi2, FpM_mul(mula,M2, p), p);
    2338      116869 :       R = FpX_roots(pol_min(mul2,p), p); /* totally split mod p */
    2339      116869 :       n = lg(R)-1;
    2340      357652 :       for (i=1; i<=n; i++)
    2341             :       {
    2342      240783 :         GEN I = RgM_Rg_sub_shallow(mula, gel(R,i));
    2343      240783 :         gel(h,c++) = FpM_image(shallowconcat(H, I), p);
    2344             :       }
    2345      116869 :       if (n == dim)
    2346      301058 :         for (i=1; i<=n; i++) gel(L,iL++) = gel(h,--c);
    2347             :     }
    2348             :     else /* A2 field ==> H maximal, f = N-r = dim(A2) */
    2349      135573 :       gel(L,iL++) = H;
    2350             :   }
    2351      213369 :   setlg(L, iL);
    2352      213368 :   return primedec_end(nf, L, p, flim);
    2353             : }
    2354             : 
    2355             : GEN
    2356     2692681 : idealprimedec_limit_f(GEN nf, GEN p, long f)
    2357             : {
    2358     2692681 :   pari_sp av = avma;
    2359             :   GEN v;
    2360     2692681 :   if (typ(p) != t_INT) pari_err_TYPE("idealprimedec",p);
    2361     2692681 :   if (f < 0) pari_err_DOMAIN("idealprimedec", "f", "<", gen_0, stoi(f));
    2362     2692681 :   v = primedec_aux(checknf(nf), p, f);
    2363     2692634 :   v = gen_sort(v, (void*)&cmp_prime_over_p, &cmp_nodata);
    2364     2692671 :   return gerepileupto(av,v);
    2365             : }
    2366             : /* true nf */
    2367             : GEN
    2368        6545 : idealprimedec_galois(GEN nf, GEN p)
    2369             : {
    2370        6545 :   pari_sp av = avma;
    2371        6545 :   GEN v = primedec_aux(nf, p, -1);
    2372        6545 :   return gerepilecopy(av,v);
    2373             : }
    2374             : /* true nf */
    2375             : GEN
    2376         364 : idealprimedec_degrees(GEN nf, GEN p)
    2377             : {
    2378         364 :   pari_sp av = avma;
    2379         364 :   GEN v = primedec_aux(nf, p, -2);
    2380         364 :   vecsmall_sort(v); return gerepileuptoleaf(av, v);
    2381             : }
    2382             : GEN
    2383      504832 : idealprimedec_limit_norm(GEN nf, GEN p, GEN B)
    2384      504832 : { return idealprimedec_limit_f(nf, p, logint(B,p)); }
    2385             : GEN
    2386     1269890 : idealprimedec(GEN nf, GEN p)
    2387     1269890 : { return idealprimedec_limit_f(nf, p, 0); }
    2388             : GEN
    2389       38528 : nf_pV_to_prV(GEN nf, GEN P)
    2390             : {
    2391             :   long i, l;
    2392       38528 :   GEN Q = cgetg_copy(P,&l);
    2393       38528 :   if (l == 1) return Q;
    2394       89523 :   for (i = 1; i < l; i++) gel(Q,i) = idealprimedec(nf, gel(P,i));
    2395       26572 :   return shallowconcat1(Q);
    2396             : }
    2397             : 
    2398             : /* return [Fp[x]: Fp] */
    2399             : static long
    2400        4109 : ffdegree(GEN x, GEN frob, GEN p)
    2401             : {
    2402        4109 :   pari_sp av = avma;
    2403        4109 :   long d, f = lg(frob)-1;
    2404        4109 :   GEN y = x;
    2405             : 
    2406       13209 :   for (d=1; d < f; d++)
    2407             :   {
    2408       10878 :     y = FpM_FpC_mul(frob, y, p);
    2409       10878 :     if (ZV_equal(y, x)) break;
    2410             :   }
    2411        4109 :   return gc_long(av,d);
    2412             : }
    2413             : 
    2414             : static GEN
    2415       91427 : lift_to_zk(GEN v, GEN c, long N)
    2416             : {
    2417       91427 :   GEN w = zerocol(N);
    2418       91427 :   long i, l = lg(c);
    2419      305123 :   for (i=1; i<l; i++) gel(w,c[i]) = gel(v,i);
    2420       91427 :   return w;
    2421             : }
    2422             : 
    2423             : /* return t = 1 mod pr, t = 0 mod p / pr^e(pr/p) */
    2424             : static GEN
    2425      963834 : anti_uniformizer(GEN nf, GEN pr)
    2426             : {
    2427      963834 :   long N = nf_get_degree(nf), e = pr_get_e(pr);
    2428             :   GEN p, b, z;
    2429             : 
    2430      963819 :   if (e * pr_get_f(pr) == N) return gen_1;
    2431      459100 :   p = pr_get_p(pr);
    2432      459098 :   b = pr_get_tau(pr); /* ZM */
    2433      459095 :   if (e != 1)
    2434             :   {
    2435       22638 :     GEN q = powiu(pr_get_p(pr), e-1);
    2436       22638 :     b = ZM_Z_divexact(ZM_powu(b,e), q);
    2437             :   }
    2438             :   /* b = tau^e / p^(e-1), v_pr(b) = 0, v_Q(b) >= e(Q/p) for other Q | p */
    2439      459095 :   z = ZM_hnfmodid(FpM_red(b,p), p); /* ideal (p) / pr^e, coprime to pr */
    2440      459116 :   z = idealaddtoone_raw(nf, pr, z);
    2441      459097 :   return Z_ZC_sub(gen_1, FpC_center(FpC_red(z,p), p, shifti(p,-1)));
    2442             : }
    2443             : 
    2444             : #define mpr_TAU 1
    2445             : #define mpr_FFP 2
    2446             : #define mpr_NFP 5
    2447             : #define SMALLMODPR 4
    2448             : #define LARGEMODPR 6
    2449             : static GEN
    2450     3521645 : modpr_TAU(GEN modpr)
    2451             : {
    2452     3521645 :   GEN tau = gel(modpr,mpr_TAU);
    2453     3521645 :   return isintzero(tau)? NULL: tau;
    2454             : }
    2455             : 
    2456             : /* prh = HNF matrix, which is identity but for the first line. Return a
    2457             :  * projector to ZK / prh ~ Z/prh[1,1] */
    2458             : GEN
    2459     1029298 : dim1proj(GEN prh)
    2460             : {
    2461     1029298 :   long i, N = lg(prh)-1;
    2462     1029298 :   GEN ffproj = cgetg(N+1, t_VEC);
    2463     1029294 :   GEN x, q = gcoeff(prh,1,1);
    2464     1029294 :   gel(ffproj,1) = gen_1;
    2465     2372589 :   for (i=2; i<=N; i++)
    2466             :   {
    2467     1343382 :     x = gcoeff(prh,1,i);
    2468     1343382 :     if (signe(x)) x = subii(q,x);
    2469     1343295 :     gel(ffproj,i) = x;
    2470             :   }
    2471     1029207 :   return ffproj;
    2472             : }
    2473             : 
    2474             : /* p not necessarily prime, but coprime to denom(basis) */
    2475             : GEN
    2476         203 : QXQV_to_FpM(GEN basis, GEN T, GEN p)
    2477             : {
    2478         203 :   long i, l = lg(basis), f = degpol(T);
    2479         203 :   GEN z = cgetg(l, t_MAT);
    2480        4515 :   for (i = 1; i < l; i++)
    2481             :   {
    2482        4312 :     GEN w = gel(basis,i);
    2483        4312 :     if (typ(w) == t_INT)
    2484           0 :       w = scalarcol_shallow(w, f);
    2485             :     else
    2486             :     {
    2487             :       GEN dx;
    2488        4312 :       w = Q_remove_denom(w, &dx);
    2489        4312 :       w = FpXQ_red(w, T, p);
    2490        4312 :       if (dx)
    2491             :       {
    2492           0 :         dx = Fp_inv(dx, p);
    2493           0 :         if (!equali1(dx)) w = FpX_Fp_mul(w, dx, p);
    2494             :       }
    2495        4312 :       w = RgX_to_RgC(w, f);
    2496             :     }
    2497        4312 :     gel(z,i) = w; /* w_i mod (T,p) */
    2498             :   }
    2499         203 :   return z;
    2500             : }
    2501             : 
    2502             : /* initialize reduction mod pr; if zk = 1, will only init data required to
    2503             :  * reduce *integral* element.  Realize (O_K/pr) as Fp[X] / (T), for a
    2504             :  * *monic* T; use variable vT for varn(T) */
    2505             : static GEN
    2506     1157304 : modprinit(GEN nf, GEN pr, int zk, long vT)
    2507             : {
    2508     1157304 :   pari_sp av = avma;
    2509             :   GEN res, tau, mul, x, p, T, pow, ffproj, nfproj, prh, c;
    2510             :   long N, i, k, f;
    2511             : 
    2512     1157304 :   nf = checknf(nf); checkprid(pr);
    2513     1157281 :   if (vT < 0) vT = nf_get_varn(nf);
    2514     1157274 :   f = pr_get_f(pr);
    2515     1157267 :   N = nf_get_degree(nf);
    2516     1157266 :   prh = pr_hnf(nf, pr);
    2517     1157309 :   tau = zk? gen_0: anti_uniformizer(nf, pr);
    2518     1157281 :   p = pr_get_p(pr);
    2519             : 
    2520     1157287 :   if (f == 1)
    2521             :   {
    2522     1011247 :     res = cgetg(SMALLMODPR, t_COL);
    2523     1011239 :     gel(res,mpr_TAU) = tau;
    2524     1011239 :     gel(res,mpr_FFP) = dim1proj(prh);
    2525     1011154 :     gel(res,3) = pr; return gerepilecopy(av, res);
    2526             :   }
    2527             : 
    2528      146040 :   c = cgetg(f+1, t_VECSMALL);
    2529      146041 :   ffproj = cgetg(N+1, t_MAT);
    2530      599802 :   for (k=i=1; i<=N; i++)
    2531             :   {
    2532      453761 :     x = gcoeff(prh, i,i);
    2533      453761 :     if (!is_pm1(x)) { c[k] = i; gel(ffproj,i) = col_ei(N, i); k++; }
    2534             :     else
    2535      127806 :       gel(ffproj,i) = ZC_neg(gel(prh,i));
    2536             :   }
    2537      146041 :   ffproj = rowpermute(ffproj, c);
    2538      146041 :   if (! dvdii(nf_get_index(nf), p))
    2539             :   {
    2540      104524 :     GEN basis = nf_get_zkprimpart(nf), D = nf_get_zkden(nf);
    2541      104524 :     if (N == f)
    2542             :     { /* pr inert */
    2543       45542 :       T = nf_get_pol(nf);
    2544       45542 :       T = FpX_red(T,p);
    2545       45542 :       ffproj = RgV_to_RgM(basis, lg(basis)-1);
    2546             :     }
    2547             :     else
    2548             :     {
    2549       58982 :       T = RgV_RgC_mul(basis, pr_get_gen(pr));
    2550       58982 :       T = FpX_normalize(FpX_red(T,p),p);
    2551       58982 :       basis = FqV_red(vecpermute(basis,c), T, p);
    2552       58982 :       basis = RgV_to_RgM(basis, lg(basis)-1);
    2553       58982 :       ffproj = ZM_mul(basis, ffproj);
    2554             :     }
    2555      104524 :     setvarn(T, vT);
    2556      104524 :     ffproj = FpM_red(ffproj, p);
    2557      104522 :     if (!equali1(D))
    2558             :     {
    2559       33354 :       D = modii(D,p);
    2560       33354 :       if (!equali1(D)) ffproj = FpM_Fp_mul(ffproj, Fp_inv(D,p), p);
    2561             :     }
    2562             : 
    2563      104518 :     res = cgetg(SMALLMODPR+1, t_COL);
    2564      104521 :     gel(res,mpr_TAU) = tau;
    2565      104521 :     gel(res,mpr_FFP) = ffproj;
    2566      104521 :     gel(res,3) = pr;
    2567      104521 :     gel(res,4) = T; return gerepilecopy(av, res);
    2568             :   }
    2569             : 
    2570       41516 :   if (uisprime(f))
    2571             :   {
    2572       39186 :     mul = ei_multable(nf, c[2]);
    2573       39186 :     mul = vecpermute(mul, c);
    2574             :   }
    2575             :   else
    2576             :   {
    2577             :     GEN v, u, u2, frob;
    2578             :     long deg,deg1,deg2;
    2579             : 
    2580             :     /* matrix of Frob: x->x^p over Z_K/pr = < w[c1], ..., w[cf] > over Fp */
    2581        2330 :     frob = cgetg(f+1, t_MAT);
    2582       11963 :     for (i=1; i<=f; i++)
    2583             :     {
    2584        9632 :       x = pow_ei_mod_p(nf,c[i],p);
    2585        9632 :       gel(frob,i) = FpM_FpC_mul(ffproj, x, p);
    2586             :     }
    2587        2331 :     u = col_ei(f,2); k = 2;
    2588        2331 :     deg1 = ffdegree(u, frob, p);
    2589        4102 :     while (deg1 < f)
    2590             :     {
    2591        1771 :       k++; u2 = col_ei(f, k);
    2592        1771 :       deg2 = ffdegree(u2, frob, p);
    2593        1771 :       deg = ulcm(deg1,deg2);
    2594        1771 :       if (deg == deg1) continue;
    2595        1764 :       if (deg == deg2) { deg1 = deg2; u = u2; continue; }
    2596           7 :       u = ZC_add(u, u2);
    2597           7 :       while (ffdegree(u, frob, p) < deg) u = ZC_add(u, u2);
    2598           7 :       deg1 = deg;
    2599             :     }
    2600        2331 :     v = lift_to_zk(u,c,N);
    2601             : 
    2602        2331 :     mul = cgetg(f+1,t_MAT);
    2603        2331 :     gel(mul,1) = v; /* assume w_1 = 1 */
    2604        9632 :     for (i=2; i<=f; i++) gel(mul,i) = zk_ei_mul(nf,v,c[i]);
    2605             :   }
    2606             : 
    2607             :   /* Z_K/pr = Fp(v), mul = mul by v */
    2608       41517 :   mul = FpM_red(mul, p);
    2609       41517 :   mul = FpM_mul(ffproj, mul, p);
    2610             : 
    2611       41517 :   pow = get_powers(mul, p);
    2612       41516 :   T = RgV_to_RgX(FpM_deplin(pow, p), vT);
    2613       41517 :   nfproj = cgetg(f+1, t_MAT);
    2614      130613 :   for (i=1; i<=f; i++) gel(nfproj,i) = lift_to_zk(gel(pow,i), c, N);
    2615             : 
    2616       41517 :   setlg(pow, f+1);
    2617       41517 :   ffproj = FpM_mul(FpM_inv(pow, p), ffproj, p);
    2618             : 
    2619       41517 :   res = cgetg(LARGEMODPR, t_COL);
    2620       41517 :   gel(res,mpr_TAU) = tau;
    2621       41517 :   gel(res,mpr_FFP) = ffproj;
    2622       41517 :   gel(res,3) = pr;
    2623       41517 :   gel(res,4) = T;
    2624       41517 :   gel(res,mpr_NFP) = nfproj; return gerepilecopy(av, res);
    2625             : }
    2626             : 
    2627             : GEN
    2628           7 : nfmodprinit(GEN nf, GEN pr) { return modprinit(nf, pr, 0, -1); }
    2629             : GEN
    2630      174935 : zkmodprinit(GEN nf, GEN pr) { return modprinit(nf, pr, 1, -1); }
    2631             : GEN
    2632          70 : nfmodprinit0(GEN nf, GEN pr, long v) { return modprinit(nf, pr, 0, v); }
    2633             : 
    2634             : /* x may be a modpr */
    2635             : static int
    2636     4341211 : ok_modpr(GEN x)
    2637     4341211 : { return typ(x) == t_COL && lg(x) >= SMALLMODPR && lg(x) <= LARGEMODPR; }
    2638             : void
    2639         210 : checkmodpr(GEN x)
    2640             : {
    2641         210 :   if (!ok_modpr(x)) pari_err_TYPE("checkmodpr [use nfmodprinit]", x);
    2642         210 :   checkprid(modpr_get_pr(x));
    2643         210 : }
    2644             : GEN
    2645      137746 : get_modpr(GEN x)
    2646      137746 : { return ok_modpr(x)? x: NULL; }
    2647             : 
    2648             : int
    2649    12233130 : checkprid_i(GEN x)
    2650             : {
    2651    11435257 :   return (typ(x) == t_VEC && lg(x) == 6
    2652    11362283 :           && typ(gel(x,2)) == t_COL && typ(gel(x,3)) == t_INT
    2653    23668387 :           && typ(gel(x,5)) != t_COL); /* tau changed to t_MAT/t_INT in 2.6 */
    2654             : }
    2655             : void
    2656    10668049 : checkprid(GEN x)
    2657    10668049 : { if (!checkprid_i(x)) pari_err_TYPE("checkprid",x); }
    2658             : GEN
    2659      918232 : get_prid(GEN x)
    2660             : {
    2661      918232 :   long lx = lg(x);
    2662      918232 :   if (lx == 3 && typ(x) == t_VEC) x = gel(x,1);
    2663      918232 :   if (checkprid_i(x)) return x;
    2664      676004 :   if (ok_modpr(x)) {
    2665      108465 :     x = modpr_get_pr(x);
    2666      108465 :     if (checkprid_i(x)) return x;
    2667             :   }
    2668      567539 :   return NULL;
    2669             : }
    2670             : 
    2671             : static GEN
    2672     3527252 : to_ff_init(GEN nf, GEN *pr, GEN *T, GEN *p, int zk)
    2673             : {
    2674     3527252 :   GEN modpr = ok_modpr(*pr)? *pr: modprinit(nf, *pr, zk, -1);
    2675     3527360 :   *T = modpr_get_T(modpr);
    2676     3527302 :   *pr = modpr_get_pr(modpr);
    2677     3527281 :   *p = pr_get_p(*pr); return modpr;
    2678             : }
    2679             : 
    2680             : /* Return an element of O_K which is set to x Mod T */
    2681             : GEN
    2682        4326 : modpr_genFq(GEN modpr)
    2683             : {
    2684        4326 :   switch(lg(modpr))
    2685             :   {
    2686         917 :     case SMALLMODPR: /* Fp */
    2687         917 :       return gen_1;
    2688        1554 :     case LARGEMODPR:  /* painful case, p \mid index */
    2689        1554 :       return gmael(modpr,mpr_NFP, 2);
    2690        1855 :     default: /* trivial case : p \nmid index */
    2691             :     {
    2692        1855 :       long v = varn( modpr_get_T(modpr) );
    2693        1855 :       return pol_x(v);
    2694             :     }
    2695             :   }
    2696             : }
    2697             : 
    2698             : GEN
    2699     3508694 : nf_to_Fq_init(GEN nf, GEN *pr, GEN *T, GEN *p) {
    2700     3508694 :   GEN modpr = to_ff_init(nf,pr,T,p,0);
    2701     3508747 :   GEN tau = modpr_TAU(modpr);
    2702     3508721 :   if (!tau) gel(modpr,mpr_TAU) = anti_uniformizer(nf, *pr);
    2703     3508721 :   return modpr;
    2704             : }
    2705             : GEN
    2706       18536 : zk_to_Fq_init(GEN nf, GEN *pr, GEN *T, GEN *p) {
    2707       18536 :   return to_ff_init(nf,pr,T,p,1);
    2708             : }
    2709             : 
    2710             : /* assume x in 'basis' form (t_COL) */
    2711             : GEN
    2712     7098035 : zk_to_Fq(GEN x, GEN modpr)
    2713             : {
    2714     7098035 :   GEN pr = modpr_get_pr(modpr), p = pr_get_p(pr);
    2715     7098018 :   GEN ffproj = gel(modpr,mpr_FFP);
    2716     7098018 :   GEN T = modpr_get_T(modpr);
    2717     7098034 :   return T? FpM_FpC_mul_FpX(ffproj,x, p, varn(T)): FpV_dotproduct(ffproj,x, p);
    2718             : }
    2719             : 
    2720             : /* REDUCTION Modulo a prime ideal */
    2721             : 
    2722             : /* nf a true nf */
    2723             : static GEN
    2724    15137667 : Rg_to_ff(GEN nf, GEN x0, GEN modpr)
    2725             : {
    2726    15137667 :   GEN x = x0, den, pr = modpr_get_pr(modpr), p = pr_get_p(pr);
    2727    15137665 :   long tx = typ(x);
    2728             : 
    2729    15137665 :   if (tx == t_POLMOD) { x = gel(x,2); tx = typ(x); }
    2730    15137665 :   switch(tx)
    2731             :   {
    2732     8121155 :     case t_INT: return modii(x, p);
    2733        8465 :     case t_FRAC: return Rg_to_Fp(x, p);
    2734      207954 :     case t_POL:
    2735      207954 :       switch(lg(x))
    2736             :       {
    2737         224 :         case 2: return gen_0;
    2738       25782 :         case 3: return Rg_to_Fp(gel(x,2), p);
    2739             :       }
    2740      181948 :       x = Q_remove_denom(x, &den);
    2741      181948 :       x = poltobasis(nf, x);
    2742             :       /* content(x) and den may not be coprime */
    2743      181777 :       break;
    2744     6800126 :     case t_COL:
    2745     6800126 :       x = Q_remove_denom(x, &den);
    2746             :       /* content(x) and den are coprime */
    2747     6800125 :       if (lg(x)-1 == nf_get_degree(nf)) break;
    2748          48 :     default: pari_err_TYPE("Rg_to_ff",x);
    2749             :       return NULL;/*LCOV_EXCL_LINE*/
    2750             :   }
    2751     6981844 :   if (den)
    2752             :   {
    2753      134966 :     long v = Z_pvalrem(den, p, &den);
    2754      134966 :     if (v)
    2755             :     {
    2756       13475 :       if (tx == t_POL) v -= ZV_pvalrem(x, p, &x);
    2757             :       /* now v = valuation(true denominator of x) */
    2758       13475 :       if (v > 0)
    2759             :       {
    2760       12901 :         GEN tau = modpr_TAU(modpr);
    2761       12901 :         if (!tau) pari_err_TYPE("zk_to_ff", x0);
    2762       12901 :         x = nfmuli(nf,x, nfpow_u(nf, tau, v));
    2763       12901 :         v -= ZV_pvalrem(x, p, &x);
    2764             :       }
    2765       13475 :       if (v > 0) pari_err_INV("Rg_to_ff", mkintmod(gen_0,p));
    2766       13447 :       if (v) return gen_0;
    2767        8638 :       if (is_pm1(den)) den = NULL;
    2768             :     }
    2769      130129 :     x = FpC_red(x, p);
    2770             :   }
    2771     6977007 :   x = zk_to_Fq(x, modpr);
    2772     6977008 :   if (den)
    2773             :   {
    2774      124796 :     GEN c = Fp_inv(den, p);
    2775      124809 :     x = typ(x) == t_INT? Fp_mul(x,c,p): FpX_Fp_mul(x,c,p);
    2776             :   }
    2777     6977021 :   return x;
    2778             : }
    2779             : 
    2780             : GEN
    2781         210 : nfreducemodpr(GEN nf, GEN x, GEN modpr)
    2782             : {
    2783         210 :   pari_sp av = avma;
    2784         210 :   nf = checknf(nf); checkmodpr(modpr);
    2785         210 :   return gerepileupto(av, algtobasis(nf, Fq_to_nf(Rg_to_ff(nf,x,modpr),modpr)));
    2786             : }
    2787             : 
    2788             : GEN
    2789         329 : nfmodpr(GEN nf, GEN x, GEN pr)
    2790             : {
    2791         329 :   pari_sp av = avma;
    2792             :   GEN T, p, modpr;
    2793         329 :   nf = checknf(nf);
    2794         329 :   modpr = nf_to_Fq_init(nf, &pr, &T, &p);
    2795         322 :   if (typ(x) == t_MAT && lg(x) == 3)
    2796             :   {
    2797          35 :     GEN y, v = famat_nfvalrem(nf, x, pr, &y);
    2798          35 :     long s = signe(v);
    2799          35 :     if (s < 0) pari_err_INV("Rg_to_ff", mkintmod(gen_0,p));
    2800          28 :     if (s > 0) return gc_const(av, gen_0);
    2801          14 :     x = FqV_factorback(nfV_to_FqV(gel(y,1), nf, modpr), gel(y,2), T, p);
    2802          14 :     return gerepileupto(av, x);
    2803             :   }
    2804         287 :   x = Rg_to_ff(nf, x, modpr);
    2805         175 :   x = Fq_to_FF(x, Tp_to_FF(T,p));
    2806         175 :   return gerepilecopy(av, x);
    2807             : }
    2808             : GEN
    2809          77 : nfmodprlift(GEN nf, GEN x, GEN pr)
    2810             : {
    2811          77 :   pari_sp av = avma;
    2812             :   GEN y, T, p, modpr;
    2813             :   long i, l, d;
    2814          77 :   nf = checknf(nf);
    2815          77 :   switch(typ(x))
    2816             :   {
    2817           7 :     case t_INT: return icopy(x);
    2818          42 :     case t_FFELT: break;
    2819          28 :     case t_VEC: case t_COL: case t_MAT:
    2820          28 :       y = cgetg_copy(x,&l);
    2821          63 :       for (i = 1; i < l; i++) gel(y,i) = nfmodprlift(nf,gel(x,i),pr);
    2822          28 :       return y;
    2823           0 :     default: pari_err_TYPE("nfmodprlit",x);
    2824             :   }
    2825          42 :   x = FF_to_FpXQ(x);
    2826          42 :   setvarn(x, nf_get_varn(nf));
    2827          42 :   d = degpol(x);
    2828          42 :   if (d <= 0) { set_avma(av); return d? gen_0: icopy(gel(x,2)); }
    2829          14 :   modpr = nf_to_Fq_init(nf, &pr, &T, &p);
    2830          14 :   return gerepilecopy(av, Fq_to_nf(x, modpr));
    2831             : }
    2832             : 
    2833             : /* lift A from residue field to nf */
    2834             : GEN
    2835     3229346 : Fq_to_nf(GEN A, GEN modpr)
    2836             : {
    2837             :   long dA;
    2838     3229346 :   if (typ(A) == t_INT || lg(modpr) < LARGEMODPR) return A;
    2839       44375 :   dA = degpol(A);
    2840       44375 :   if (dA <= 0) return dA ? gen_0: gel(A,2);
    2841       40335 :   return ZM_ZX_mul(gel(modpr,mpr_NFP), A);
    2842             : }
    2843             : GEN
    2844           0 : FqV_to_nfV(GEN x, GEN modpr)
    2845           0 : { pari_APPLY_same(Fq_to_nf(gel(x,i), modpr)) }
    2846             : GEN
    2847        9924 : FqM_to_nfM(GEN A, GEN modpr)
    2848             : {
    2849        9924 :   long i,j,h,l = lg(A);
    2850        9924 :   GEN B = cgetg(l, t_MAT);
    2851             : 
    2852        9924 :   if (l == 1) return B;
    2853        9232 :   h = lgcols(A);
    2854       41731 :   for (j=1; j<l; j++)
    2855             :   {
    2856       32499 :     GEN Aj = gel(A,j), Bj = cgetg(h,t_COL); gel(B,j) = Bj;
    2857      214229 :     for (i=1; i<h; i++) gel(Bj,i) = Fq_to_nf(gel(Aj,i), modpr);
    2858             :   }
    2859        9232 :   return B;
    2860             : }
    2861             : GEN
    2862       11549 : FqX_to_nfX(GEN A, GEN modpr)
    2863             : {
    2864             :   long i, l;
    2865             :   GEN B;
    2866             : 
    2867       11549 :   if (typ(A)!=t_POL) return icopy(A); /* scalar */
    2868       11549 :   B = cgetg_copy(A, &l); B[1] = A[1];
    2869       51321 :   for (i=2; i<l; i++) gel(B,i) = Fq_to_nf(gel(A,i), modpr);
    2870       11549 :   return B;
    2871             : }
    2872             : 
    2873             : /* reduce A to residue field */
    2874             : GEN
    2875    15137158 : nf_to_Fq(GEN nf, GEN A, GEN modpr)
    2876             : {
    2877    15137158 :   pari_sp av = avma;
    2878    15137158 :   return gerepileupto(av, Rg_to_ff(checknf(nf), A, modpr));
    2879             : }
    2880             : /* A t_VEC/t_COL */
    2881             : GEN
    2882      167403 : nfV_to_FqV(GEN A, GEN nf,GEN modpr)
    2883             : {
    2884      167403 :   long i,l = lg(A);
    2885      167403 :   GEN B = cgetg(l,typ(A));
    2886      916739 :   for (i=1; i<l; i++) gel(B,i) = nf_to_Fq(nf,gel(A,i), modpr);
    2887      167403 :   return B;
    2888             : }
    2889             : /* A  t_MAT */
    2890             : GEN
    2891        5382 : nfM_to_FqM(GEN A, GEN nf,GEN modpr)
    2892             : {
    2893        5382 :   long i,j,h,l = lg(A);
    2894        5382 :   GEN B = cgetg(l,t_MAT);
    2895             : 
    2896        5382 :   if (l == 1) return B;
    2897        5382 :   h = lgcols(A);
    2898      153344 :   for (j=1; j<l; j++)
    2899             :   {
    2900      147962 :     GEN Aj = gel(A,j), Bj = cgetg(h,t_COL); gel(B,j) = Bj;
    2901     1044416 :     for (i=1; i<h; i++) gel(Bj,i) = nf_to_Fq(nf, gel(Aj,i), modpr);
    2902             :   }
    2903        5382 :   return B;
    2904             : }
    2905             : /* A t_POL */
    2906             : GEN
    2907        9751 : nfX_to_FqX(GEN A, GEN nf,GEN modpr)
    2908             : {
    2909        9751 :   long i,l = lg(A);
    2910        9751 :   GEN B = cgetg(l,t_POL); B[1] = A[1];
    2911       54620 :   for (i=2; i<l; i++) gel(B,i) = nf_to_Fq(nf,gel(A,i),modpr);
    2912        9744 :   return normalizepol_lg(B, l);
    2913             : }
    2914             : 
    2915             : /*******************************************************************/
    2916             : /*                                                                 */
    2917             : /*                       RELATIVE ROUND 2                          */
    2918             : /*                                                                 */
    2919             : /*******************************************************************/
    2920             : /* Shallow functions */
    2921             : /* FIXME: use a bb_field and export the nfX_* routines */
    2922             : static GEN
    2923        4403 : nfX_sub(GEN nf, GEN x, GEN y)
    2924             : {
    2925        4403 :   long i, lx = lg(x), ly = lg(y);
    2926             :   GEN z;
    2927        4403 :   if (ly <= lx) {
    2928        4403 :     z = cgetg(lx,t_POL); z[1] = x[1];
    2929       27237 :     for (i=2; i < ly; i++) gel(z,i) = nfsub(nf,gel(x,i),gel(y,i));
    2930        4403 :     for (   ; i < lx; i++) gel(z,i) = gel(x,i);
    2931        4403 :     z = normalizepol_lg(z, lx);
    2932             :   } else {
    2933           0 :     z = cgetg(ly,t_POL); z[1] = y[1];
    2934           0 :     for (i=2; i < lx; i++) gel(z,i) = nfsub(nf,gel(x,i),gel(y,i));
    2935           0 :     for (   ; i < ly; i++) gel(z,i) = gneg(gel(y,i));
    2936           0 :     z = normalizepol_lg(z, ly);
    2937             :   }
    2938        4403 :   return z;
    2939             : }
    2940             : /* FIXME: quadratic multiplication */
    2941             : static GEN
    2942       65148 : nfX_mul(GEN nf, GEN a, GEN b)
    2943             : {
    2944       65148 :   long da = degpol(a), db = degpol(b), dc, lc, k;
    2945             :   GEN c;
    2946       65148 :   if (da < 0 || db < 0) return gen_0;
    2947       65148 :   dc = da + db;
    2948       65148 :   if (dc == 0) return nfmul(nf, gel(a,2),gel(b,2));
    2949       65148 :   lc = dc+3;
    2950       65148 :   c = cgetg(lc, t_POL); c[1] = a[1];
    2951      522973 :   for (k = 0; k <= dc; k++)
    2952             :   {
    2953      457825 :     long i, I = minss(k, da);
    2954      457825 :     GEN d = NULL;
    2955     1559201 :     for (i = maxss(k-db, 0); i <= I; i++)
    2956             :     {
    2957     1101376 :       GEN e = nfmul(nf, gel(a, i+2), gel(b, k-i+2));
    2958     1101376 :       d = d? nfadd(nf, d, e): e;
    2959             :     }
    2960      457825 :     gel(c, k+2) = d;
    2961             :   }
    2962       65148 :   return normalizepol_lg(c, lc);
    2963             : }
    2964             : /* assume b monic */
    2965             : static GEN
    2966       60745 : nfX_rem(GEN nf, GEN a, GEN b)
    2967             : {
    2968       60745 :   long da = degpol(a), db = degpol(b);
    2969       60745 :   if (da < 0) return gen_0;
    2970       60745 :   a = leafcopy(a);
    2971      145801 :   while (da >= db)
    2972             :   {
    2973       85056 :     long i, k = da;
    2974       85056 :     GEN A = gel(a, k+2);
    2975      620631 :     for (i = db-1, k--; i >= 0; i--, k--)
    2976      535575 :       gel(a,k+2) = nfsub(nf, gel(a,k+2), nfmul(nf, A, gel(b,i+2)));
    2977       85056 :     a = normalizepol_lg(a, lg(a)-1);
    2978       85056 :     da = degpol(a);
    2979             :   }
    2980       60745 :   return a;
    2981             : }
    2982             : static GEN
    2983       60745 : nfXQ_mul(GEN nf, GEN a, GEN b, GEN T)
    2984             : {
    2985       60745 :   GEN c = nfX_mul(nf, a, b);
    2986       60745 :   if (typ(c) != t_POL) return c;
    2987       60745 :   return nfX_rem(nf, c, T);
    2988             : }
    2989             : 
    2990             : static void
    2991       12584 : fill(long l, GEN H, GEN Hx, GEN I, GEN Ix)
    2992             : {
    2993             :   long i;
    2994       12584 :   if (typ(Ix) == t_VEC) /* standard */
    2995       46428 :     for (i=1; i<l; i++) { gel(H,i) = gel(Hx,i); gel(I,i) = gel(Ix,i); }
    2996             :   else /* constant ideal */
    2997       12739 :     for (i=1; i<l; i++) { gel(H,i) = gel(Hx,i); gel(I,i) = Ix; }
    2998       12584 : }
    2999             : 
    3000             : /* given MODULES x and y by their pseudo-bases, returns a pseudo-basis of the
    3001             :  * module generated by x and y. */
    3002             : static GEN
    3003        6292 : rnfjoinmodules_i(GEN nf, GEN Hx, GEN Ix, GEN Hy, GEN Iy)
    3004             : {
    3005        6292 :   long lx = lg(Hx), ly = lg(Hy), l = lx+ly-1;
    3006        6292 :   GEN H = cgetg(l, t_MAT), I = cgetg(l, t_VEC);
    3007        6292 :   fill(lx, H     , Hx, I     , Ix);
    3008        6292 :   fill(ly, H+lx-1, Hy, I+lx-1, Iy); return nfhnf(nf, mkvec2(H, I));
    3009             : }
    3010             : static GEN
    3011        2134 : rnfjoinmodules(GEN nf, GEN x, GEN y)
    3012             : {
    3013        2134 :   if (!x) return y;
    3014        1330 :   if (!y) return x;
    3015        1330 :   return rnfjoinmodules_i(nf, gel(x,1), gel(x,2), gel(y,1), gel(y,2));
    3016             : }
    3017             : 
    3018             : typedef struct {
    3019             :   GEN multab, T,p;
    3020             :   long h;
    3021             : } rnfeltmod_muldata;
    3022             : 
    3023             : static GEN
    3024       69641 : _sqr(void *data, GEN x)
    3025             : {
    3026       69641 :   rnfeltmod_muldata *D = (rnfeltmod_muldata *) data;
    3027       48551 :   GEN z = x? tablesqr(D->multab,x)
    3028       69641 :            : tablemul_ei_ej(D->multab,D->h,D->h);
    3029       69641 :   return FqV_red(z,D->T,D->p);
    3030             : }
    3031             : static GEN
    3032       11347 : _msqr(void *data, GEN x)
    3033             : {
    3034       11347 :   GEN x2 = _sqr(data, x), z;
    3035       11347 :   rnfeltmod_muldata *D = (rnfeltmod_muldata *) data;
    3036       11347 :   z = tablemul_ei(D->multab, x2, D->h);
    3037       11347 :   return FqV_red(z,D->T,D->p);
    3038             : }
    3039             : 
    3040             : /* Compute W[h]^n mod (T,p) in the extension, assume n >= 0. T a ZX */
    3041             : static GEN
    3042       21090 : rnfeltid_powmod(GEN multab, long h, GEN n, GEN T, GEN p)
    3043             : {
    3044       21090 :   pari_sp av = avma;
    3045             :   GEN y;
    3046             :   rnfeltmod_muldata D;
    3047             : 
    3048       21090 :   if (!signe(n)) return gen_1;
    3049             : 
    3050       21090 :   D.multab = multab;
    3051       21090 :   D.h = h;
    3052       21090 :   D.T = T;
    3053       21090 :   D.p = p;
    3054       21090 :   y = gen_pow_fold(NULL, n, (void*)&D, &_sqr, &_msqr);
    3055       21090 :   return gerepilecopy(av, y);
    3056             : }
    3057             : 
    3058             : /* P != 0 has at most degpol(P) roots. Look for an element in Fq which is not
    3059             :  * a root, cf repres() */
    3060             : static GEN
    3061          21 : FqX_non_root(GEN P, GEN T, GEN p)
    3062             : {
    3063          21 :   long dP = degpol(P), f, vT;
    3064             :   long i, j, k, pi, pp;
    3065             :   GEN v;
    3066             : 
    3067          21 :   if (dP == 0) return gen_1;
    3068          21 :   pp = is_bigint(p) ? dP+1: itos(p);
    3069          21 :   v = cgetg(dP + 2, t_VEC);
    3070          21 :   gel(v,1) = gen_0;
    3071          21 :   if (T)
    3072           0 :   { f = degpol(T); vT = varn(T); }
    3073             :   else
    3074          21 :   { f = 1; vT = 0; }
    3075          42 :   for (i=pi=1; i<=f; i++,pi*=pp)
    3076             :   {
    3077          21 :     GEN gi = i == 1? gen_1: pol_xn(i-1, vT), jgi = gi;
    3078          42 :     for (j=1; j<pp; j++)
    3079             :     {
    3080          42 :       for (k=1; k<=pi; k++)
    3081             :       {
    3082          21 :         GEN z = Fq_add(gel(v,k), jgi, T,p);
    3083          21 :         if (!gequal0(FqX_eval(P, z, T,p))) return z;
    3084          21 :         gel(v, j*pi+k) = z;
    3085             :       }
    3086          21 :       if (j < pp-1) jgi = Fq_add(jgi, gi, T,p); /* j*g[i] */
    3087             :     }
    3088             :   }
    3089          21 :   return NULL;
    3090             : }
    3091             : 
    3092             : /* Relative Dedekind criterion over (true) nf, applied to the order defined by a
    3093             :  * root of monic irreducible polynomial P, modulo the prime ideal pr. Assume
    3094             :  * vdisc = v_pr( disc(P) ).
    3095             :  * Return NULL if nf[X]/P is pr-maximal. Otherwise, return [flag, O, v]:
    3096             :  *   O = enlarged order, given by a pseudo-basis
    3097             :  *   flag = 1 if O is proven pr-maximal (may be 0 and O nevertheless pr-maximal)
    3098             :  *   v = v_pr(disc(O)). */
    3099             : static GEN
    3100        4438 : rnfdedekind_i(GEN nf, GEN P, GEN pr, long vdisc, long only_maximal)
    3101             : {
    3102             :   GEN Ppr, A, I, p, tau, g, h, k, base, T, gzk, hzk, prinvp, pal, nfT, modpr;
    3103             :   long m, vt, r, d, i, j, mpr;
    3104             : 
    3105        4438 :   if (vdisc < 0) pari_err_TYPE("rnfdedekind [non integral pol]", P);
    3106        4431 :   if (vdisc == 1) return NULL; /* pr-maximal */
    3107        4431 :   if (!only_maximal && !gequal1(leading_coeff(P)))
    3108           0 :     pari_err_IMPL( "the full Dedekind criterion in the nonmonic case");
    3109             :   /* either monic OR only_maximal = 1 */
    3110        4431 :   m = degpol(P);
    3111        4431 :   nfT = nf_get_pol(nf);
    3112        4431 :   modpr = nf_to_Fq_init(nf,&pr, &T, &p);
    3113        4431 :   Ppr = nfX_to_FqX(P, nf, modpr);
    3114        4424 :   mpr = degpol(Ppr);
    3115        4424 :   if (mpr < m) /* nonmonic => only_maximal = 1 */
    3116             :   {
    3117          21 :     if (mpr < 0) return NULL;
    3118          21 :     if (! RgX_valrem(Ppr, &Ppr))
    3119             :     { /* nonzero constant coefficient */
    3120           0 :       Ppr = RgX_shift_shallow(RgX_recip_i(Ppr), m - mpr);
    3121           0 :       P = RgX_recip_i(P);
    3122             :     }
    3123             :     else
    3124             :     {
    3125          21 :       GEN z = FqX_non_root(Ppr, T, p);
    3126          21 :       if (!z) pari_err_IMPL( "Dedekind in the difficult case");
    3127           0 :       z = Fq_to_nf(z, modpr);
    3128           0 :       if (typ(z) == t_INT)
    3129           0 :         P = RgX_translate(P, z);
    3130             :       else
    3131           0 :         P = RgXQX_translate(P, z, T);
    3132           0 :       P = RgX_recip_i(P);
    3133           0 :       Ppr = nfX_to_FqX(P, nf, modpr); /* degpol(P) = degpol(Ppr) = m */
    3134             :     }
    3135             :   }
    3136        4403 :   A = gel(FqX_factor(Ppr,T,p),1);
    3137        4403 :   r = lg(A); /* > 1 */
    3138        4403 :   g = gel(A,1);
    3139        8813 :   for (i=2; i<r; i++) g = FqX_mul(g, gel(A,i), T, p);
    3140        4403 :   h = FqX_div(Ppr,g, T, p);
    3141        4403 :   gzk = FqX_to_nfX(g, modpr);
    3142        4403 :   hzk = FqX_to_nfX(h, modpr);
    3143        4403 :   k = nfX_sub(nf, P, nfX_mul(nf, gzk,hzk));
    3144        4403 :   tau = pr_get_tau(pr);
    3145        4403 :   switch(typ(tau))
    3146             :   {
    3147        1883 :     case t_INT: k = gdiv(k, p); break;
    3148        2520 :     case t_MAT: k = RgX_Rg_div(tablemulvec(NULL,tau, k), p); break;
    3149             :   }
    3150        4403 :   k = nfX_to_FqX(k, nf, modpr);
    3151        4403 :   k = FqX_normalize(FqX_gcd(FqX_gcd(g,h,  T,p), k, T,p), T,p);
    3152        4403 :   d = degpol(k);  /* <= m */
    3153        4403 :   if (!d) return NULL; /* pr-maximal */
    3154        2757 :   if (only_maximal) return gen_0; /* not maximal */
    3155             : 
    3156        2736 :   A = cgetg(m+d+1,t_MAT);
    3157        2736 :   I = cgetg(m+d+1,t_VEC); base = mkvec2(A, I);
    3158             :  /* base[2] temporarily multiplied by p, for the final nfhnfmod,
    3159             :   * which requires integral ideals */
    3160        2736 :   prinvp = pr_inv_p(pr); /* again multiplied by p */
    3161       15642 :   for (j=1; j<=m; j++)
    3162             :   {
    3163       12906 :     gel(A,j) = col_ei(m, j);
    3164       12906 :     gel(I,j) = p;
    3165             :   }
    3166        2736 :   pal = FqX_to_nfX(FqX_div(Ppr,k, T,p), modpr);
    3167        5857 :   for (   ; j<=m+d; j++)
    3168             :   {
    3169        3121 :     gel(A,j) = RgX_to_RgC(pal,m);
    3170        3121 :     gel(I,j) = prinvp;
    3171        3121 :     if (j < m+d) pal = RgXQX_rem(RgX_shift_shallow(pal,1),P,nfT);
    3172             :   }
    3173             :   /* the modulus is integral */
    3174        2736 :   base = nfhnfmod(nf,base, idealmulpowprime(nf, powiu(p,m), pr, utoineg(d)));
    3175        2736 :   gel(base,2) = gdiv(gel(base,2), p); /* cancel the factor p */
    3176        2736 :   vt = vdisc - 2*d;
    3177        2736 :   return mkvec3(vt < 2? gen_1: gen_0, base, stoi(vt));
    3178             : }
    3179             : 
    3180             : /* [L:K] = n */
    3181             : static GEN
    3182        1016 : triv_order(long n)
    3183             : {
    3184        1016 :   GEN z = cgetg(3, t_VEC);
    3185        1016 :   gel(z,1) = matid(n);
    3186        1016 :   gel(z,2) = const_vec(n, gen_1); return z;
    3187             : }
    3188             : 
    3189             : /* if flag is set, return gen_1 (resp. gen_0) if the order K[X]/(P)
    3190             :  * is pr-maximal (resp. not pr-maximal). */
    3191             : GEN
    3192          91 : rnfdedekind(GEN nf, GEN P, GEN pr, long flag)
    3193             : {
    3194          91 :   pari_sp av = avma;
    3195             :   GEN z, dP;
    3196             :   long v;
    3197             : 
    3198          91 :   nf = checknf(nf);
    3199          91 :   P = RgX_nffix("rnfdedekind", nf_get_pol(nf), P, 1);
    3200          91 :   dP = nfX_disc(nf, P);
    3201          91 :   if (gequal0(dP))
    3202           7 :     pari_err_DOMAIN("rnfdedekind","issquarefree(pol)","=",gen_0,P);
    3203          84 :   if (!pr)
    3204             :   {
    3205          21 :     GEN fa = idealfactor(nf, dP);
    3206          21 :     GEN Q = gel(fa,1), E = gel(fa,2);
    3207          21 :     pari_sp av2 = avma;
    3208          21 :     long i, l = lg(Q);
    3209          21 :     for (i = 1; i < l; i++, set_avma(av2))
    3210             :     {
    3211          21 :       v = itos(gel(E,i));
    3212          21 :       if (rnfdedekind_i(nf,P,gel(Q,i),v,1)) { set_avma(av); return gen_0; }
    3213           0 :       set_avma(av2);
    3214             :     }
    3215           0 :     set_avma(av); return gen_1;
    3216             :   }
    3217          63 :   else if (typ(pr) == t_VEC)
    3218             :   { /* flag = 1 is implicit */
    3219          63 :     if (lg(pr) == 1) { set_avma(av); return gen_1; }
    3220          63 :     if (typ(gel(pr,1)) == t_VEC)
    3221             :     { /* list of primes */
    3222          14 :       GEN Q = pr;
    3223          14 :       pari_sp av2 = avma;
    3224          14 :       long i, l = lg(Q);
    3225          14 :       for (i = 1; i < l; i++, set_avma(av2))
    3226             :       {
    3227          14 :         v = nfval(nf, dP, gel(Q,i));
    3228          14 :         if (rnfdedekind_i(nf,P,gel(Q,i),v,1)) { set_avma(av); return gen_0; }
    3229             :       }
    3230           0 :       set_avma(av); return gen_1;
    3231             :     }
    3232             :   }
    3233             :   /* single prime */
    3234          49 :   v = nfval(nf, dP, pr);
    3235          49 :   z = rnfdedekind_i(nf, P, pr, v, flag);
    3236          42 :   if (z)
    3237             :   {
    3238          21 :     if (flag) { set_avma(av); return gen_0; }
    3239          14 :     z = gerepilecopy(av, z);
    3240             :   }
    3241             :   else
    3242             :   {
    3243          21 :     set_avma(av); if (flag) return gen_1;
    3244           7 :     z = cgetg(4, t_VEC);
    3245           7 :     gel(z,1) = gen_1;
    3246           7 :     gel(z,2) = triv_order(degpol(P));
    3247           7 :     gel(z,3) = stoi(v);
    3248             :   }
    3249          21 :   return z;
    3250             : }
    3251             : 
    3252             : static int
    3253       26815 : ideal_is1(GEN x) {
    3254       26815 :   switch(typ(x))
    3255             :   {
    3256       11542 :     case t_INT: return is_pm1(x);
    3257       14433 :     case t_MAT: return RgM_isidentity(x);
    3258             :   }
    3259         840 :   return 0;
    3260             : }
    3261             : 
    3262             : /* return a in ideal A such that v_pr(a) = v_pr(A) */
    3263             : static GEN
    3264       15119 : minval(GEN nf, GEN A, GEN pr)
    3265             : {
    3266       15119 :   GEN ab = idealtwoelt(nf,A), a = gel(ab,1), b = gel(ab,2);
    3267       15119 :   if (nfval(nf,a,pr) > nfval(nf,b,pr)) a = b;
    3268       15119 :   return a;
    3269             : }
    3270             : 
    3271             : /* nf a true nf. Return NULL if power order is pr-maximal */
    3272             : static GEN
    3273        4354 : rnfmaxord(GEN nf, GEN pol, GEN pr, long vdisc)
    3274             : {
    3275        4354 :   pari_sp av = avma, av1;
    3276             :   long i, j, k, n, nn, vpol, cnt, sep;
    3277             :   GEN q, q1, p, T, modpr, W, I, p1;
    3278             :   GEN prhinv, mpi, Id;
    3279             : 
    3280        4354 :   if (DEBUGLEVEL>1) err_printf(" treating %Ps^%ld\n", pr, vdisc);
    3281        4354 :   modpr = nf_to_Fq_init(nf,&pr,&T,&p);
    3282        4354 :   av1 = avma;
    3283        4354 :   p1 = rnfdedekind_i(nf, pol, modpr, vdisc, 0);
    3284        4347 :   if (!p1) return gc_NULL(av);
    3285        2722 :   if (is_pm1(gel(p1,1))) return gerepilecopy(av,gel(p1,2));
    3286        1238 :   sep = itos(gel(p1,3));
    3287        1238 :   W = gmael(p1,2,1);
    3288        1238 :   I = gmael(p1,2,2);
    3289        1238 :   gerepileall(av1, 2, &W, &I);
    3290             : 
    3291        1238 :   mpi = zk_multable(nf, pr_get_gen(pr));
    3292        1238 :   n = degpol(pol); nn = n*n;
    3293        1238 :   vpol = varn(pol);
    3294        1238 :   q1 = q = pr_norm(pr);
    3295        1721 :   while (abscmpiu(q1,n) < 0) q1 = mulii(q1,q);
    3296        1238 :   Id = matid(n);
    3297        1238 :   prhinv = pr_inv(pr);
    3298        1238 :   av1 = avma;
    3299        1238 :   for(cnt=1;; cnt++)
    3300        4144 :   {
    3301        5382 :     GEN I0 = leafcopy(I), W0 = leafcopy(W);
    3302             :     GEN Wa, Winv, Ip, A, MW, MWmod, F, pseudo, C, G;
    3303        5382 :     GEN Tauinv = cgetg(n+1, t_VEC), Tau = cgetg(n+1, t_VEC);
    3304             : 
    3305        5382 :     if (DEBUGLEVEL>1) err_printf("    pass no %ld\n",cnt);
    3306       31854 :     for (j=1; j<=n; j++)
    3307             :     {
    3308             :       GEN tau, tauinv;
    3309       26472 :       if (ideal_is1(gel(I,j)))
    3310             :       {
    3311       11353 :         gel(I,j) = gel(Tau,j) = gel(Tauinv,j) = gen_1;
    3312       11353 :         continue;
    3313             :       }
    3314       15119 :       gel(Tau,j) = tau = minval(nf, gel(I,j), pr);
    3315       15119 :       gel(Tauinv,j) = tauinv = nfinv(nf, tau);
    3316       15119 :       gel(W,j) = nfC_nf_mul(nf, gel(W,j), tau);
    3317       15119 :       gel(I,j) = idealmul(nf, tauinv, gel(I,j)); /* v_pr(I[j]) = 0 */
    3318             :     }
    3319             :     /* W = (Z_K/pr)-basis of O/pr. O = (W0,I0) ~ (W, I) */
    3320             : 
    3321             :    /* compute MW: W_i*W_j = sum MW_k,(i,j) W_k */
    3322        5382 :     Wa = RgM_to_RgXV(W,vpol);
    3323        5382 :     Winv = nfM_inv(nf, W);
    3324        5382 :     MW = cgetg(nn+1, t_MAT);
    3325             :     /* W_1 = 1 */
    3326       31854 :     for (j=1; j<=n; j++) gel(MW, j) = gel(MW, (j-1)*n+1) = gel(Id,j);
    3327       26472 :     for (i=2; i<=n; i++)
    3328       81835 :       for (j=i; j<=n; j++)
    3329             :       {
    3330       60745 :         GEN z = nfXQ_mul(nf, gel(Wa,i), gel(Wa,j), pol);
    3331       60745 :         if (typ(z) != t_POL)
    3332           0 :           z = nfC_nf_mul(nf, gel(Winv,1), z);
    3333             :         else
    3334             :         {
    3335       60745 :           z = RgX_to_RgC(z, lg(Winv)-1);
    3336       60745 :           z = nfM_nfC_mul(nf, Winv, z);
    3337             :         }
    3338       60745 :         gel(MW, (i-1)*n+j) = gel(MW, (j-1)*n+i) = z;
    3339             :       }
    3340             : 
    3341             :     /* compute Ip =  pr-radical [ could use Ker(trace) if q large ] */
    3342        5382 :     MWmod = nfM_to_FqM(MW,nf,modpr);
    3343        5382 :     F = cgetg(n+1, t_MAT); gel(F,1) = gel(Id,1);
    3344       26472 :     for (j=2; j<=n; j++) gel(F,j) = rnfeltid_powmod(MWmod, j, q1, T,p);
    3345        5382 :     Ip = FqM_ker(F,T,p);
    3346        5382 :     if (lg(Ip) == 1) { W = W0; I = I0; break; }
    3347             : 
    3348             :     /* Fill C: W_k A_j = sum_i C_(i,j),k A_i */
    3349        4962 :     A = FqM_to_nfM(FqM_suppl(Ip,T,p), modpr);
    3350       13088 :     for (j = lg(Ip); j<=n; j++) gel(A,j) = nfC_multable_mul(gel(A,j), mpi);
    3351        4962 :     MW = nfM_mul(nf, nfM_inv(nf,A), MW);
    3352        4962 :     C = cgetg(n+1, t_MAT);
    3353       29684 :     for (k=1; k<=n; k++)
    3354             :     {
    3355       24722 :       GEN mek = vecslice(MW, (k-1)*n+1, k*n), Ck;
    3356       24722 :       gel(C,k) = Ck = cgetg(nn+1, t_COL);
    3357      163416 :       for (j=1; j<=n; j++)
    3358             :       {
    3359      138694 :         GEN z = nfM_nfC_mul(nf, mek, gel(A,j));
    3360      975186 :         for (i=1; i<=n; i++) gel(Ck, (j-1)*n+i) = nf_to_Fq(nf,gel(z,i),modpr);
    3361             :       }
    3362             :     }
    3363        4962 :     G = FqM_to_nfM(FqM_ker(C,T,p), modpr);
    3364             : 
    3365        4962 :     pseudo = rnfjoinmodules_i(nf, G,prhinv, Id,I);
    3366             :     /* express W in terms of the power basis */
    3367        4962 :     W = nfM_mul(nf, W, gel(pseudo,1));
    3368        4962 :     I = gel(pseudo,2);
    3369             :     /* restore the HNF property W[i,i] = 1. NB: W upper triangular, with
    3370             :      * W[i,i] = Tau[i] */
    3371       29684 :     for (j=1; j<=n; j++)
    3372       24722 :       if (gel(Tau,j) != gen_1)
    3373             :       {
    3374       13985 :         gel(W,j) = nfC_nf_mul(nf, gel(W,j), gel(Tauinv,j));
    3375       13985 :         gel(I,j) = idealmul(nf, gel(Tau,j), gel(I,j));
    3376             :       }
    3377        4962 :     if (DEBUGLEVEL>3) err_printf(" new order:\n%Ps\n%Ps\n", W, I);
    3378        4962 :     if (sep <= 3 || gequal(I,I0)) break;
    3379             : 
    3380        4144 :     if (gc_needed(av1,2))
    3381             :     {
    3382           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"rnfmaxord");
    3383           0 :       gerepileall(av1,2, &W,&I);
    3384             :     }
    3385             :   }
    3386        1238 :   return gerepilecopy(av, mkvec2(W, I));
    3387             : }
    3388             : 
    3389             : GEN
    3390      842260 : Rg_nffix(const char *f, GEN T, GEN c, int lift)
    3391             : {
    3392      842260 :   switch(typ(c))
    3393             :   {
    3394      461831 :     case t_INT: case t_FRAC: return c;
    3395       67342 :     case t_POL:
    3396       67342 :       if (lg(c) >= lg(T)) c = RgX_rem(c,T);
    3397       67342 :       break;
    3398      313081 :     case t_POLMOD:
    3399      313081 :       if (!RgX_equal_var(gel(c,1), T)) pari_err_MODULUS(f, gel(c,1),T);
    3400      312500 :       c = gel(c,2);
    3401      312500 :       switch(typ(c))
    3402             :       {
    3403      275998 :         case t_POL: break;
    3404       36502 :         case t_INT: case t_FRAC: return c;
    3405           0 :         default: pari_err_TYPE(f, c);
    3406             :       }
    3407      275998 :       break;
    3408           6 :     default: pari_err_TYPE(f,c);
    3409             :   }
    3410             :   /* typ(c) = t_POL */
    3411      343340 :   if (varn(c) != varn(T)) pari_err_VAR(f, c,T);
    3412      343326 :   switch(lg(c))
    3413             :   {
    3414       13311 :     case 2: return gen_0;
    3415       27957 :     case 3:
    3416       27957 :       c = gel(c,2); if (is_rational_t(typ(c))) return c;
    3417           0 :       pari_err_TYPE(f,c);
    3418             :   }
    3419      302058 :   RgX_check_QX(c, f);
    3420      302036 :   return lift? c: mkpolmod(c, T);
    3421             : }
    3422             : /* check whether P is a polynomials with coeffs in number field Q[y]/(T) */
    3423             : GEN
    3424      279456 : RgX_nffix(const char *f, GEN T, GEN P, int lift)
    3425             : {
    3426      279456 :   long i, l, vT = varn(T);
    3427      279456 :   GEN Q = cgetg_copy(P, &l);
    3428      279456 :   if (typ(P) != t_POL) pari_err_TYPE(stack_strcat(f," [t_POL expected]"), P);
    3429      279456 :   if (varncmp(varn(P), vT) >= 0) pari_err_PRIORITY(f, P, ">=", vT);
    3430      279435 :   Q[1] = P[1];
    3431     1078013 :   for (i=2; i<l; i++) gel(Q,i) = Rg_nffix(f, T, gel(P,i), lift);
    3432      279421 :   return normalizepol_lg(Q, l);
    3433             : }
    3434             : GEN
    3435          49 : RgV_nffix(const char *f, GEN T, GEN P, int lift)
    3436             : {
    3437             :   long i, l;
    3438          49 :   GEN Q = cgetg_copy(P, &l);
    3439         119 :   for (i=1; i<l; i++) gel(Q,i) = Rg_nffix(f, T, gel(P,i), lift);
    3440          42 :   return Q;
    3441             : }
    3442             : 
    3443             : static GEN
    3444        2240 : get_d(GEN nf, GEN d)
    3445             : {
    3446        2240 :   GEN b = idealredmodpower(nf, d, 2, 100000);
    3447        2240 :   return nfmul(nf, d, nfsqr(nf,b));
    3448             : }
    3449             : 
    3450             : /* true nf */
    3451             : static GEN
    3452        3535 : pr_factorback(GEN nf, GEN fa)
    3453             : {
    3454        3535 :   GEN P = gel(fa,1), E = gel(fa,2), z = gen_1;
    3455        3535 :   long i, l = lg(P);
    3456        8330 :   for (i = 1; i < l; i++) z = idealmulpowprime(nf, z, gel(P,i), gel(E,i));
    3457        3535 :   return z;
    3458             : }
    3459             : /* true nf */
    3460             : static GEN
    3461        3535 : pr_factorback_scal(GEN nf, GEN fa)
    3462             : {
    3463        3535 :   GEN D = pr_factorback(nf,fa);
    3464        3535 :   if (typ(D) == t_MAT && RgM_isscalar(D,NULL)) D = gcoeff(D,1,1);
    3465        3535 :   return D;
    3466             : }
    3467             : 
    3468             : /* nf = base field K
    3469             :  * pol= monic polynomial in Z_K[X] defining a relative extension L = K[X]/(pol).
    3470             :  * Returns a pseudo-basis [A,I] of Z_L, set *pD to [D,d] and *pf to the
    3471             :  * index-ideal; rnf is used when lim != 0 and may be NULL */
    3472             : GEN
    3473        2184 : rnfallbase(GEN nf, GEN pol, GEN lim, GEN rnf, GEN *pD, GEN *pf, GEN *pDKP)
    3474             : {
    3475             :   long i, j, jf, l;
    3476             :   GEN fa, E, P, Ef, Pf, z, disc;
    3477             : 
    3478        2184 :   nf = checknf(nf); pol = liftpol_shallow(pol);
    3479        2184 :   if (!gequal1(leading_coeff(pol)))
    3480           7 :     pari_err_IMPL("nonmonic relative polynomials in rnfallbase");
    3481        2177 :   disc = nf_to_scalar_or_basis(nf, nfX_disc(nf, pol));
    3482        2177 :   if (gequal0(disc))
    3483           7 :     pari_err_DOMAIN("rnfpseudobasis","issquarefree(pol)","=",gen_0, pol);
    3484        2170 :   if (lim)
    3485             :   {
    3486             :     GEN rnfeq, zknf, dzknf, U, vU, dA, A, MB, dB, BdB, vj, B, Tabs;
    3487         406 :     GEN D = idealhnf_shallow(nf, disc), extendP = NULL;
    3488         406 :     long rU, m = nf_get_degree(nf), n = degpol(pol), N = n*m;
    3489             :     nfmaxord_t S;
    3490             : 
    3491         406 :     if (typ(lim) == t_INT)
    3492         112 :       P = ZV_union_shallow(nf_get_ramified_primes(nf),
    3493         112 :                            gel(Z_factor_limit(gcoeff(D,1,1), itou(lim)), 1));
    3494             :     else
    3495             :     {
    3496         294 :       P = cgetg_copy(lim, &l);
    3497         882 :       for (i = 1; i < l; i++)
    3498             :       {
    3499         588 :         GEN p = gel(lim,i);
    3500         588 :         if (typ(p) != t_INT) p = pr_get_p(p);
    3501         588 :         gel(P,i) = p;
    3502             :       }
    3503         294 :       P = ZV_sort_uniq_shallow(P);
    3504             :     }
    3505         406 :     if (rnf)
    3506             :     {
    3507         357 :       rnfeq = rnf_get_map(rnf);
    3508         357 :       zknf = rnf_get_nfzk(rnf);
    3509             :     }
    3510             :     else
    3511             :     {
    3512          49 :       rnfeq = nf_rnfeq(nf, pol);
    3513          49 :       zknf = nf_nfzk(nf, rnfeq);
    3514             :     }
    3515         406 :     dzknf = gel(zknf,1);
    3516         406 :     if (gequal1(dzknf)) dzknf = NULL;
    3517         273 : RESTART:
    3518         413 :     if (extendP)
    3519             :     {
    3520           7 :       GEN oldP = P;
    3521           7 :       if (typ(extendP)==t_POL)
    3522             :       {
    3523           0 :         long l = lg(extendP);
    3524           0 :         for (i = 2; i < l; i++)
    3525             :         {
    3526           0 :           GEN q = gel(extendP,i);
    3527           0 :           if (typ(q) == t_FRAC) P = ZV_cba_extend(P, gel(q,2));
    3528             :         }
    3529             :       } else /*t_FRAC*/
    3530           7 :         P = ZV_cba_extend(P, gel(extendP,2));
    3531           7 :       if (ZV_equal(P, oldP))
    3532           0 :         pari_err(e_MISC, "rnfpseudobasis fails, try increasing B");
    3533           7 :       extendP = NULL;
    3534             :     }
    3535         413 :     Tabs = gel(rnfeq,1);
    3536         413 :     nfmaxord(&S, mkvec2(Tabs,P), 0);
    3537         413 :     B = RgXV_unscale(S.basis, S.unscale);
    3538         413 :     BdB = Q_remove_denom(B, &dB);
    3539         413 :     MB = RgXV_to_RgM(BdB, N); /* HNF */
    3540             : 
    3541         413 :     vU = cgetg(N+1, t_VEC);
    3542         413 :     vj = cgetg(N+1, t_VECSMALL);
    3543         413 :     gel(vU,1) = U = cgetg(m+1, t_MAT);
    3544         413 :     gel(U,1) = col_ei(N, 1);
    3545         413 :     A = dB? (dzknf? gdiv(dB,dzknf): dB): NULL;
    3546         413 :     if (A)
    3547             :     {
    3548         364 :       if (typ(A) != t_INT) { extendP = A; goto RESTART; }
    3549         357 :       if (equali1(A)) A = NULL;
    3550             :     }
    3551         784 :     for (j = 2; j <= m; j++)
    3552             :     {
    3553         378 :       GEN t = gel(zknf,j);
    3554         378 :       if (!RgX_is_ZX(t)) { extendP = t; goto RESTART; }
    3555         378 :       if (A) t = ZX_Z_mul(t, A);
    3556         378 :       gel(U,j) = hnf_solve(MB, RgX_to_RgC(t, N));
    3557             :     }
    3558        2450 :     for (i = 2; i <= N; i++)
    3559             :     {
    3560        2044 :       GEN b = gel(BdB,i);
    3561        2044 :       gel(vU,i) = U = cgetg(m+1, t_MAT);
    3562        2044 :       gel(U,1) = hnf_solve(MB, RgX_to_RgC(b, N));
    3563        4354 :       for (j = 2; j <= m; j++)
    3564             :       {
    3565        2310 :         GEN t = ZX_rem(ZX_mul(b, gel(zknf,j)), Tabs);
    3566        2310 :         if (dzknf)
    3567             :         {
    3568        1694 :           t = RgX_Rg_div(t, dzknf);
    3569        1694 :           if (!RgX_is_ZX(t)) { extendP = t; goto RESTART; }
    3570             :         }
    3571        2310 :         gel(U,j) = hnf_solve(MB, RgX_to_RgC(t, N));
    3572             :       }
    3573             :     }
    3574         406 :     vj[1] = 1; U = gel(vU,1); rU = m;
    3575         945 :     for (i = j = 2; i <= N; i++)
    3576             :     {
    3577         938 :       GEN V = shallowconcat(U, gel(vU,i));
    3578         938 :       if (ZM_rank(V) != rU)
    3579             :       {
    3580         938 :         U = V; rU += m; vj[j++] = i;
    3581         938 :         if (rU == N) break;
    3582             :       }
    3583             :     }
    3584         406 :     if (dB) for(;;)
    3585         420 :     {
    3586         777 :       GEN c = gen_1, H = ZM_hnfmodid(U, dB);
    3587         777 :       long ic = 0;
    3588        6650 :       for (i = 1; i <= N; i++)
    3589        5873 :         if (cmpii(gcoeff(H,i,i), c) > 0) { c = gcoeff(H,i,i); ic = i; }
    3590         777 :       if (!ic) break;
    3591         420 :       vj[j++] = ic;
    3592         420 :       U = shallowconcat(H, gel(vU, ic));
    3593             :     }
    3594         406 :     setlg(vj, j);
    3595         406 :     B = vecpermute(B, vj);
    3596             : 
    3597         406 :     l = lg(B);
    3598         406 :     A = cgetg(l,t_MAT);
    3599        2170 :     for (j = 1; j < l; j++)
    3600             :     {
    3601        1764 :       GEN t = eltabstorel_lift(rnfeq, gel(B,j));
    3602        1764 :       gel(A,j) = Rg_to_RgC(t, n);
    3603             :     }
    3604         406 :     A = RgM_to_nfM(nf, A);
    3605         406 :     A = Q_remove_denom(A, &dA);
    3606         406 :     if (!dA)
    3607             :     { /* order is maximal */
    3608          56 :       z = triv_order(n);
    3609          56 :       if (pf) *pf = gen_1;
    3610             :     }
    3611             :     else
    3612             :     {
    3613             :       GEN fi;
    3614             :       /* the first n columns of A are probably in HNF already */
    3615         350 :       A = shallowconcat(vecslice(A,n+1,lg(A)-1), vecslice(A,1,n));
    3616         350 :       A = mkvec2(A, const_vec(l-1,gen_1));
    3617         350 :       if (DEBUGLEVEL > 2) err_printf("rnfallbase: nfhnf in dim %ld\n", l-1);
    3618         350 :       z = nfhnfmod(nf, A, nfdetint(nf,A));
    3619         350 :       gel(z,2) = gdiv(gel(z,2), dA);
    3620         350 :       fi = idealprod(nf,gel(z,2));
    3621         350 :       D = idealmul(nf, D, idealsqr(nf, fi));
    3622         350 :       if (pf) *pf = idealinv(nf, fi);
    3623             :     }
    3624         406 :     if (RgM_isscalar(D,NULL)) D = gcoeff(D,1,1);
    3625         406 :     if (pDKP) *pDKP = S.dKP;
    3626         406 :     *pD = mkvec2(D, get_d(nf, disc)); return z;
    3627             :   }
    3628        1764 :   fa = idealfactor(nf, disc);
    3629        1764 :   P = gel(fa,1); l = lg(P); z = NULL;
    3630        1764 :   E = gel(fa,2);
    3631        1764 :   Pf = cgetg(l, t_COL);
    3632        1764 :   Ef = cgetg(l, t_COL);
    3633        5895 :   for (i = j = jf = 1; i < l; i++)
    3634             :   {
    3635        4138 :     GEN pr = gel(P,i);
    3636        4138 :     long e = itos(gel(E,i));
    3637        4138 :     if (e > 1)
    3638             :     {
    3639        3171 :       GEN vD = rnfmaxord(nf, pol, pr, e);
    3640        3164 :       if (vD)
    3641             :       {
    3642        2134 :         long ef = idealprodval(nf, gel(vD,2), pr);
    3643        2134 :         z = rnfjoinmodules(nf, z, vD);
    3644        2134 :         if (ef) { gel(Pf, jf) = pr; gel(Ef, jf++) = stoi(-ef); }
    3645        2134 :         e += 2 * ef;
    3646             :       }
    3647             :     }
    3648        4131 :     if (e) { gel(P, j) = pr; gel(E, j++) = stoi(e); }
    3649             :   }
    3650        1757 :   setlg(P,j);
    3651        1757 :   setlg(E,j);
    3652        1757 :   if (pDKP) *pDKP = prV_primes(P);
    3653        1757 :   if (pf)
    3654             :   {
    3655        1701 :     setlg(Pf, jf);
    3656        1701 :     setlg(Ef, jf); *pf = pr_factorback_scal(nf, mkmat2(Pf,Ef));
    3657             :   }
    3658        1757 :   *pD = mkvec2(pr_factorback_scal(nf,fa), get_d(nf, disc));
    3659        1757 :   return z? z: triv_order(degpol(pol));
    3660             : }
    3661             : 
    3662             : static GEN
    3663        1631 : RgX_to_algX(GEN nf, GEN x)
    3664             : {
    3665             :   long i, l;
    3666        1631 :   GEN y = cgetg_copy(x, &l); y[1] = x[1];
    3667        8687 :   for (i=2; i<l; i++) gel(y,i) = nf_to_scalar_or_alg(nf, gel(x,i));
    3668        1631 :   return y;
    3669             : }
    3670             : 
    3671             : GEN
    3672        1645 : nfX_to_monic(GEN nf, GEN T, GEN *pL)
    3673             : {
    3674             :   GEN lT, g, a;
    3675        1645 :   long i, l = lg(T);
    3676        1645 :   if (l == 2) return pol_0(varn(T));
    3677        1645 :   if (l == 3) return pol_1(varn(T));
    3678        1645 :   nf = checknf(nf);
    3679        1645 :   T = Q_primpart(RgX_to_nfX(nf, T));
    3680        1645 :   lT = leading_coeff(T); if (pL) *pL = lT;
    3681        1645 :   if (isint1(T)) return T;
    3682        1645 :   g = cgetg_copy(T, &l); g[1] = T[1]; a = lT;
    3683        1645 :   gel(g, l-1) = gen_1;
    3684        1645 :   gel(g, l-2) = gel(T,l-2);
    3685        1645 :   if (l == 4) { gel(g,l-2) = nf_to_scalar_or_alg(nf, gel(g,l-2)); return g; }
    3686        1631 :   if (typ(lT) == t_INT)
    3687             :   {
    3688        1617 :     gel(g, l-3) = gmul(a, gel(T,l-3));
    3689        3773 :     for (i = l-4; i > 1; i--) { a = mulii(a,lT); gel(g,i) = gmul(a, gel(T,i)); }
    3690             :   }
    3691             :   else
    3692             :   {
    3693          14 :     gel(g, l-3) = nfmul(nf, a, gel(T,l-3));
    3694          35 :     for (i = l-3; i > 1; i--)
    3695             :     {
    3696          21 :       a = nfmul(nf,a,lT);
    3697          21 :       gel(g,i) = nfmul(nf, a, gel(T,i));
    3698             :     }
    3699             :   }
    3700        1631 :   return RgX_to_algX(nf, g);
    3701             : }
    3702             : 
    3703             : GEN
    3704         854 : rnfdisc_factored(GEN nf, GEN pol, GEN *pd)
    3705             : {
    3706             :   long i, j, l;
    3707             :   GEN fa, E, P, disc, lim;
    3708             : 
    3709         854 :   pol = rnfdisc_get_T(nf, pol, &lim);
    3710         854 :   disc = nf_to_scalar_or_basis(nf, nfX_disc(nf, pol));
    3711         854 :   if (gequal0(disc))
    3712           0 :     pari_err_DOMAIN("rnfdisc","issquarefree(pol)","=",gen_0, pol);
    3713         854 :   pol = nfX_to_monic(nf, pol, NULL);
    3714         854 :   fa = idealfactor_partial(nf, disc, lim);
    3715         854 :   P = gel(fa,1); l = lg(P);
    3716         854 :   E = gel(fa,2);
    3717        2310 :   for (i = j = 1; i < l; i++)
    3718             :   {
    3719        1456 :     long e = itos(gel(E,i));
    3720        1456 :     GEN pr = gel(P,i);
    3721        1456 :     if (e > 1)
    3722             :     {
    3723        1183 :       GEN vD = rnfmaxord(nf, pol, pr, e);
    3724        1183 :       if (vD) e += 2*idealprodval(nf, gel(vD,2), pr);
    3725             :     }
    3726        1456 :     if (e) { gel(P, j) = pr; gel(E, j++) = stoi(e); }
    3727             :   }
    3728         854 :   if (pd) *pd = get_d(nf, disc);
    3729         854 :   setlg(P, j);
    3730         854 :   setlg(E, j); return fa;
    3731             : }
    3732             : GEN
    3733          77 : rnfdiscf(GEN nf, GEN pol)
    3734             : {
    3735          77 :   pari_sp av = avma;
    3736             :   GEN d, fa;
    3737          77 :   nf = checknf(nf); fa = rnfdisc_factored(nf, pol, &d);
    3738          77 :   return gerepilecopy(av, mkvec2(pr_factorback_scal(nf,fa), d));
    3739             : }
    3740             : 
    3741             : GEN
    3742          35 : gen_if_principal(GEN bnf, GEN x)
    3743             : {
    3744          35 :   pari_sp av = avma;
    3745          35 :   GEN z = bnfisprincipal0(bnf,x, nf_GEN_IF_PRINCIPAL | nf_FORCE);
    3746          35 :   return isintzero(z)? gc_NULL(av): z;
    3747             : }
    3748             : 
    3749             : /* given bnf and a HNF pseudo-basis of a proj. module, simplify the HNF as
    3750             :  * much as possible. The resulting matrix will be upper triangular but the
    3751             :  * diagonal coefficients will not be equal to 1. The ideals are integral and
    3752             :  * primitive. */
    3753             : GEN
    3754           0 : rnfsimplifybasis(GEN bnf, GEN M)
    3755             : {
    3756           0 :   pari_sp av = avma;
    3757             :   long i, l;
    3758             :   GEN y, Az, Iz, nf, A, I;
    3759             : 
    3760           0 :   bnf = checkbnf(bnf); nf = bnf_get_nf(bnf);
    3761           0 :   if (!check_ZKmodule_i(M)) pari_err_TYPE("rnfsimplifybasis",M);
    3762           0 :   A = gel(M,1);
    3763           0 :   I = gel(M,2); l = lg(I);
    3764           0 :   Az = cgetg(l, t_MAT);
    3765           0 :   Iz = cgetg(l, t_VEC); y = mkvec2(Az, Iz);
    3766           0 :   for (i = 1; i < l; i++)
    3767             :   {
    3768             :     GEN c, d;
    3769           0 :     if (ideal_is1(gel(I,i)))
    3770             :     {
    3771           0 :       gel(Iz,i) = gen_1;
    3772           0 :       gel(Az,i) = gel(A,i); continue;
    3773             :     }
    3774             : 
    3775           0 :     gel(Iz,i) = Q_primitive_part(gel(I,i), &c);
    3776           0 :     gel(Az,i) = c? RgC_Rg_mul(gel(A,i),c): gel(A,i);
    3777           0 :     if (c && ideal_is1(gel(Iz,i))) continue;
    3778             : 
    3779           0 :     d = gen_if_principal(bnf, gel(Iz,i));
    3780           0 :     if (d)
    3781             :     {
    3782           0 :       gel(Iz,i) = gen_1;
    3783           0 :       gel(Az,i) = nfC_nf_mul(nf, gel(Az,i), d);
    3784             :     }
    3785             :   }
    3786           0 :   return gerepilecopy(av, y);
    3787             : }
    3788             : 
    3789             : static GEN
    3790          63 : get_module(GEN nf, GEN O, const char *s)
    3791             : {
    3792          63 :   if (typ(O) == t_POL) return rnfpseudobasis(nf, O);
    3793          56 :   if (!check_ZKmodule_i(O)) pari_err_TYPE(s, O);
    3794          56 :   return shallowcopy(O);
    3795             : }
    3796             : 
    3797             : GEN
    3798          14 : rnfdet(GEN nf, GEN M)
    3799             : {
    3800          14 :   pari_sp av = avma;
    3801             :   GEN D;
    3802          14 :   nf = checknf(nf);
    3803          14 :   M = get_module(nf, M, "rnfdet");
    3804          14 :   D = idealmul(nf, nfM_det(nf, gel(M,1)), idealprod(nf, gel(M,2)));
    3805          14 :   return gerepileupto(av, D);
    3806             : }
    3807             : 
    3808             : /* Given two fractional ideals a and b, gives x in a, y in b, z in b^-1,
    3809             :    t in a^-1 such that xt-yz=1. In the present version, z is in Z. */
    3810             : static void
    3811          63 : nfidealdet1(GEN nf, GEN a, GEN b, GEN *px, GEN *py, GEN *pz, GEN *pt)
    3812             : {
    3813             :   GEN x, uv, y, da, db;
    3814             : 
    3815          63 :   a = idealinv(nf,a);
    3816          63 :   a = Q_remove_denom(a, &da);
    3817          63 :   b = Q_remove_denom(b, &db);
    3818          63 :   x = idealcoprime(nf,a,b);
    3819          63 :   uv = idealaddtoone(nf, idealmul(nf,x,a), b);
    3820          63 :   y = gel(uv,2);
    3821          63 :   if (da) x = gmul(x,da);
    3822          63 :   if (db) y = gdiv(y,db);
    3823          63 :   *px = x;
    3824          63 :   *py = y;
    3825          63 :   *pz = db ? negi(db): gen_m1;
    3826          63 :   *pt = nfdiv(nf, gel(uv,1), x);
    3827          63 : }
    3828             : 
    3829             : /* given a pseudo-basis of a proj. module in HNF [A,I] (or [A,I,D,d]), gives
    3830             :  * an n x n matrix (not HNF) of a pseudo-basis and an ideal vector
    3831             :  * [1,...,1,I] such that M ~ Z_K^(n-1) x I. Uses the approximation theorem.*/
    3832             : GEN
    3833          28 : rnfsteinitz(GEN nf, GEN M)
    3834             : {
    3835          28 :   pari_sp av = avma;
    3836             :   long i, n;
    3837             :   GEN A, I;
    3838             : 
    3839          28 :   nf = checknf(nf);
    3840          28 :   M = get_module(nf, M, "rnfsteinitz");
    3841          28 :   A = RgM_to_nfM(nf, gel(M,1));
    3842          28 :   I = leafcopy(gel(M,2)); n = lg(A)-1;
    3843         189 :   for (i = 1; i < n; i++)
    3844             :   {
    3845         161 :     GEN c1, c2, b, a = gel(I,i);
    3846         161 :     gel(I,i) = gen_1;
    3847         161 :     if (ideal_is1(a)) continue;
    3848             : 
    3849          63 :     c1 = gel(A,i);
    3850          63 :     c2 = gel(A,i+1);
    3851          63 :     b = gel(I,i+1);
    3852          63 :     if (ideal_is1(b))
    3853             :     {
    3854           0 :       gel(A,i) = c2;
    3855           0 :       gel(A,i+1) = gneg(c1);
    3856           0 :       gel(I,i+1) = a;
    3857             :     }
    3858             :     else
    3859             :     {
    3860          63 :       pari_sp av2 = avma;
    3861             :       GEN x, y, z, t, c;
    3862          63 :       nfidealdet1(nf,a,b, &x,&y,&z,&t);
    3863          63 :       x = RgC_add(nfC_nf_mul(nf, c1, x), nfC_nf_mul(nf, c2, y));
    3864          63 :       y = RgC_add(nfC_nf_mul(nf, c1, z), nfC_nf_mul(nf, c2, t));
    3865          63 :       gerepileall(av2, 2, &x,&y);
    3866          63 :       gel(A,i) = x;
    3867          63 :       gel(A,i+1) = y;
    3868          63 :       gel(I,i+1) = Q_primitive_part(idealmul(nf,a,b), &c);
    3869          63 :       if (c) gel(A,i+1) = nfC_nf_mul(nf, gel(A,i+1), c);
    3870             :     }
    3871             :   }
    3872          28 :   gel(M,1) = A;
    3873          28 :   gel(M,2) = I; return gerepilecopy(av, M);
    3874             : }
    3875             : 
    3876             : /* Given bnf and a proj. module (or a t_POL -> rnfpseudobasis), and outputs a
    3877             :  * basis if it is free, an n+1-generating set if it is not */
    3878             : GEN
    3879          21 : rnfbasis(GEN bnf, GEN M)
    3880             : {
    3881          21 :   pari_sp av = avma;
    3882             :   long j, n;
    3883             :   GEN nf, A, I, cl, col, a;
    3884             : 
    3885          21 :   bnf = checkbnf(bnf); nf = bnf_get_nf(bnf);
    3886          21 :   M = get_module(nf, M, "rnfbasis");
    3887          21 :   I = gel(M,2); n = lg(I)-1;
    3888          98 :   j = 1; while (j < n && ideal_is1(gel(I,j))) j++;
    3889          21 :   if (j < n) { M = rnfsteinitz(nf,M); I = gel(M,2); }
    3890          21 :   A = gel(M,1);
    3891          21 :   col= gel(A,n); A = vecslice(A, 1, n-1);
    3892          21 :   cl = gel(I,n);
    3893          21 :   a = gen_if_principal(bnf, cl);
    3894          21 :   if (!a)
    3895             :   {
    3896           7 :     GEN v = idealtwoelt(nf, cl);
    3897           7 :     A = vec_append(A, gmul(gel(v,1), col));
    3898           7 :     a = gel(v,2);
    3899             :   }
    3900          21 :   A = vec_append(A, nfC_nf_mul(nf, col, a));
    3901          21 :   return gerepilecopy(av, A);
    3902             : }
    3903             : 
    3904             : /* Given a Z_K-module M (or a polynomial => rnfpseudobasis) outputs a
    3905             :  * Z_K-basis in HNF if it exists, zero if not */
    3906             : GEN
    3907           7 : rnfhnfbasis(GEN bnf, GEN M)
    3908             : {
    3909           7 :   pari_sp av = avma;
    3910             :   long j, l;
    3911             :   GEN nf, A, I, a;
    3912             : 
    3913           7 :   bnf = checkbnf(bnf); nf = bnf_get_nf(bnf);
    3914           7 :   if (typ(M) == t_POL) M = rnfpseudobasis(nf, M);
    3915             :   else
    3916             :   {
    3917           7 :     if (typ(M) != t_VEC) pari_err_TYPE("rnfhnfbasis", M);
    3918           7 :     if (lg(M) == 5) M = mkvec2(gel(M,1), gel(M,2));
    3919           7 :     M = nfhnf(nf, M); /* in case M is not in HNF */
    3920             :   }
    3921           7 :   A = shallowcopy(gel(M,1));
    3922           7 :   I = gel(M,2); l = lg(A);
    3923          42 :   for (j = 1; j < l; j++)
    3924             :   {
    3925          35 :     if (ideal_is1(gel(I,j))) continue;
    3926          14 :     a = gen_if_principal(bnf, gel(I,j));
    3927          14 :     if (!a) return gc_const(av, gen_0);
    3928          14 :     gel(A,j) = nfC_nf_mul(nf, gel(A,j), a);
    3929             :   }
    3930           7 :   return gerepilecopy(av,A);
    3931             : }
    3932             : 
    3933             : long
    3934           7 : rnfisfree(GEN bnf, GEN M)
    3935             : {
    3936           7 :   pari_sp av = avma;
    3937             :   GEN nf, P, I;
    3938             :   long l, j;
    3939             : 
    3940           7 :   bnf = checkbnf(bnf);
    3941           7 :   if (is_pm1( bnf_get_no(bnf) )) return 1;
    3942           0 :   nf = bnf_get_nf(bnf);
    3943           0 :   M = get_module(nf, M, "rnfisfree");
    3944           0 :   I = gel(M,2); l = lg(I); P = NULL;
    3945           0 :   for (j = 1; j < l; j++)
    3946           0 :     if (!ideal_is1(gel(I,j))) P = P? idealmul(nf, P, gel(I,j)): gel(I,j);
    3947           0 :   return gc_long(av, P? gequal0( isprincipal(bnf,P) ): 1);
    3948             : }
    3949             : 
    3950             : /**********************************************************************/
    3951             : /**                                                                  **/
    3952             : /**                   COMPOSITUM OF TWO NUMBER FIELDS                **/
    3953             : /**                                                                  **/
    3954             : /**********************************************************************/
    3955             : static GEN
    3956       26470 : compositum_fix(GEN nf, GEN A)
    3957             : {
    3958             :   int ok;
    3959       26470 :   if (nf)
    3960             :   {
    3961         952 :     A = Q_primpart(liftpol_shallow(A)); RgX_check_ZXX(A,"polcompositum");
    3962         952 :     ok = nfissquarefree(nf,A);
    3963             :   }
    3964             :   else
    3965             :   {
    3966       25518 :     A = Q_primpart(A); RgX_check_ZX(A,"polcompositum");
    3967       25512 :     ok = ZX_is_squarefree(A);
    3968             :   }
    3969       26465 :   if (!ok) pari_err_DOMAIN("polcompositum","issquarefree(arg)","=",gen_0,A);
    3970       26458 :   return A;
    3971             : }
    3972             : #define next_lambda(a) (a>0 ? -a : 1-a)
    3973             : 
    3974             : static long
    3975         497 : nfcompositum_lambda(GEN nf, GEN A, GEN B, long lambda)
    3976             : {
    3977         497 :   pari_sp av = avma;
    3978             :   forprime_t S;
    3979         497 :   GEN T = nf_get_pol(nf);
    3980         497 :   long vT = varn(T);
    3981             :   ulong p;
    3982         497 :   init_modular_big(&S);
    3983         497 :   p = u_forprime_next(&S);
    3984             :   while (1)
    3985          42 :   {
    3986             :     GEN Hp, Tp, a;
    3987         539 :     if (DEBUGLEVEL>4) err_printf("Trying lambda = %ld\n", lambda);
    3988         539 :     a = ZXX_to_FlxX(RgX_rescale(A, stoi(-lambda)), p, vT);
    3989         539 :     Tp = ZX_to_Flx(T, p);
    3990         539 :     Hp = FlxqX_composedsum(a, ZXX_to_FlxX(B, p, vT), Tp, p);
    3991         539 :     if (!FlxqX_is_squarefree(Hp, Tp, p))
    3992          42 :       { lambda = next_lambda(lambda); continue; }
    3993         497 :     if (DEBUGLEVEL>4) err_printf("Final lambda = %ld\n", lambda);
    3994         497 :     return gc_long(av, lambda);
    3995             :   }
    3996             : }
    3997             : 
    3998             : /* modular version */
    3999             : GEN
    4000       13347 : nfcompositum(GEN nf, GEN A, GEN B, long flag)
    4001             : {
    4002       13347 :   pari_sp av = avma;
    4003             :   int same;
    4004             :   long v, k;
    4005             :   GEN C, D, LPRS;
    4006             : 
    4007       13347 :   if (typ(A)!=t_POL) pari_err_TYPE("polcompositum",A);
    4008       13347 :   if (typ(B)!=t_POL) pari_err_TYPE("polcompositum",B);
    4009       13347 :   if (degpol(A)<=0 || degpol(B)<=0) pari_err_CONSTPOL("polcompositum");
    4010       13347 :   v = varn(A);
    4011       13347 :   if (varn(B) != v) pari_err_VAR("polcompositum", A,B);
    4012       13347 :   if (nf)
    4013             :   {
    4014         539 :     nf = checknf(nf);
    4015         532 :     if (varncmp(v,nf_get_varn(nf))>=0) pari_err_PRIORITY("polcompositum", nf, ">=",  v);
    4016             :   }
    4017       13305 :   same = (A == B || RgX_equal(A,B));
    4018       13305 :   A = compositum_fix(nf,A);
    4019       13299 :   B = same ? A: compositum_fix(nf,B);
    4020             : 
    4021       13295 :   D = LPRS = NULL; /* -Wall */
    4022       13295 :   k = same? -1: 1;
    4023       13295 :   if (nf)
    4024             :   {
    4025         497 :     long v0 = fetch_var();
    4026         497 :     GEN q, T = nf_get_pol(nf);
    4027         497 :     A = liftpol_shallow(A);
    4028         497 :     B = liftpol_shallow(B);
    4029         497 :     k = nfcompositum_lambda(nf, A, B, k);
    4030         497 :     if (flag&1)
    4031             :     {
    4032             :       GEN H0, H1;
    4033         189 :       GEN chgvar = deg1pol_shallow(stoi(k),pol_x(v0),v);
    4034         189 :       GEN B1 = poleval(QXQX_to_mod_shallow(B, T), chgvar);
    4035         189 :       C = RgX_resultant_all(QXQX_to_mod_shallow(A, T), B1, &q);
    4036         189 :       C = gsubst(C,v0,pol_x(v));
    4037         189 :       C = lift_if_rational(C);
    4038         189 :       H0 = gsubst(gel(q,2),v0,pol_x(v));
    4039         189 :       H1 = gsubst(gel(q,3),v0,pol_x(v));
    4040         189 :       if (typ(H0) != t_POL) H0 = scalarpol_shallow(H0,v);
    4041         189 :       if (typ(H1) != t_POL) H1 = scalarpol_shallow(H1,v);
    4042         189 :       H0 = lift_if_rational(H0);
    4043         189 :       H1 = lift_if_rational(H1);
    4044         189 :       LPRS = mkvec2(H0,H1);
    4045             :     }
    4046             :     else
    4047             :     {
    4048         308 :       C = nf_direct_compositum(nf, RgX_rescale(A,stoi(-k)), B);
    4049         308 :       setvarn(C, v); C = QXQX_to_mod_shallow(C, T);
    4050             :     }
    4051         497 :     C = RgX_normalize(C);
    4052             :   }
    4053             :   else
    4054             :   {
    4055       12798 :     B = leafcopy(B); setvarn(B,fetch_var_higher());
    4056        3080 :     C = (flag&1)? ZX_ZXY_resultant_all(A, B, &k, &LPRS)
    4057       12798 :                 : ZX_compositum(A, B, &k);
    4058       12802 :     setvarn(C, v);
    4059             :   }
    4060             :   /* C = Res_Y (A(Y), B(X + kY)) guaranteed squarefree */
    4061       13299 :   if (flag & 2)
    4062       10254 :     C = mkvec(C);
    4063             :   else
    4064             :   {
    4065        3045 :     if (same)
    4066             :     {
    4067         105 :       D = RgX_rescale(A, stoi(1 - k));
    4068         105 :       if (nf) D = RgX_normalize(QXQX_to_mod_shallow(D, nf_get_pol(nf)));
    4069         105 :       C = RgX_div(C, D);
    4070         105 :       if (degpol(C) <= 0)
    4071           0 :         C = mkvec(D);
    4072             :       else
    4073         105 :         C = shallowconcat(nf? gel(nffactor(nf,C),1): ZX_DDF(C), D);
    4074             :     }
    4075             :     else
    4076        2940 :       C = nf? gel(nffactor(nf,C),1): ZX_DDF(C);
    4077             :   }
    4078       13297 :   gen_sort_inplace(C, (void*)(nf?&cmp_RgX: &cmpii), &gen_cmp_RgX, NULL);
    4079       13297 :   if (flag&1)
    4080             :   { /* a,b,c root of A,B,C = compositum, c = b - k a */
    4081        3269 :     long i, l = lg(C);
    4082        3269 :     GEN a, b, mH0 = RgX_neg(gel(LPRS,1)), H1 = gel(LPRS,2);
    4083        3269 :     setvarn(mH0,v);
    4084        3269 :     setvarn(H1,v);
    4085        6615 :     for (i=1; i<l; i++)
    4086             :     {
    4087        3346 :       GEN D = gel(C,i);
    4088        3346 :       a = RgXQ_mul(mH0, nf? RgXQ_inv(H1,D): QXQ_inv(H1,D), D);
    4089        3346 :       b = gadd(pol_x(v), gmulsg(k,a));
    4090        3346 :       if (degpol(D) == 1) b = RgX_rem(b,D);
    4091        3346 :       gel(C,i) = mkvec4(D, mkpolmod(a,D), mkpolmod(b,D), stoi(-k));
    4092             :     }
    4093             :   }
    4094       13297 :   (void)delete_var();
    4095       13297 :   settyp(C, t_VEC);
    4096       13297 :   if (flag&2) C = gel(C,1);
    4097       13297 :   return gerepilecopy(av, C);
    4098             : }
    4099             : GEN
    4100       12808 : polcompositum0(GEN A, GEN B, long flag)
    4101       12808 : { return nfcompositum(NULL,A,B,flag); }
    4102             : 
    4103             : GEN
    4104          91 : compositum(GEN pol1,GEN pol2) { return polcompositum0(pol1,pol2,0); }
    4105             : GEN
    4106        2828 : compositum2(GEN pol1,GEN pol2) { return polcompositum0(pol1,pol2,1); }

Generated by: LCOV version 1.14