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-bordeaux1.fr machine (x86_64 architecture), and agregate them in the final report:

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

LCOV - code coverage report
Current view: top level - basemath - base2.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.8.0 lcov report (development 16962-5a32637) Lines: 1810 2116 85.5 %
Date: 2014-10-29 Functions: 144 154 93.5 %
Legend: Lines: hit not hit | Branches: + taken - not taken # not executed Branches: 986 1362 72.4 %

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

Generated by: LCOV version 1.9