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.8.0 lcov report (development 19207-2ed2f69) Lines: 1895 2212 85.7 %
Date: 2016-07-25 07:10:32 Functions: 149 159 93.7 %
Legend: Lines: hit not hit

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

Generated by: LCOV version 1.11