Code coverage tests

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

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

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

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

Generated by: LCOV version 1.11