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 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.10.0 lcov report (development 22303-eb3e11d) Lines: 2010 2125 94.6 %
Date: 2018-04-21 06:16:28 Functions: 160 163 98.2 %
Legend: Lines: hit not hit

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

Generated by: LCOV version 1.11