Code coverage tests

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

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

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

LCOV - code coverage report
Current view: top level - basemath - base5.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.10.0 lcov report (development 21343-6216058) Lines: 1031 1114 92.5 %
Date: 2017-11-19 06:21:17 Functions: 72 76 94.7 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2000  The PARI group.
       2             : 
       3             : This file is part of the PARI/GP package.
       4             : 
       5             : PARI/GP is free software; you can redistribute it and/or modify it under the
       6             : terms of the GNU General Public License as published by the Free Software
       7             : Foundation. It is distributed in the hope that it will be useful, but WITHOUT
       8             : ANY WARRANTY WHATSOEVER.
       9             : 
      10             : Check the License for details. You should have received a copy of it, along
      11             : with the package; see the file 'COPYING'. If not, write to the Free Software
      12             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      13             : 
      14             : /*******************************************************************/
      15             : /*                                                                 */
      16             : /*                     RNF STRUCTURE AND OPERATIONS                */
      17             : /*                                                                 */
      18             : /*******************************************************************/
      19             : #include "pari.h"
      20             : #include "paripriv.h"
      21             : 
      22             : /* must return a t_POL */
      23             : GEN
      24       21413 : eltreltoabs(GEN rnfeq, GEN x)
      25             : {
      26             :   long i, k, v;
      27       21413 :   pari_sp av = avma;
      28             :   GEN T, pol, teta, a, s;
      29             : 
      30       21413 :   pol = gel(rnfeq,1);
      31       21413 :   a = gel(rnfeq,2);
      32       21413 :   k = itos(gel(rnfeq,3));
      33       21413 :   T = gel(rnfeq,4);
      34             : 
      35       21413 :   v = varn(pol);
      36       21413 :   if (varncmp(gvar(x), v) > 0) x = scalarpol(x,v);
      37       21413 :   x = RgX_nffix("eltreltoabs", T, x, 1);
      38             :   /* Mod(X - k a, pol(X)), a root of the polynomial defining base */
      39       21406 :   teta = gadd(pol_x(v), gmulsg(-k,a));
      40       21406 :   s = gen_0;
      41       79317 :   for (i=lg(x)-1; i>1; i--)
      42             :   {
      43       57911 :     GEN c = gel(x,i);
      44       57911 :     if (typ(c) == t_POL) c = RgX_RgXQ_eval(c, a, pol);
      45       57911 :     s = RgX_rem(gadd(c, gmul(teta,s)), pol);
      46             :   }
      47       21406 :   return gerepileupto(av, s);
      48             : }
      49             : GEN
      50       41930 : rnfeltreltoabs(GEN rnf,GEN x)
      51             : {
      52       41930 :   const char *f = "rnfeltreltoabs";
      53             :   GEN pol;
      54       41930 :   checkrnf(rnf);
      55       41930 :   pol = rnf_get_polabs(rnf);
      56       41930 :   switch(typ(x))
      57             :   {
      58        6272 :     case t_INT: return icopy(x);
      59          21 :     case t_FRAC: return gcopy(x);
      60             :     case t_POLMOD:
      61       34993 :       if (RgX_equal_var(gel(x,1), pol))
      62             :       { /* already in 'abs' form, unless possibly if nf = Q */
      63       13923 :         if (rnf_get_nfdegree(rnf) == 1)
      64             :         {
      65       13902 :           GEN y = gel(x,2);
      66       13902 :           pari_sp av = avma;
      67       13902 :           y = simplify_shallow(liftpol_shallow(y));
      68       13902 :           return gerepilecopy(av, mkpolmod(y, pol));
      69             :         }
      70          21 :         return gcopy(x);
      71             :       }
      72       21070 :       x = polmod_nffix(f,rnf,x,0);
      73       21056 :       if (typ(x) == t_POLMOD) return rnfeltup(rnf,x);
      74       17402 :       retmkpolmod(eltreltoabs(rnf_get_map(rnf), x), ZX_copy(pol));
      75             :     case t_POL:
      76         595 :       if (varn(x) == rnf_get_nfvarn(rnf)) return rnfeltup(rnf,x);
      77         539 :       retmkpolmod(eltreltoabs(rnf_get_map(rnf), x), ZX_copy(pol));
      78             :   }
      79          49 :   pari_err_TYPE(f,x); return NULL;
      80             : }
      81             : 
      82             : GEN
      83       20174 : eltabstorel_lift(GEN rnfeq, GEN P)
      84             : {
      85       20174 :   GEN k, T = gel(rnfeq,4), relpol = gel(rnfeq,5);
      86       20174 :   if (is_scalar_t(typ(P))) return P;
      87       19348 :   k = gel(rnfeq,3);
      88       19348 :   P = lift_shallow(P);
      89       19348 :   if (signe(k)) P = RgXQX_translate(P, deg1pol_shallow(k, gen_0, varn(T)), T);
      90       19348 :   P = RgXQX_rem(P, relpol, T);
      91       19348 :   return QXQX_to_mod_shallow(P, T);
      92             : }
      93             : /* rnfeq = [pol,a,k,T,relpol], P a t_POL or scalar
      94             :  * Return Mod(P(x + k Mod(y, T(y))), pol(x)) */
      95             : GEN
      96       20090 : eltabstorel(GEN rnfeq, GEN P)
      97             : {
      98       20090 :   GEN T = gel(rnfeq,4), relpol = gel(rnfeq,5);
      99       20090 :   return mkpolmod(eltabstorel_lift(rnfeq,P), QXQX_to_mod_shallow(relpol,T));
     100             : }
     101             : GEN
     102       45038 : rnfeltabstorel(GEN rnf,GEN x)
     103             : {
     104       45038 :   const char *f = "rnfeltabstorel";
     105       45038 :   pari_sp av = avma;
     106             :   GEN pol, T, P, NF;
     107       45038 :   checkrnf(rnf);
     108       45038 :   T = rnf_get_nfpol(rnf);
     109       45038 :   P = rnf_get_pol(rnf);
     110       45038 :   pol = rnf_get_polabs(rnf);
     111       45038 :   switch(typ(x))
     112             :   {
     113        4844 :     case t_INT: return icopy(x);
     114          49 :     case t_FRAC: return gcopy(x);
     115             :     case t_POLMOD:
     116       38829 :       if (RgX_equal_var(P, gel(x,1)))
     117             :       {
     118       13972 :         x = polmod_nffix(f, rnf, x, 0);
     119       13972 :         P = QXQX_to_mod_shallow(P,T);
     120       13972 :         return gerepilecopy(av, mkpolmod(x,P));
     121             :       }
     122       24857 :       if (RgX_equal_var(T, gel(x,1))) { x = Rg_nffix(f, T, x, 0); goto END; }
     123       24773 :       if (!RgX_equal_var(pol, gel(x,1))) pari_err_MODULUS(f, gel(x,1),pol);
     124       24731 :       x = gel(x,2); break;
     125         875 :     case t_POL: break;
     126             :     case t_COL:
     127         441 :       NF = obj_check(rnf, rnf_NFABS);
     128         441 :       if (!NF) pari_err_TYPE("rnfeltabstorel, apply nfinit(rnf)",x);
     129         294 :       x = nf_to_scalar_or_alg(NF,x); break;
     130             :     default:
     131           0 :       pari_err_TYPE(f,x);
     132           0 :       return NULL;
     133             :   }
     134       25900 :   switch(typ(x))
     135             :   {
     136       10024 :     case t_INT: return icopy(x);
     137          56 :     case t_FRAC: return gcopy(x);
     138       15820 :     case t_POL: break;
     139           0 :     default: pari_err_TYPE(f, x);
     140             :   }
     141       15820 :   RgX_check_QX(x,f);
     142       15743 :   if (varn(x) != varn(pol))
     143             :   {
     144          70 :     if (varn(x) == varn(T)) { x = Rg_nffix(f,T,x,0); goto END; }
     145          28 :     pari_err_VAR(f, x,pol);
     146             :   }
     147       15673 :   switch(lg(x))
     148             :   {
     149           0 :     case 2: avma = av; return gen_0;
     150          70 :     case 3: return gerepilecopy(av, gel(x,2));
     151             :   }
     152             : END:
     153       15729 :   return gerepilecopy(av, eltabstorel(rnf_get_map(rnf), x));
     154             : }
     155             : 
     156             : /* x a t_VEC of rnf elements in 'alg' form (t_POL). Assume maximal rank or 0 */
     157             : static GEN
     158        1407 : modulereltoabs(GEN rnf, GEN x)
     159             : {
     160        1407 :   GEN W=gel(x,1), I=gel(x,2), rnfeq = rnf_get_map(rnf), polabs = gel(rnfeq,1);
     161        1407 :   long i, j, k, m, N = lg(W)-1;
     162             :   GEN zknf, czknf, M;
     163             : 
     164        1407 :   if (!N) return cgetg(1, t_VEC);
     165        1344 :   rnf_get_nfzk(rnf, &zknf,&czknf);
     166        1344 :   m = rnf_get_nfdegree(rnf);
     167        1344 :   M = cgetg(N*m+1, t_VEC);
     168        4865 :   for (k=i=1; i<=N; i++)
     169             :   {
     170        3521 :     GEN c0, cid, w = gel(W,i), id = gel(I,i);
     171             : 
     172        3521 :     if (lg(id) == 1) continue; /* must be a t_MAT */
     173        3472 :     id = Q_primitive_part(id, &cid);
     174        3472 :     w = Q_primitive_part(eltreltoabs(rnfeq,w), &c0);
     175        3472 :     c0 = mul_content(c0, mul_content(cid,czknf));
     176        3472 :     if (typ(id) == t_INT)
     177        6587 :       for (j=1; j<=m; j++)
     178             :       {
     179        4340 :         GEN z = RgX_rem(gmul(w, gel(zknf,j)), polabs);
     180        4340 :         if (c0) z = RgX_Rg_mul(z, c0);
     181        4340 :         gel(M,k++) = z;
     182             :       }
     183             :     else
     184        4214 :       for (j=1; j<=m; j++)
     185             :       {
     186        2989 :         GEN c, z = Q_primitive_part(RgV_RgC_mul(zknf,gel(id,j)), &c);
     187        2989 :         z = RgX_rem(gmul(w, z), polabs);
     188        2989 :         c = mul_content(c, c0); if (c) z = RgX_Rg_mul(z, c);
     189        2989 :         gel(M,k++) = z;
     190             :       }
     191             :   }
     192        1344 :   setlg(M, k); return M;
     193             : }
     194             : 
     195             : /* Z-basis for absolute maximal order: [NF.pol, NF.zk] */
     196             : GEN
     197         959 : rnf_zkabs(GEN rnf)
     198             : {
     199         959 :   GEN d, M = modulereltoabs(rnf, rnf_get_zk(rnf));
     200         959 :   GEN T = rnf_get_polabs(rnf);
     201         959 :   long n = degpol(T);
     202         959 :   M = Q_remove_denom(M, &d); /* t_VEC of t_POL */
     203         959 :   if (d)
     204             :   {
     205         665 :     M = RgXV_to_RgM(M,n);
     206         665 :     M = ZM_hnfmodall(M, d, hnf_MODID|hnf_CENTER);
     207         665 :     M = RgM_Rg_div(M, d);
     208             :   }
     209             :   else
     210         294 :     M = matid(n);
     211         959 :   return mkvec2(T, RgM_to_RgXV(M, varn(T)));
     212             : }
     213             : 
     214             : static GEN
     215         904 : mknfabs(GEN rnf, long prec)
     216             : {
     217             :   GEN NF;
     218         904 :   if ((NF = obj_check(rnf,rnf_NFABS)))
     219           1 :   { if (nf_get_prec(NF) < prec) NF = nfnewprec_shallow(NF,prec); }
     220             :   else
     221         903 :     NF = nfinit(rnf_zkabs(rnf), prec);
     222         904 :   return NF;
     223             : }
     224             : 
     225             : static GEN
     226         903 : mkupdown(GEN rnf)
     227             : {
     228         903 :   GEN NF = obj_check(rnf, rnf_NFABS), M, zknf, czknf;
     229             :   long i, m;
     230         903 :   rnf_get_nfzk(rnf, &zknf, &czknf);
     231         903 :   if (isint1(czknf)) czknf = NULL;
     232         903 :   m = lg(zknf)-1; M = cgetg(m+1, t_MAT);
     233         903 :   gel(M,1) = vec_ei(nf_get_degree(NF), 1);
     234        1645 :   for (i = 2; i <= m; i++)
     235             :   {
     236         742 :     GEN c = poltobasis(NF, gel(zknf,i));
     237         742 :     if (czknf) c = gmul(c, czknf);
     238         742 :     gel(M,i) = c;
     239             :   }
     240         903 :   return Qevproj_init(M);
     241             : }
     242             : GEN
     243       22463 : rnf_build_nfabs(GEN rnf, long prec)
     244             : {
     245       22463 :   GEN NF = obj_checkbuild_prec(rnf, rnf_NFABS, &mknfabs, &nf_get_prec, prec);
     246       22463 :   (void)obj_checkbuild(rnf, rnf_MAPS, &mkupdown);
     247       22463 :   return NF;
     248             : }
     249             : 
     250             : void
     251        1533 : rnfcomplete(GEN rnf)
     252        1533 : { (void)rnf_build_nfabs(rnf, nf_get_prec(rnf_get_nf(rnf))); }
     253             : 
     254             : void
     255        1127 : nf_nfzk(GEN nf, GEN rnfeq, GEN *zknf, GEN *czknf)
     256             : {
     257        1127 :   GEN pol = gel(rnfeq,1), a = gel(rnfeq,2);
     258        1127 :   GEN C, zk = nf_get_zkprimpart(nf), D = nf_get_zkden(nf);
     259        1127 :   zk = QXV_QXQ_eval(zk, a, pol);
     260        1127 :   *zknf = Q_primitive_part(zk, &C);
     261        1127 :   if (!C) C = gen_1;
     262        1127 :   *czknf = gdiv(C, D);
     263        1127 : }
     264             : 
     265             : GEN
     266        1127 : rnfinit0(GEN nf, GEN polrel, long flag)
     267             : {
     268        1127 :   pari_sp av = avma;
     269             :   GEN rnf, bas, D,d,f, B, rnfeq, zknf,czknf;
     270        1127 :   nf = checknf(nf);
     271        1127 :   bas = rnfallbase(nf,&polrel, &D,&d, &f);
     272        1120 :   B = matbasistoalg(nf,gel(bas,1));
     273        1120 :   gel(bas,1) = lift_if_rational( RgM_to_RgXV(B,varn(polrel)) );
     274        1120 :   rnfeq = nf_rnfeq(nf,polrel);
     275        1120 :   nf_nfzk(nf, rnfeq, &zknf, &czknf);
     276        1120 :   rnf = obj_init(11, 2);
     277        1120 :   gel(rnf,1) = polrel;
     278        1120 :   gel(rnf,2) = mkvec2(zknf, czknf);
     279        1120 :   gel(rnf,3) = mkvec2(D, d);
     280        1120 :   gel(rnf,4) = f;
     281        1120 :   gel(rnf,5) = cgetg(1, t_VEC); /* dummy */
     282        1120 :   gel(rnf,6) = cgetg(1, t_VEC); /* dummy */
     283        1120 :   gel(rnf,7) = bas;
     284        1120 :   gel(rnf,8) = lift_if_rational( RgM_inv_upper(B) );
     285        1120 :   gel(rnf,9) = typ(f) == t_INT? gen_1: RgM_det_triangular(f);
     286        1120 :   gel(rnf,10)= nf;
     287        1120 :   gel(rnf,11)= rnfeq;
     288        1120 :   rnf = gerepilecopy(av, rnf);
     289        1120 :   if (flag) rnfcomplete(rnf);
     290        1120 :   return rnf;
     291             : }
     292             : GEN
     293         413 : rnfinit(GEN nf, GEN T) { return rnfinit0(nf,T,0); }
     294             : 
     295             : GEN
     296        5586 : rnfeltup0(GEN rnf, GEN x, long flag)
     297             : {
     298        5586 :   pari_sp av = avma;
     299             :   GEN zknf, czknf, nf, NF, POL;
     300        5586 :   long tx = typ(x);
     301        5586 :   checkrnf(rnf);
     302        5586 :   if (flag) rnfcomplete(rnf);
     303        5586 :   NF = obj_check(rnf,rnf_NFABS);
     304        5586 :   POL = rnf_get_polabs(rnf);
     305        5586 :   if (tx == t_POLMOD && RgX_equal_var(gel(x,1), POL))
     306             :   {
     307          28 :     if (flag) x = nf_to_scalar_or_basis(NF,x);
     308          28 :     return gerepilecopy(av, x);
     309             :   }
     310        5558 :   if (NF && tx == t_COL && lg(x)-1 == nf_get_degree(NF))
     311             :   {
     312           0 :     x = flag? nf_to_scalar_or_basis(NF,x)
     313           0 :             : mkpolmod(nf_to_scalar_or_alg(NF,x), POL);
     314           0 :     return gerepilecopy(av, x);
     315             :   }
     316        5558 :   nf = rnf_get_nf(rnf);
     317        5558 :   if (NF)
     318             :   {
     319             :     GEN d, proj;
     320        5327 :     x = nf_to_scalar_or_basis(nf, x);
     321        5327 :     if (typ(x) != t_COL) return gerepilecopy(av, x);
     322        4571 :     proj = obj_check(rnf,rnf_MAPS);
     323        4571 :     x = Q_remove_denom(x,&d);
     324        4571 :     x = ZM_ZC_mul(gel(proj,1), x);
     325        4571 :     if (d) x = gdiv(x,d);
     326        4571 :     if (!flag) x = basistoalg(NF,x);
     327             :   }
     328             :   else
     329             :   {
     330         231 :     rnf_get_nfzk(rnf, &zknf, &czknf);
     331         231 :     x = nfeltup(nf, x, zknf, czknf);
     332         112 :     if (typ(x) == t_POL) x = mkpolmod(x, POL);
     333             :   }
     334        4683 :   return gerepilecopy(av, x);
     335             : }
     336             : GEN
     337        5320 : rnfeltup(GEN rnf, GEN x) { return rnfeltup0(rnf,x,0); }
     338             : 
     339             : GEN
     340         259 : nfeltup(GEN nf, GEN x, GEN zknf, GEN czknf)
     341             : {
     342             :   GEN c;
     343         259 :   x = nf_to_scalar_or_basis(nf, x);
     344         140 :   if (typ(x) != t_COL) return x;
     345          42 :   x = Q_primitive_part(x, &c);
     346          42 :   if (!RgV_is_ZV(x)) pari_err_TYPE("rnfeltup", x);
     347          42 :   c = mul_content(c, czknf);
     348          42 :   x = RgV_RgC_mul(zknf, x); if (c) x = RgX_Rg_mul(x, c);
     349          42 :   return x;
     350             : }
     351             : 
     352             : static void
     353          49 : fail(const char *f, GEN x)
     354          49 : { pari_err_DOMAIN(f,"element","not in", strtoGENstr("the base field"),x); }
     355             : /* x t_COL of length degabs */
     356             : static GEN
     357           0 : eltdown(GEN rnf, GEN x, long flag)
     358             : {
     359           0 :   GEN z,y, d, proj = obj_check(rnf,rnf_MAPS);
     360           0 :   GEN M= gel(proj,1), iM=gel(proj,2), diM=gel(proj,3), perm=gel(proj,4);
     361           0 :   x = Q_remove_denom(x,&d);
     362           0 :   if (!RgV_is_ZV(x)) pari_err_TYPE("rnfeltdown", x);
     363           0 :   y = ZM_ZC_mul(iM, vecpermute(x, perm));
     364           0 :   z = ZM_ZC_mul(M,y);
     365           0 :   if (!isint1(diM)) z = ZC_Z_mul(z,diM);
     366           0 :   if (!ZV_equal(z,x)) fail("rnfeltdown",x);
     367             : 
     368           0 :   d = mul_denom(d, diM);
     369           0 :   if (d) y = gdiv(y,d);
     370           0 :   if (!flag) y = basistoalg(rnf_get_nf(rnf), y);
     371           0 :   return y;
     372             : }
     373             : GEN
     374        1785 : rnfeltdown0(GEN rnf, GEN x, long flag)
     375             : {
     376        1785 :   const char *f = "rnfeltdown";
     377        1785 :   pari_sp av = avma;
     378             :   GEN z, T, NF, nf;
     379             :   long v;
     380             : 
     381        1785 :   checkrnf(rnf);
     382        1785 :   NF = obj_check(rnf,rnf_NFABS);
     383        1785 :   nf = rnf_get_nf(rnf);
     384        1785 :   T = nf_get_pol(nf);
     385        1785 :   v = varn(T);
     386        1785 :   switch(typ(x))
     387             :   { /* directly belonging to base field ? */
     388         434 :     case t_INT: return icopy(x);
     389          56 :     case t_FRAC:return gcopy(x);
     390             :     case t_POLMOD:
     391        1183 :       if (RgX_equal_var(gel(x,1), rnf_get_polabs(rnf)))
     392             :       {
     393         168 :         if (degpol(T) == 1)
     394             :         {
     395         140 :           x = simplify_shallow(liftpol_shallow(gel(x,2)));
     396         140 :           if (typ(x) != t_POL) return gerepilecopy(av,x);
     397             :         }
     398          35 :         break;
     399             :       }
     400        1015 :       x = polmod_nffix(f,rnf,x,0);
     401             :       /* x was defined mod the relative polynomial & non constant => fail */
     402        1001 :       if (typ(x) == t_POL) fail(f,x);
     403         980 :       if (flag) x = nf_to_scalar_or_basis(nf,x);
     404         980 :       return gerepilecopy(av, x);
     405             : 
     406             :     case t_POL:
     407          63 :       if (varn(x) != v) break;
     408          21 :       x = Rg_nffix(f,T,x,0);
     409          14 :       if (flag) x = nf_to_scalar_or_basis(nf,x);
     410          14 :       return gerepilecopy(av, x);
     411             :     case t_COL:
     412             :     {
     413          49 :       long n = lg(x)-1;
     414          49 :       if (n == degpol(T) && RgV_is_QV(x))
     415             :       {
     416           7 :         if (RgV_isscalar(x)) return gcopy(gel(x,1));
     417           0 :         if (!flag) return gcopy(x);
     418           0 :         return basistoalg(nf,x);
     419             :       }
     420          42 :       if (NF) break;
     421             :     }
     422          42 :     default: pari_err_TYPE(f, x);
     423             :   }
     424             :   /* x defined mod the absolute equation */
     425          77 :   if (NF)
     426             :   {
     427           0 :     x = nf_to_scalar_or_basis(NF, x);
     428           0 :     if (typ(x) == t_COL) x = eltdown(rnf,x,flag);
     429           0 :     return gerepilecopy(av, x);
     430             :   }
     431          77 :   z = rnfeltabstorel(rnf,x);
     432          56 :   switch(typ(z))
     433             :   {
     434             :     case t_INT:
     435          14 :     case t_FRAC: return z;
     436             :   }
     437             :   /* typ(z) = t_POLMOD, varn of both components is rnf_get_varn(rnf) */
     438          42 :   z = gel(z,2);
     439          42 :   if (typ(z) == t_POL)
     440             :   {
     441          42 :     if (lg(z) != 3) fail(f,x);
     442          14 :     z = gel(z,2);
     443             :   }
     444          14 :   return gerepilecopy(av, z);
     445             : }
     446             : GEN
     447        1540 : rnfeltdown(GEN rnf, GEN x) { return rnfeltdown0(rnf,x,0); }
     448             : 
     449             : /* vector of rnf elt -> matrix of nf elts */
     450             : static GEN
     451         476 : rnfV_to_nfM(GEN rnf, GEN x)
     452             : {
     453         476 :   long i, l = lg(x);
     454         476 :   GEN y = cgetg(l, t_MAT);
     455         476 :   for (i = 1; i < l; i++) gel(y,i) = rnfalgtobasis(rnf,gel(x,i));
     456         476 :   return y;
     457             : }
     458             : 
     459             : static GEN
     460         770 : rnfprincipaltohnf(GEN rnf,GEN x)
     461             : {
     462         770 :   pari_sp av = avma;
     463         770 :   GEN bas = rnf_get_zk(rnf), nf = rnf_get_nf(rnf);
     464         770 :   x = rnfbasistoalg(rnf,x);
     465         434 :   x = gmul(x, gmodulo(gel(bas,1), rnf_get_pol(rnf)));
     466         434 :   return gerepileupto(av, nfhnf(nf, mkvec2(rnfV_to_nfM(rnf,x), gel(bas,2))));
     467             : }
     468             : 
     469             : /* pseudo-basis for the 0 ideal */
     470             : static GEN
     471         154 : rnfideal0(void) { retmkvec2(cgetg(1,t_MAT),cgetg(1,t_VEC)); }
     472             : 
     473             : GEN
     474        1316 : rnfidealhnf(GEN rnf, GEN x)
     475             : {
     476             :   GEN z, nf, bas;
     477             : 
     478        1316 :   checkrnf(rnf); nf = rnf_get_nf(rnf);
     479        1316 :   switch(typ(x))
     480             :   {
     481             :     case t_INT: case t_FRAC:
     482         168 :       if (isintzero(x)) return rnfideal0();
     483         112 :       bas = rnf_get_zk(rnf); z = cgetg(3,t_VEC);
     484         112 :       gel(z,1) = matid(rnf_get_degree(rnf));
     485         112 :       gel(z,2) = gmul(x, gel(bas,2)); return z;
     486             : 
     487             :     case t_VEC:
     488         266 :       if (lg(x) == 3 && typ(gel(x,1)) == t_MAT) return nfhnf(nf, x);
     489             :     case t_MAT:
     490         252 :       return rnfidealabstorel(rnf, x);
     491             : 
     492             :     case t_POLMOD: case t_POL: case t_COL:
     493         770 :       return rnfprincipaltohnf(rnf,x);
     494             :   }
     495           0 :   pari_err_TYPE("rnfidealhnf",x);
     496             :   return NULL; /* LCOV_EXCL_LINE */
     497             : }
     498             : 
     499             : static GEN
     500         105 : prodidnorm(GEN nf, GEN I)
     501             : {
     502         105 :   long i, l = lg(I);
     503             :   GEN z;
     504         105 :   if (l == 1) return gen_1;
     505         105 :   z = idealnorm(nf, gel(I,1));
     506         105 :   for (i=2; i<l; i++) z = gmul(z, idealnorm(nf, gel(I,i)));
     507         105 :   return z;
     508             : }
     509             : 
     510             : GEN
     511         196 : rnfidealnormrel(GEN rnf, GEN id)
     512             : {
     513         196 :   pari_sp av = avma;
     514         196 :   GEN nf, z = gel(rnfidealhnf(rnf,id), 2);
     515         126 :   if (lg(z) == 1) return cgetg(1, t_MAT);
     516          98 :   nf = rnf_get_nf(rnf); z = idealprod(nf, z);
     517          98 :   return gerepileupto(av, idealmul(nf,z, rnf_get_index(rnf)));
     518             : }
     519             : 
     520             : GEN
     521         203 : rnfidealnormabs(GEN rnf, GEN id)
     522             : {
     523         203 :   pari_sp av = avma;
     524         203 :   GEN nf, z = gel(rnfidealhnf(rnf,id), 2);
     525         133 :   if (lg(z) == 1) return gen_0;
     526         105 :   nf = rnf_get_nf(rnf); z = prodidnorm(nf, z);
     527         105 :   return gerepileupto(av, gmul(z, gel(rnf,9)));
     528             : }
     529             : 
     530             : static GEN
     531         490 : rnfidealreltoabs_i(GEN rnf, GEN x)
     532             : {
     533             :   long i, l;
     534             :   GEN w;
     535         490 :   x = rnfidealhnf(rnf,x);
     536         350 :   w = gel(x,1); l = lg(w); settyp(w, t_VEC);
     537         350 :   for (i=1; i<l; i++) gel(w,i) = lift_shallow( rnfbasistoalg(rnf, gel(w,i)) );
     538         350 :   return modulereltoabs(rnf, x);
     539             : }
     540             : GEN
     541           0 : rnfidealreltoabs(GEN rnf, GEN x)
     542             : {
     543           0 :   pari_sp av = avma;
     544           0 :   return gerepilecopy(av, rnfidealreltoabs_i(rnf,x));
     545             : }
     546             : GEN
     547         238 : rnfidealreltoabs0(GEN rnf, GEN x, long flag)
     548             : {
     549         238 :   pari_sp av = avma;
     550             :   long i, l;
     551             :   GEN NF;
     552             : 
     553         238 :   x = rnfidealreltoabs_i(rnf, x);
     554         168 :   if (!flag) return gerepilecopy(av,x);
     555          35 :   rnfcomplete(rnf);
     556          35 :   NF = obj_check(rnf,rnf_NFABS);
     557          35 :   l = lg(x); settyp(x, t_MAT);
     558          35 :   for (i=1; i<l; i++) gel(x,i) = algtobasis(NF, gel(x,i));
     559          35 :   return gerepileupto(av, idealhnf(NF,x));
     560             : }
     561             : 
     562             : GEN
     563         455 : rnfidealabstorel(GEN rnf, GEN x)
     564             : {
     565         455 :   long n, N, j, tx = typ(x);
     566         455 :   pari_sp av = avma;
     567             :   GEN A, I, invbas;
     568             : 
     569         455 :   checkrnf(rnf);
     570         455 :   invbas = rnf_get_invzk(rnf);
     571         455 :   if (tx != t_VEC && tx != t_MAT) pari_err_TYPE("rnfidealabstorel",x);
     572         315 :   N = lg(x)-1;
     573         315 :   if (N != rnf_get_absdegree(rnf))
     574             :   {
     575         196 :     if (!N) return rnfideal0();
     576         105 :     pari_err_DIM("rnfidealabstorel");
     577             :   }
     578         119 :   n = rnf_get_degree(rnf);
     579         119 :   A = cgetg(N+1,t_MAT);
     580         119 :   I = cgetg(N+1,t_VEC);
     581         833 :   for (j=1; j<=N; j++)
     582             :   {
     583         714 :     GEN t = lift_shallow( rnfeltabstorel(rnf, gel(x,j)) );
     584         714 :     if (typ(t) == t_POL)
     585         595 :       t = RgM_RgX_mul(invbas, t);
     586             :     else
     587         119 :       t = scalarcol_shallow(t, n);
     588         714 :     gel(A,j) = t;
     589         714 :     gel(I,j) = gen_1;
     590             :   }
     591         119 :   return gerepileupto(av, nfhnf(rnf_get_nf(rnf), mkvec2(A,I)));
     592             : }
     593             : 
     594             : GEN
     595         217 : rnfidealdown(GEN rnf,GEN x)
     596             : {
     597         217 :   pari_sp av = avma;
     598             :   GEN I;
     599         217 :   if (typ(x) == t_MAT)
     600             :   {
     601             :     GEN d;
     602          28 :     x = Q_remove_denom(x,&d);
     603          28 :     if (RgM_is_ZM(x))
     604             :     {
     605          28 :       GEN NF = obj_check(rnf,rnf_NFABS);
     606          28 :       if (NF)
     607             :       {
     608          28 :         GEN z, proj = obj_check(rnf,rnf_MAPS), ZK = gel(proj,1);
     609             :         long i, lz, l;
     610          28 :         x = idealhnf(NF,x);
     611          42 :         if (lg(x) == 1) { avma = av; return cgetg(1,t_MAT); }
     612          14 :         z = ZM_lll(shallowconcat(ZK,x), 0.99, LLL_KER);
     613          14 :         lz = lg(z); l = lg(ZK);
     614          14 :         for (i = 1; i < lz; i++) setlg(gel(z,i), l);
     615          14 :         z = ZM_hnfmodid(z, gcoeff(x,1,1));
     616          14 :         if (d) z = gdiv(z,d);
     617          14 :         return gerepileupto(av, z);
     618             :       }
     619             :     }
     620             :   }
     621         189 :   x = rnfidealhnf(rnf,x); I = gel(x,2);
     622         126 :   if (lg(I) == 1) { avma = av; return cgetg(1,t_MAT); }
     623         105 :   return gerepilecopy(av, gel(I,1));
     624             : }
     625             : 
     626             : /* lift ideal x to the relative extension, returns a Z-basis */
     627             : GEN
     628         217 : rnfidealup(GEN rnf,GEN x)
     629             : {
     630         217 :   pari_sp av = avma;
     631             :   long i, n;
     632             :   GEN nf, bas, bas2, I, x2;
     633             : 
     634         217 :   checkrnf(rnf); nf = rnf_get_nf(rnf);
     635         217 :   n = rnf_get_degree(rnf);
     636         217 :   bas = rnf_get_zk(rnf); bas2 = gel(bas,2);
     637             : 
     638         217 :   (void)idealtyp(&x, &I); /* I is junk */
     639         203 :   x2 = idealtwoelt(nf,x);
     640          98 :   I = cgetg(n+1,t_VEC);
     641         287 :   for (i=1; i<=n; i++)
     642             :   {
     643         189 :     GEN c = gel(bas2,i), d;
     644         189 :     if (typ(c) == t_MAT)
     645             :     {
     646           0 :       c = Q_remove_denom(c,&d);
     647           0 :       c = idealHNF_mul(nf,c,x2);
     648           0 :       if (d) c = gdiv(c,d);
     649             :     }
     650             :     else
     651         189 :       c = idealmul(nf,c,x);
     652         189 :     gel(I,i) = c;
     653             :   }
     654          98 :   return gerepilecopy(av, modulereltoabs(rnf, mkvec2(gel(bas,1), I)));
     655             : }
     656             : GEN
     657        1043 : rnfidealup0(GEN rnf,GEN x, long flag)
     658             : {
     659        1043 :   pari_sp av = avma;
     660             :   GEN NF, nf, proj, d, x2;
     661             : 
     662        1043 :   if (!flag) return rnfidealup(rnf,x);
     663         826 :   checkrnf(rnf); nf = rnf_get_nf(rnf);
     664         826 :   rnfcomplete(rnf);
     665         826 :   proj = obj_check(rnf,rnf_MAPS);
     666         826 :   NF = obj_check(rnf,rnf_NFABS);
     667             : 
     668         826 :   (void)idealtyp(&x, &d); /* d is junk */
     669         826 :   x2 = idealtwoelt(nf,x);
     670         826 :   x2 = Q_remove_denom(x2,&d);
     671         826 :   gel(x2,2) = ZM_ZC_mul(gel(proj,1),gel(x2,2));
     672         826 :   x2 = idealhnf_two(NF, x2);
     673         826 :   if (d) x2 = gdiv(x2,d);
     674         826 :   return gerepileupto(av, x2);
     675             : }
     676             : 
     677             : /* x a relative HNF => vector of 2 generators (relative polmods) */
     678             : GEN
     679         252 : rnfidealtwoelement(GEN rnf, GEN x)
     680             : {
     681         252 :   pari_sp av = avma;
     682             :   GEN y, cy, z, NF;
     683             : 
     684         252 :   y = rnfidealreltoabs_i(rnf,x);
     685         182 :   rnfcomplete(rnf);
     686         182 :   NF = obj_check(rnf,rnf_NFABS);
     687         182 :   y = matalgtobasis(NF, y); settyp(y, t_MAT);
     688         182 :   y = Q_primitive_part(y, &cy);
     689         182 :   y = ZM_hnf(y);
     690         182 :   if (lg(y) == 1) { avma = av; return mkvec2(gen_0, gen_0); }
     691         147 :   y = idealtwoelt(NF, y);
     692         140 :   if (cy) y = RgV_Rg_mul(y, cy);
     693         140 :   z = gel(y,2);
     694         140 :   if (typ(z) == t_COL) z = rnfeltabstorel(rnf, nf_to_scalar_or_alg(NF, z));
     695         140 :   return gerepilecopy(av, mkvec2(gel(y,1), z));
     696             : }
     697             : 
     698             : GEN
     699          49 : rnfidealmul(GEN rnf,GEN x,GEN y)
     700             : {
     701          49 :   pari_sp av = avma;
     702             :   GEN nf, z, x1, x2, p1, p2, bas;
     703             : 
     704          49 :   y = rnfidealtwoelement(rnf,y);
     705          49 :   if (isintzero(gel(y,1))) { avma = av; return rnfideal0(); }
     706          42 :   nf = rnf_get_nf(rnf);
     707          42 :   bas = rnf_get_zk(rnf);
     708          42 :   x = rnfidealhnf(rnf,x);
     709          42 :   x1 = gmodulo(gmul(gel(bas,1), matbasistoalg(nf,gel(x,1))), rnf_get_pol(rnf));
     710          42 :   x2 = gel(x,2);
     711          42 :   p1 = gmul(gel(y,1), gel(x,1));
     712          42 :   p2 = rnfV_to_nfM(rnf, gmul(gel(y,2), x1));
     713          42 :   z = mkvec2(shallowconcat(p1, p2), shallowconcat(x2, x2));
     714          42 :   return gerepileupto(av, nfhnf(nf,z));
     715             : }
     716             : 
     717             : static GEN
     718          35 : rnfidealprimedec_1(GEN rnf, GEN SL, GEN prK)
     719             : {
     720          35 :   GEN v, piL = rnfeltup0(rnf, pr_get_gen(prK), 1);
     721             :   long i, c, l;
     722          35 :   if (typ(piL) != t_COL) return SL; /* p inert in K/Q */
     723          35 :   v = cgetg_copy(SL, &l);
     724          70 :   for (i = c = 1; i < l; i++)
     725             :   {
     726          35 :     GEN P = gel(SL,i);
     727          35 :     if (ZC_prdvd(piL, P)) gel(v,c++) = P;
     728             :   }
     729          35 :   setlg(v, c); return v;
     730             : }
     731             : GEN
     732          35 : rnfidealprimedec(GEN rnf, GEN pr)
     733             : {
     734          35 :   pari_sp av = avma;
     735             :   GEN p, z, NF, nf, SL;
     736          35 :   checkrnf(rnf);
     737          35 :   rnfcomplete(rnf);
     738          35 :   NF = obj_check(rnf,rnf_NFABS);
     739          35 :   nf = rnf_get_nf(rnf);
     740          35 :   if (typ(pr) == t_INT)
     741             :   {
     742          28 :     p = pr;
     743          28 :     pr = NULL;
     744             :   }
     745             :   else
     746             :   {
     747           7 :     checkprid(pr);
     748           7 :     p = pr_get_p(pr);
     749             :   }
     750          35 :   SL = idealprimedec(NF, p);
     751          35 :   if (pr) z = rnfidealprimedec_1(rnf, SL, pr);
     752             :   else
     753             :   {
     754          28 :     GEN vK = idealprimedec(nf, p), vL;
     755          28 :     long l = lg(vK), i;
     756          28 :     vL = cgetg(l, t_VEC);
     757          28 :     for (i = 1; i < l; i++) gel(vL,i) = rnfidealprimedec_1(rnf, SL, gel(vK,i));
     758          28 :     z = mkvec2(vK, vL);
     759             :   }
     760          35 :   return gerepilecopy(av, z);
     761             : }
     762             : 
     763             : GEN
     764          35 : rnfidealfactor(GEN rnf, GEN x)
     765             : {
     766          35 :   pari_sp av = avma;
     767             :   GEN NF;
     768          35 :   checkrnf(rnf);
     769          35 :   rnfcomplete(rnf);
     770          35 :   NF = obj_check(rnf,rnf_NFABS);
     771          35 :   return gerepileupto(av, idealfactor(NF, rnfidealreltoabs0(rnf, x, 1)));
     772             : }
     773             : 
     774             : GEN
     775        1813 : rnfequationall(GEN A, GEN B, long *pk, GEN *pLPRS)
     776             : {
     777             :   long lA, lB;
     778             :   GEN nf, C;
     779             : 
     780        1813 :   A = get_nfpol(A, &nf); lA = lg(A);
     781        1813 :   if (!nf) {
     782         588 :     if (lA<=3) pari_err_CONSTPOL("rnfequation");
     783         588 :     RgX_check_ZX(A,"rnfequation");
     784             :   }
     785        1813 :   B = RgX_nffix("rnfequation", A,B,1); lB = lg(B);
     786        1813 :   if (lB<=3) pari_err_CONSTPOL("rnfequation");
     787        1813 :   B = Q_primpart(B);
     788             : 
     789        1813 :   if (!nfissquarefree(A,B))
     790           0 :     pari_err_DOMAIN("rnfequation","issquarefree(B)","=",gen_0,B);
     791             : 
     792        1813 :   *pk = 0; C = ZX_ZXY_resultant_all(A, B, pk, pLPRS);
     793        1813 :   if (gsigne(leading_coeff(C)) < 0) C = RgX_neg(C);
     794        1813 :   *pk = -*pk; return Q_primpart(C);
     795             : }
     796             : 
     797             : GEN
     798        1757 : rnfequation0(GEN A, GEN B, long flall)
     799             : {
     800        1757 :   pari_sp av = avma;
     801             :   GEN LPRS, C;
     802             :   long k;
     803             : 
     804        1757 :   C = rnfequationall(A, B, &k, flall? &LPRS: NULL);
     805        1757 :   if (flall)
     806             :   { /* a,b,c root of A,B,C = compositum, c = b + k a */
     807        1323 :     GEN a, mH0 = RgX_neg(gel(LPRS,1)), H1 = gel(LPRS,2);
     808        1323 :     a = RgXQ_mul(mH0, QXQ_inv(H1, C), C);
     809        1323 :     C = mkvec3(C, mkpolmod(a, C), stoi(k));
     810             :   }
     811        1757 :   return gerepilecopy(av, C);
     812             : }
     813             : GEN
     814         420 : rnfequation(GEN nf, GEN pol) { return rnfequation0(nf,pol,0); }
     815             : GEN
     816        1211 : rnfequation2(GEN nf, GEN pol) { return rnfequation0(nf,pol,1); }
     817             : GEN
     818        1197 : nf_rnfeq(GEN nf, GEN relpol)
     819             : {
     820             :   GEN pol, a, k, junk, eq;
     821        1197 :   relpol = liftpol_shallow(relpol);
     822        1197 :   eq = rnfequation2(nf, relpol);
     823        1197 :   pol = gel(eq,1);
     824        1197 :   a = gel(eq,2); if (typ(a) == t_POLMOD) a = gel(a,2);
     825        1197 :   k = gel(eq,3);
     826        1197 :   return mkvec5(pol,a,k,get_nfpol(nf, &junk),relpol);
     827             : }
     828             : /* only allow abstorel */
     829             : GEN
     830          21 : nf_rnfeqsimple(GEN nf, GEN relpol)
     831             : {
     832             :   long sa;
     833             :   GEN junk, pol;
     834          21 :   relpol = liftpol_shallow(relpol);
     835          21 :   pol = rnfequationall(nf, relpol, &sa, NULL);
     836          21 :   return mkvec5(pol,gen_0/*dummy*/,stoi(sa),get_nfpol(nf, &junk),relpol);
     837             : }
     838             : 
     839             : /*******************************************************************/
     840             : /*                                                                 */
     841             : /*                            RELATIVE LLL                         */
     842             : /*                                                                 */
     843             : /*******************************************************************/
     844             : static GEN
     845         196 : nftau(long r1, GEN x)
     846             : {
     847         196 :   long i, l = lg(x);
     848         196 :   GEN s = r1? gel(x,1): gmul2n(real_i(gel(x,1)),1);
     849         196 :   for (i=2; i<=r1; i++) s = gadd(s, gel(x,i));
     850         196 :   for (   ; i < l; i++) s = gadd(s, gmul2n(real_i(gel(x,i)),1));
     851         196 :   return s;
     852             : }
     853             : 
     854             : static GEN
     855          28 : initmat(long l)
     856             : {
     857          28 :   GEN x = cgetg(l, t_MAT);
     858             :   long i;
     859          28 :   for (i = 1; i < l; i++) gel(x,i) = cgetg(l, t_COL);
     860          28 :   return x;
     861             : }
     862             : 
     863             : static GEN
     864        1022 : nftocomplex(GEN nf, GEN x)
     865             : {
     866        1022 :   GEN M = nf_get_M(nf);
     867        1022 :   x = nf_to_scalar_or_basis(nf,x);
     868        1022 :   if (typ(x) != t_COL) return const_col(nbrows(M), x);
     869         161 :   return RgM_RgC_mul(M, x);
     870             : }
     871             : /* assume x a square t_MAT, return a t_VEC of embeddings of its columns */
     872             : static GEN
     873          14 : mattocomplex(GEN nf, GEN x)
     874             : {
     875          14 :   long i,j, l = lg(x);
     876          14 :   GEN v = cgetg(l, t_VEC);
     877          98 :   for (j=1; j<l; j++)
     878             :   {
     879          84 :     GEN c = gel(x,j), b = cgetg(l, t_MAT);
     880          84 :     for (i=1; i<l; i++) gel(b,i) = nftocomplex(nf, gel(c,i));
     881          84 :     b = shallowtrans(b); settyp(b, t_COL);
     882          84 :     gel(v,j) = b;
     883             :   }
     884          14 :   return v;
     885             : }
     886             : 
     887             : static GEN
     888          14 : nf_all_roots(GEN nf, GEN x, long prec)
     889             : {
     890          14 :   long i, j, l = lg(x), ru = lg(nf_get_roots(nf));
     891          14 :   GEN y = cgetg(l, t_POL), v, z;
     892             : 
     893          14 :   x = RgX_to_nfX(nf, x);
     894          14 :   y[1] = x[1];
     895          14 :   for (i=2; i<l; i++) gel(y,i) = nftocomplex(nf, gel(x,i));
     896          14 :   i = gprecision(y); if (i && i <= 3) return NULL;
     897             : 
     898          14 :   v = cgetg(ru, t_VEC);
     899          14 :   z = cgetg(l, t_POL); z[1] = x[1];
     900          42 :   for (i=1; i<ru; i++)
     901             :   {
     902          28 :     for (j = 2; j < l; j++) gel(z,j) = gmael(y,j,i);
     903          28 :     gel(v,i) = cleanroots(z, prec);
     904             :   }
     905          14 :   return v;
     906             : }
     907             : 
     908             : static GEN
     909         357 : rnfscal(GEN m, GEN x, GEN y)
     910             : {
     911         357 :   long i, l = lg(m);
     912         357 :   GEN z = cgetg(l, t_COL);
     913        1071 :   for (i = 1; i < l; i++)
     914         714 :     gel(z,i) = gmul(gconj(shallowtrans(gel(x,i))), gmul(gel(m,i), gel(y,i)));
     915         357 :   return z;
     916             : }
     917             : 
     918             : /* x ideal in HNF */
     919             : static GEN
     920         364 : findmin(GEN nf, GEN x, GEN muf)
     921             : {
     922         364 :   pari_sp av = avma;
     923             :   long e;
     924         364 :   GEN cx, y, m, M = nf_get_M(nf);
     925             : 
     926         364 :   x = Q_primitive_part(x, &cx);
     927         364 :   if (gequal1(gcoeff(x,1,1))) y = M;
     928             :   else
     929             :   {
     930         210 :     GEN G = nf_get_G(nf);
     931         210 :     m = lllfp(RgM_mul(G,x), 0.75, 0);
     932         210 :     if (typ(m) != t_MAT)
     933             :     {
     934           0 :       x = ZM_lll(x, 0.75, LLL_INPLACE);
     935           0 :       m = lllfp(RgM_mul(G,x), 0.75, 0);
     936           0 :       if (typ(m) != t_MAT) pari_err_PREC("rnflllgram");
     937             :     }
     938         210 :     x = ZM_mul(x, m);
     939         210 :     y = RgM_mul(M, x);
     940             :   }
     941         364 :   m = RgM_solve_realimag(y, muf);
     942         364 :   if (!m) return NULL; /* precision problem */
     943         364 :   if (cx) m = RgC_Rg_div(m, cx);
     944         364 :   m = grndtoi(m, &e);
     945         364 :   if (e >= 0) return NULL; /* precision problem */
     946         364 :   m = ZM_ZC_mul(x, m);
     947         364 :   if (cx) m = RgC_Rg_mul(m, cx);
     948         364 :   return gerepileupto(av, m);
     949             : }
     950             : 
     951             : static int
     952         364 : RED(long k, long l, GEN U, GEN mu, GEN MC, GEN nf, GEN I, GEN *Ik_inv)
     953             : {
     954             :   GEN x, xc, ideal;
     955             :   long i;
     956             : 
     957         364 :   if (!*Ik_inv) *Ik_inv = idealinv(nf, gel(I,k));
     958         364 :   ideal = idealmul(nf,gel(I,l), *Ik_inv);
     959         364 :   x = findmin(nf, ideal, gcoeff(mu,k,l));
     960         364 :   if (!x) return 0;
     961         364 :   if (gequal0(x)) return 1;
     962             : 
     963         294 :   xc = nftocomplex(nf,x);
     964         294 :   gel(MC,k) = gsub(gel(MC,k), vecmul(xc,gel(MC,l)));
     965         294 :   gel(U,k) = gsub(gel(U,k), gmul(coltoalg(nf,x), gel(U,l)));
     966         294 :   gcoeff(mu,k,l) = gsub(gcoeff(mu,k,l), xc);
     967        1029 :   for (i=1; i<l; i++)
     968         735 :     gcoeff(mu,k,i) = gsub(gcoeff(mu,k,i), vecmul(xc,gcoeff(mu,l,i)));
     969         294 :   return 1;
     970             : }
     971             : 
     972             : static int
     973          84 : check_0(GEN B)
     974             : {
     975          84 :   long i, l = lg(B);
     976         252 :   for (i = 1; i < l; i++)
     977         168 :     if (gsigne(gel(B,i)) <= 0) return 1;
     978          84 :   return 0;
     979             : }
     980             : 
     981             : static int
     982          98 : do_SWAP(GEN I, GEN MC, GEN MCS, GEN h, GEN mu, GEN B, long kmax, long k,
     983             :         const long alpha, long r1)
     984             : {
     985             :   GEN p1, p2, muf, mufc, Bf, temp;
     986             :   long i, j;
     987             : 
     988         196 :   p1 = nftau(r1, gadd(gel(B,k),
     989         196 :                       gmul(gnorml2(gcoeff(mu,k,k-1)), gel(B,k-1))));
     990          98 :   p2 = nftau(r1, gel(B,k-1));
     991          98 :   if (gcmp(gmulsg(alpha,p1), gmulsg(alpha-1,p2)) > 0) return 0;
     992             : 
     993          14 :   swap(gel(MC,k-1),gel(MC,k));
     994          14 :   swap(gel(h,k-1), gel(h,k));
     995          14 :   swap(gel(I,k-1), gel(I,k));
     996          14 :   for (j=1; j<=k-2; j++) swap(gcoeff(mu,k-1,j),gcoeff(mu,k,j));
     997          14 :   muf = gcoeff(mu,k,k-1);
     998          14 :   mufc = gconj(muf);
     999          14 :   Bf = gadd(gel(B,k), vecmul(real_i(vecmul(muf,mufc)), gel(B,k-1)));
    1000          14 :   if (check_0(Bf)) return 1; /* precision problem */
    1001             : 
    1002          14 :   p1 = vecdiv(gel(B,k-1),Bf);
    1003          14 :   gcoeff(mu,k,k-1) = vecmul(mufc,p1);
    1004          14 :   temp = gel(MCS,k-1);
    1005          14 :   gel(MCS,k-1) = gadd(gel(MCS,k), vecmul(muf,gel(MCS,k-1)));
    1006          42 :   gel(MCS,k) = gsub(vecmul(vecdiv(gel(B,k),Bf), temp),
    1007          28 :                     vecmul(gcoeff(mu,k,k-1), gel(MCS,k)));
    1008          14 :   gel(B,k) = vecmul(gel(B,k),p1);
    1009          14 :   gel(B,k-1) = Bf;
    1010          14 :   for (i=k+1; i<=kmax; i++)
    1011             :   {
    1012           0 :     temp = gcoeff(mu,i,k);
    1013           0 :     gcoeff(mu,i,k) = gsub(gcoeff(mu,i,k-1), vecmul(muf, gcoeff(mu,i,k)));
    1014           0 :     gcoeff(mu,i,k-1) = gadd(temp, vecmul(gcoeff(mu,k,k-1),gcoeff(mu,i,k)));
    1015             :   }
    1016          14 :   return 1;
    1017             : }
    1018             : 
    1019             : static GEN
    1020          14 : rel_T2(GEN nf, GEN pol, long lx, long prec)
    1021             : {
    1022             :   long ru, i, j, k, l;
    1023             :   GEN T2, s, unro, roorder, powreorder;
    1024             : 
    1025          14 :   roorder = nf_all_roots(nf, pol, prec);
    1026          14 :   if (!roorder) return NULL;
    1027          14 :   ru = lg(roorder);
    1028          14 :   unro = cgetg(lx,t_COL); for (i=1; i<lx; i++) gel(unro,i) = gen_1;
    1029          14 :   powreorder = cgetg(lx,t_MAT); gel(powreorder,1) = unro;
    1030          14 :   T2 = cgetg(ru, t_VEC);
    1031          42 :   for (i = 1; i < ru; i++)
    1032             :   {
    1033          28 :     GEN ro = gel(roorder,i);
    1034          28 :     GEN m = initmat(lx);
    1035         168 :     for (k=2; k<lx; k++)
    1036             :     {
    1037         140 :       GEN c = cgetg(lx, t_COL); gel(powreorder,k) = c;
    1038        1232 :       for (j=1; j < lx; j++)
    1039        1092 :         gel(c,j) = gmul(gel(ro,j), gmael(powreorder,k-1,j));
    1040             :     }
    1041         196 :     for (l = 1; l < lx; l++)
    1042         882 :       for (k = 1; k <= l; k++)
    1043             :       {
    1044         714 :         s = gen_0;
    1045        6636 :         for (j = 1; j < lx; j++)
    1046        5922 :           s = gadd(s, gmul(gconj(gmael(powreorder,k,j)),
    1047        5922 :                                  gmael(powreorder,l,j)));
    1048         714 :         if (l == k)
    1049         168 :           gcoeff(m, l, l) = real_i(s);
    1050             :         else
    1051             :         {
    1052         546 :           gcoeff(m, k, l) = s;
    1053         546 :           gcoeff(m, l, k) = gconj(s);
    1054             :         }
    1055             :       }
    1056          28 :     gel(T2,i) = m;
    1057             :   }
    1058          14 :   return T2;
    1059             : }
    1060             : 
    1061             : /* given a base field nf (e.g main variable y), a polynomial pol with
    1062             :  * coefficients in nf    (e.g main variable x), and an order as output
    1063             :  * by rnfpseudobasis, outputs a reduced order. */
    1064             : GEN
    1065          14 : rnflllgram(GEN nf, GEN pol, GEN order,long prec)
    1066             : {
    1067          14 :   pari_sp av = avma;
    1068          14 :   long j, k, l, kmax, r1, lx, count = 0;
    1069             :   GEN M, I, h, H, mth, MC, MPOL, MCS, B, mu;
    1070          14 :   const long alpha = 10, MAX_COUNT = 4;
    1071             : 
    1072          14 :   nf = checknf(nf); r1 = nf_get_r1(nf);
    1073          14 :   check_ZKmodule(order, "rnflllgram");
    1074          14 :   M = gel(order,1);
    1075          14 :   I = gel(order,2); lx = lg(I);
    1076          14 :   if (lx < 3) return gcopy(order);
    1077          14 :   if (lx-1 != degpol(pol)) pari_err_DIM("rnflllgram");
    1078          14 :   I = leafcopy(I);
    1079          14 :   H = NULL;
    1080          14 :   MPOL = matbasistoalg(nf, M);
    1081          14 :   MCS = matid(lx-1); /* dummy for gerepile */
    1082             : PRECNF:
    1083          14 :   if (count == MAX_COUNT)
    1084             :   {
    1085           0 :     prec = precdbl(prec); count = 0;
    1086           0 :     if (DEBUGLEVEL) pari_warn(warnprec,"rnflllgram",prec);
    1087           0 :     nf = nfnewprec_shallow(nf,prec);
    1088             :   }
    1089          14 :   mth = rel_T2(nf, pol, lx, prec);
    1090          14 :   if (!mth) { count = MAX_COUNT; goto PRECNF; }
    1091          14 :   h = NULL;
    1092             : PRECPB:
    1093          14 :   if (h)
    1094             :   { /* precision problem, recompute. If no progress, increase nf precision */
    1095           0 :     if (++count == MAX_COUNT || RgM_isidentity(h)) {count = MAX_COUNT; goto PRECNF;}
    1096           0 :     H = H? gmul(H, h): h;
    1097           0 :     MPOL = gmul(MPOL, h);
    1098             :   }
    1099          14 :   h = matid(lx-1);
    1100          14 :   MC = mattocomplex(nf, MPOL);
    1101          14 :   mu = cgetg(lx,t_MAT);
    1102          14 :   B  = cgetg(lx,t_COL);
    1103          98 :   for (j=1; j<lx; j++)
    1104             :   {
    1105          84 :     gel(mu,j) = zerocol(lx - 1);
    1106          84 :     gel(B,j) = gen_0;
    1107             :   }
    1108          14 :   if (DEBUGLEVEL) err_printf("k = ");
    1109          14 :   gel(B,1) = real_i(rnfscal(mth,gel(MC,1),gel(MC,1)));
    1110          14 :   gel(MCS,1) = gel(MC,1);
    1111          14 :   kmax = 1; k = 2;
    1112             :   do
    1113             :   {
    1114          98 :     GEN Ik_inv = NULL;
    1115          98 :     if (DEBUGLEVEL) err_printf("%ld ",k);
    1116          98 :     if (k > kmax)
    1117             :     { /* Incremental Gram-Schmidt */
    1118          70 :       kmax = k; gel(MCS,k) = gel(MC,k);
    1119         343 :       for (j=1; j<k; j++)
    1120             :       {
    1121         546 :         gcoeff(mu,k,j) = vecdiv(rnfscal(mth,gel(MCS,j),gel(MC,k)),
    1122         273 :                                 gel(B,j));
    1123         273 :         gel(MCS,k) = gsub(gel(MCS,k), vecmul(gcoeff(mu,k,j),gel(MCS,j)));
    1124             :       }
    1125          70 :       gel(B,k) = real_i(rnfscal(mth,gel(MCS,k),gel(MCS,k)));
    1126          70 :       if (check_0(gel(B,k))) goto PRECPB;
    1127             :     }
    1128          98 :     if (!RED(k, k-1, h, mu, MC, nf, I, &Ik_inv)) goto PRECPB;
    1129          98 :     if (do_SWAP(I,MC,MCS,h,mu,B,kmax,k,alpha, r1))
    1130             :     {
    1131          14 :       if (!B[k]) goto PRECPB;
    1132          14 :       if (k > 2) k--;
    1133             :     }
    1134             :     else
    1135             :     {
    1136         350 :       for (l=k-2; l; l--)
    1137         266 :         if (!RED(k, l, h, mu, MC, nf, I, &Ik_inv)) goto PRECPB;
    1138          84 :       k++;
    1139             :     }
    1140          98 :     if (gc_needed(av,2))
    1141             :     {
    1142           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"rnflllgram");
    1143           0 :       gerepileall(av, H?10:9, &nf,&mth,&h,&MPOL,&B,&MC,&MCS,&mu,&I,&H);
    1144             :     }
    1145             :   }
    1146          98 :   while (k < lx);
    1147          14 :   MPOL = gmul(MPOL,h);
    1148          14 :   if (H) h = gmul(H, h);
    1149          14 :   if (DEBUGLEVEL) err_printf("\n");
    1150          14 :   MPOL = RgM_to_nfM(nf,MPOL);
    1151          14 :   h = RgM_to_nfM(nf,h);
    1152          14 :   return gerepilecopy(av, mkvec2(mkvec2(MPOL,I), h));
    1153             : }
    1154             : 
    1155             : GEN
    1156           7 : rnfpolred(GEN nf, GEN pol, long prec)
    1157             : {
    1158           7 :   pari_sp av = avma;
    1159           7 :   long i, j, n, v = varn(pol);
    1160             :   GEN id, w, I, O, bnf, nfpol;
    1161             : 
    1162           7 :   if (typ(pol)!=t_POL) pari_err_TYPE("rnfpolred",pol);
    1163           7 :   bnf = nf; nf = checknf(bnf);
    1164           7 :   bnf = (nf == bnf)? NULL: checkbnf(bnf);
    1165           7 :   if (degpol(pol) <= 1) { w = cgetg(2, t_VEC); gel(w,1) = pol_x(v); return w; }
    1166           7 :   nfpol = nf_get_pol(nf);
    1167             : 
    1168           7 :   id = rnfpseudobasis(nf,pol);
    1169           7 :   if (bnf && is_pm1( bnf_get_no(bnf) )) /* if bnf is principal */
    1170             :   {
    1171             :     GEN newI, newO;
    1172           0 :     O = gel(id,1);
    1173           0 :     I = gel(id,2); n = lg(I)-1;
    1174           0 :     newI = cgetg(n+1,t_VEC);
    1175           0 :     newO = cgetg(n+1,t_MAT);
    1176           0 :     for (j=1; j<=n; j++)
    1177             :     {
    1178           0 :       GEN al = gen_if_principal(bnf,gel(I,j));
    1179           0 :       gel(newI,j) = gen_1;
    1180           0 :       gel(newO,j) = nfC_nf_mul(nf, gel(O,j), al);
    1181             :     }
    1182           0 :     id = mkvec2(newO, newI);
    1183             :   }
    1184             : 
    1185           7 :   id = gel(rnflllgram(nf,pol,id,prec),1);
    1186           7 :   O = gel(id,1);
    1187           7 :   I = gel(id,2); n = lg(I)-1;
    1188           7 :   w = cgetg(n+1,t_VEC);
    1189           7 :   pol = lift_shallow(pol);
    1190          70 :   for (j=1; j<=n; j++)
    1191             :   {
    1192          63 :     GEN newpol, L, a, Ij = gel(I,j);
    1193          63 :     a = RgC_Rg_mul(gel(O,j), (typ(Ij) == t_MAT)? gcoeff(Ij,1,1): Ij);
    1194          63 :     for (i=n; i; i--) gel(a,i) = nf_to_scalar_or_alg(nf, gel(a,i));
    1195          63 :     a = RgV_to_RgX(a, v);
    1196          63 :     newpol = RgXQX_red(RgXQ_charpoly(a, pol, v), nfpol);
    1197          63 :     newpol = Q_primpart(newpol);
    1198             : 
    1199          63 :     (void)nfgcd_all(newpol, RgX_deriv(newpol), nfpol, nf_get_index(nf), &newpol);
    1200          63 :     L = leading_coeff(newpol);
    1201         126 :     gel(w,j) = (typ(L) == t_POL)? RgXQX_div(newpol, L, nfpol)
    1202          63 :                                 : RgX_Rg_div(newpol, L);
    1203             :   }
    1204           7 :   return gerepilecopy(av,w);
    1205             : }
    1206             : 
    1207             : /*******************************************************************/
    1208             : /*                                                                 */
    1209             : /*                  LINEAR ALGEBRA OVER Z_K  (HNF,SNF)             */
    1210             : /*                                                                 */
    1211             : /*******************************************************************/
    1212             : /* A torsion-free module M over Z_K is given by [A,I].
    1213             :  * I=[a_1,...,a_k] is a row vector of k fractional ideals given in HNF.
    1214             :  * A is an n x k matrix (same k) such that if A_j is the j-th column of A then
    1215             :  * M=a_1 A_1+...+a_k A_k. We say that [A,I] is a pseudo-basis if k=n */
    1216             : 
    1217             : /* Given an element x and an ideal I in HNF, gives an r such that x-r is in H
    1218             :  * and r is small */
    1219             : GEN
    1220           7 : nfreduce(GEN nf, GEN x, GEN I)
    1221             : {
    1222           7 :   pari_sp av = avma;
    1223             :   GEN aI;
    1224           7 :   x = nf_to_scalar_or_basis(checknf(nf), x);
    1225           7 :   if (idealtyp(&I,&aI) != id_MAT || lg(I)==1) pari_err_TYPE("nfreduce",I);
    1226           7 :   if (typ(x) != t_COL) x = scalarcol( gmod(x, gcoeff(I,1,1)), lg(I)-1 );
    1227           7 :   else x = reducemodinvertible(x, I);
    1228           7 :   return gerepileupto(av, x);
    1229             : }
    1230             : /* Given an element x and an ideal in HNF, gives an a in ideal such that
    1231             :  * x-a is small. No checks */
    1232             : static GEN
    1233       30520 : element_close(GEN nf, GEN x, GEN ideal)
    1234             : {
    1235       30520 :   pari_sp av = avma;
    1236       30520 :   GEN y = gcoeff(ideal,1,1);
    1237       30520 :   x = nf_to_scalar_or_basis(nf, x);
    1238       30520 :   if (typ(y) == t_INT && is_pm1(y)) return ground(x);
    1239       29295 :   if (typ(x) == t_COL)
    1240       12964 :     x = closemodinvertible(x, ideal);
    1241             :   else
    1242       16331 :     x = gmul(y, gdivround(x,y));
    1243       29295 :   return gerepileupto(av, x);
    1244             : }
    1245             : 
    1246             : /* A + v B */
    1247             : static GEN
    1248      115528 : colcomb1(GEN nf, GEN v, GEN A, GEN B)
    1249             : {
    1250      115528 :   if (isintzero(v)) return A;
    1251       75866 :   return RgC_to_nfC(nf, RgC_add(A, nfC_nf_mul(nf,B,v)));
    1252             : }
    1253             : /* u A + v B */
    1254             : static GEN
    1255       98882 : colcomb(GEN nf, GEN u, GEN v, GEN A, GEN B)
    1256             : {
    1257       98882 :   if (isintzero(u)) return nfC_nf_mul(nf,B,v);
    1258       84987 :   if (u != gen_1) A = nfC_nf_mul(nf,A,u);
    1259       84987 :   return colcomb1(nf, v, A, B);
    1260             : }
    1261             : 
    1262             : /* return m[i,1..lim] * x */
    1263             : static GEN
    1264         231 : element_mulvecrow(GEN nf, GEN x, GEN m, long i, long lim)
    1265             : {
    1266         231 :   long j, l = minss(lg(m), lim+1);
    1267         231 :   GEN dx, y = cgetg(l, t_VEC);
    1268         231 :   x = nf_to_scalar_or_basis(nf, x);
    1269         231 :   if (typ(x) == t_COL)
    1270             :   {
    1271          91 :     x = zk_multable(nf, Q_remove_denom(x, &dx));
    1272         350 :     for (j=1; j<l; j++)
    1273             :     {
    1274         259 :       GEN t = gcoeff(m,i,j);
    1275         259 :       if (!isintzero(t))
    1276             :       {
    1277         112 :         if (typ(t) == t_COL)
    1278          28 :           t = RgM_RgC_mul(x, t);
    1279             :         else
    1280          84 :           t = RgC_Rg_mul(gel(x,1), t);
    1281         112 :         if (dx) t = gdiv(t, dx);
    1282         112 :         t = nf_to_scalar_or_basis(nf,t);
    1283             :       }
    1284         259 :       gel(y,j) = t;
    1285             :     }
    1286             :   }
    1287             :   else
    1288             :   {
    1289         140 :     for (j=1; j<l; j++) gel(y,j) = gmul(x, gcoeff(m,i,j));
    1290             :   }
    1291         231 :   return y;
    1292             : }
    1293             : 
    1294             : /* u Z[s,] + v Z[t,], limitied to the first lim entries */
    1295             : static GEN
    1296         154 : rowcomb(GEN nf, GEN u, GEN v, long s, long t, GEN Z, long lim)
    1297             : {
    1298             :   GEN z;
    1299         154 :   if (gequal0(u))
    1300           7 :     z = element_mulvecrow(nf,v,Z,t, lim);
    1301             :   else
    1302             :   {
    1303         147 :     z = element_mulvecrow(nf,u,Z,s, lim);
    1304         147 :     if (!gequal0(v)) z = gadd(z, element_mulvecrow(nf,v,Z,t, lim));
    1305             :   }
    1306         154 :   return z;
    1307             : }
    1308             : 
    1309             : /* nfbezout(0,b,A,B). Either bB = NULL or b*B */
    1310             : static GEN
    1311       54229 : zero_nfbezout(GEN nf,GEN bB, GEN b, GEN A,GEN B,GEN *u,GEN *v,GEN *w,GEN *di)
    1312             : {
    1313             :   GEN d;
    1314       54229 :   if (isint1(b))
    1315             :   {
    1316       52822 :     *v = gen_1;
    1317       52822 :     *w = A;
    1318       52822 :     d = B;
    1319       52822 :     *di = idealinv(nf,d);
    1320             :   }
    1321             :   else
    1322             :   {
    1323        1407 :     *v = nfinv(nf,b);
    1324        1407 :     *w = idealmul(nf,A,*v);
    1325        1407 :     d = bB? bB: idealmul(nf,b,B);
    1326        1407 :     *di = idealHNF_inv(nf,d);
    1327             :   }
    1328       54229 :   *u = gen_0; return d;
    1329             : }
    1330             : 
    1331             : /* Given elements a,b and ideals A, B, outputs d = a.A+b.B and gives
    1332             :  * di=d^-1, w=A.B.di, u, v such that au+bv=1 and u in A.di, v in B.di.
    1333             :  * Assume A, B non-zero, but a or b can be zero (not both) */
    1334             : static GEN
    1335       58254 : nfbezout(GEN nf,GEN a,GEN b, GEN A,GEN B, GEN *pu,GEN *pv,GEN *pw,GEN *pdi)
    1336             : {
    1337             :   GEN w, u, v, d, di, aA, bB;
    1338             : 
    1339       58254 :   if (isintzero(a)) return zero_nfbezout(nf,NULL,b,A,B,pu,pv,pw,pdi);
    1340       58254 :   if (isintzero(b)) return zero_nfbezout(nf,NULL,a,B,A,pv,pu,pw,pdi);
    1341             : 
    1342       58254 :   if (a != gen_1) /* frequently called with a = gen_1 */
    1343             :   {
    1344       36603 :     a = nf_to_scalar_or_basis(nf,a);
    1345       36603 :     if (isint1(a)) a = gen_1;
    1346             :   }
    1347       58254 :   aA = (a == gen_1)? idealhnf_shallow(nf,A): idealmul(nf,a,A);
    1348       58254 :   bB = idealmul(nf,b,B);
    1349       58254 :   d = idealadd(nf,aA,bB);
    1350       58254 :   if (gequal(aA, d)) return zero_nfbezout(nf,aA, a,B,A,pv,pu,pw,pdi);
    1351       24927 :   if (gequal(bB, d)) return zero_nfbezout(nf,bB, b,A,B,pu,pv,pw,pdi);
    1352             :   /* general case is slow */
    1353        4025 :   di = idealHNF_inv(nf,d);
    1354        4025 :   aA = idealmul(nf,aA,di); /* integral */
    1355        4025 :   bB = idealmul(nf,bB,di); /* integral */
    1356             : 
    1357        4025 :   u = idealaddtoone_i(nf, aA, bB);
    1358        4025 :   w = idealmul(nf,aA,B);
    1359        4025 :   v = nfdiv(nf, nfsub(nf, gen_1, u), b);
    1360        4025 :   if (a != gen_1)
    1361             :   {
    1362        1603 :     GEN inva = nfinv(nf, a);
    1363        1603 :     u =  nfmul(nf,u,inva);
    1364        1603 :     w = idealmul(nf, inva, w); /* AB/d */
    1365             :   }
    1366        4025 :   *pu = u;
    1367        4025 :   *pv = v;
    1368        4025 :   *pw = w;
    1369        4025 :   *pdi = di; return d;
    1370             : }
    1371             : /* v a vector of ideals, simplify in place the ones generated by elts of Q */
    1372             : static void
    1373        8274 : idV_simplify(GEN v)
    1374             : {
    1375        8274 :   long i, l = lg(v);
    1376       46655 :   for (i = 1; i < l; i++)
    1377             :   {
    1378       38381 :     GEN M = gel(v,i);
    1379       38381 :     if (typ(M)==t_MAT && RgM_isscalar(M,NULL))
    1380       11858 :       gel(v,i) = Q_abs_shallow(gcoeff(M,1,1));
    1381             :   }
    1382        8274 : }
    1383             : /* Given a torsion-free module x outputs a pseudo-basis for x in HNF */
    1384             : GEN
    1385        6496 : nfhnf0(GEN nf, GEN x, long flag)
    1386             : {
    1387             :   long i, j, def, idef, m, n;
    1388        6496 :   pari_sp av0 = avma, av;
    1389             :   GEN y, A, I, J, U;
    1390             : 
    1391        6496 :   nf = checknf(nf);
    1392        6496 :   check_ZKmodule(x, "nfhnf");
    1393        6496 :   A = gel(x,1); RgM_dimensions(A, &m, &n);
    1394        6496 :   I = gel(x,2);
    1395        6496 :   if (!n) {
    1396          49 :     if (!flag) return gcopy(x);
    1397           0 :     retmkvec2(gcopy(x), cgetg(1,t_MAT));
    1398             :   }
    1399        6447 :   U = flag? matid(n): NULL;
    1400        6447 :   idef = (n < m)? m-n : 0;
    1401        6447 :   av = avma;
    1402        6447 :   A = RgM_to_nfM(nf,A);
    1403        6447 :   I = leafcopy(I);
    1404        6447 :   J = zerovec(n); def = n;
    1405       36015 :   for (i=m; i>idef; i--)
    1406             :   {
    1407       29568 :     GEN d, di = NULL;
    1408             : 
    1409       29568 :     j=def; while (j>=1 && isintzero(gcoeff(A,i,j))) j--;
    1410       29568 :     if (!j)
    1411             :     { /* no pivot on line i */
    1412           7 :       if (idef) idef--;
    1413           7 :       continue;
    1414             :     }
    1415       29561 :     if (j==def) j--;
    1416             :     else {
    1417         490 :       swap(gel(A,j), gel(A,def));
    1418         490 :       swap(gel(I,j), gel(I,def));
    1419         490 :       if (U) swap(gel(U,j), gel(U,def));
    1420             :     }
    1421      175861 :     for (  ; j; j--)
    1422             :     {
    1423      146300 :       GEN a,b, u,v,w, S, T, S0, T0 = gel(A,j);
    1424      146300 :       b = gel(T0,i); if (isintzero(b)) continue;
    1425             : 
    1426       40978 :       S0 = gel(A,def); a = gel(S0,i);
    1427       40978 :       d = nfbezout(nf, a,b, gel(I,def),gel(I,j), &u,&v,&w,&di);
    1428       40978 :       S = colcomb(nf, u,v, S0,T0);
    1429       40978 :       T = colcomb(nf, a,gneg(b), T0,S0);
    1430       40978 :       gel(A,def) = S; gel(A,j) = T;
    1431       40978 :       gel(I,def) = d; gel(I,j) = w;
    1432       40978 :       if (U)
    1433             :       {
    1434          42 :         S0 = gel(U,def);
    1435          42 :         T0 = gel(U,j);
    1436          42 :         gel(U,def) = colcomb(nf, u,v, S0,T0);
    1437          42 :         gel(U,j) = colcomb(nf, a,gneg(b), T0,S0);
    1438             :       }
    1439             :     }
    1440       29561 :     y = gcoeff(A,i,def);
    1441       29561 :     if (!isint1(y))
    1442             :     {
    1443         651 :       GEN yi = nfinv(nf,y);
    1444         651 :       gel(A,def) = nfC_nf_mul(nf, gel(A,def), yi);
    1445         651 :       gel(I,def) = idealmul(nf, y, gel(I,def));
    1446         651 :       if (U) gel(U,def) = nfC_nf_mul(nf, gel(U,def), yi);
    1447         651 :       di = NULL;
    1448             :     }
    1449       29561 :     if (!di) di = idealinv(nf,gel(I,def));
    1450       29561 :     d = gel(I,def);
    1451       29561 :     gel(J,def) = di;
    1452       94647 :     for (j=def+1; j<=n; j++)
    1453             :     {
    1454       65086 :       GEN mc, c = gcoeff(A,i,j); if (isintzero(c)) continue;
    1455       23226 :       c = element_close(nf, c, idealmul(nf,d,gel(J,j)));
    1456       23226 :       mc = gneg(c);
    1457       23226 :       gel(A,j) = colcomb1(nf, mc, gel(A,j),gel(A,def));
    1458       23226 :       if (U) gel(U,j) = colcomb1(nf, mc, gel(U,j),gel(U,def));
    1459             :     }
    1460       29561 :     def--;
    1461       29561 :     if (gc_needed(av,2))
    1462             :     {
    1463           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"nfhnf, i = %ld", i);
    1464           0 :       gerepileall(av,U?4:3, &A,&I,&J,&U);
    1465             :     }
    1466             :   }
    1467        6447 :   n -= def;
    1468        6447 :   A += def; A[0] = evaltyp(t_MAT)|evallg(n+1);
    1469        6447 :   I += def; I[0] = evaltyp(t_VEC)|evallg(n+1);
    1470        6447 :   idV_simplify(I);
    1471        6447 :   x = mkvec2(A,I);
    1472        6447 :   if (U) x = mkvec2(x,U);
    1473        6447 :   return gerepilecopy(av0, x);
    1474             : }
    1475             : 
    1476             : GEN
    1477        6482 : nfhnf(GEN nf, GEN x) { return nfhnf0(nf, x, 0); }
    1478             : 
    1479             : static GEN
    1480           0 : RgV_find_denom(GEN x)
    1481             : {
    1482           0 :   long i, l = lg(x);
    1483           0 :   for (i = 1; i < l; i++)
    1484           0 :     if (Q_denom(gel(x,i)) != gen_1) return gel(x,i);
    1485           0 :   return NULL;
    1486             : }
    1487             : /* A torsion module M over Z_K will be given by a row vector [A,I,J] with
    1488             :  * three components. I=[b_1,...,b_n] is a row vector of n fractional ideals
    1489             :  * given in HNF, J=[a_1,...,a_n] is a row vector of n fractional ideals in
    1490             :  * HNF. A is an nxn matrix (same n) such that if A_j is the j-th column of A
    1491             :  * and e_n is the canonical basis of K^n, then
    1492             :  * M=(b_1e_1+...+b_ne_n)/(a_1A_1+...a_nA_n) */
    1493             : 
    1494             : /* x=[A,I,J] a torsion module as above. Output the
    1495             :  * smith normal form as K=[c_1,...,c_n] such that x = Z_K/c_1+...+Z_K/c_n */
    1496             : GEN
    1497          21 : nfsnf0(GEN nf, GEN x, long flag)
    1498             : {
    1499             :   long i, j, k, l, n, m;
    1500             :   pari_sp av;
    1501             :   GEN z,u,v,w,d,dinv,A,I,J, U,V;
    1502             : 
    1503          21 :   nf = checknf(nf);
    1504          21 :   if (typ(x)!=t_VEC || lg(x)!=4) pari_err_TYPE("nfsnf",x);
    1505          21 :   A = gel(x,1);
    1506          21 :   I = gel(x,2);
    1507          21 :   J = gel(x,3);
    1508          21 :   if (typ(A)!=t_MAT) pari_err_TYPE("nfsnf",A);
    1509          21 :   n = lg(A)-1;
    1510          21 :   if (typ(I)!=t_VEC) pari_err_TYPE("nfsnf",I);
    1511          21 :   if (typ(J)!=t_VEC) pari_err_TYPE("nfsnf",J);
    1512          21 :   if (lg(I)!=n+1 || lg(J)!=n+1) pari_err_DIM("nfsnf");
    1513          21 :   RgM_dimensions(A, &m, &n);
    1514          21 :   if (!n || n != m) pari_err_IMPL("nfsnf for empty or non square matrices");
    1515             : 
    1516          21 :   av = avma;
    1517          21 :   if (!flag) U = V = NULL;
    1518             :   else
    1519             :   {
    1520           7 :     U = matid(m);
    1521           7 :     V = matid(n);
    1522             :   }
    1523          21 :   A = RgM_to_nfM(nf, A);
    1524          21 :   I = leafcopy(I);
    1525          21 :   J = leafcopy(J);
    1526          21 :   for (i = 1; i <= n; i++) gel(J,i) = idealinv(nf, gel(J,i));
    1527          21 :   z = zerovec(n);
    1528         126 :   for (i=n; i>=1; i--)
    1529             :   {
    1530             :     GEN Aii, a, b, db;
    1531         105 :     long c = 0;
    1532         238 :     for (j=i-1; j>=1; j--)
    1533             :     {
    1534         133 :       GEN S, T, S0, T0 = gel(A,j);
    1535         133 :       b = gel(T0,i); if (gequal0(b)) continue;
    1536             : 
    1537          49 :       S0 = gel(A,i); a = gel(S0,i);
    1538          49 :       d = nfbezout(nf, a,b, gel(J,i),gel(J,j), &u,&v,&w,&dinv);
    1539          49 :       S = colcomb(nf, u,v, S0,T0);
    1540          49 :       T = colcomb(nf, a,gneg(b), T0,S0);
    1541          49 :       gel(A,i) = S; gel(A,j) = T;
    1542          49 :       gel(J,i) = d; gel(J,j) = w;
    1543          49 :       if (V)
    1544             :       {
    1545          21 :         T0 = gel(V,j);
    1546          21 :         S0 = gel(V,i);
    1547          21 :         gel(V,i) = colcomb(nf, u,v, S0,T0);
    1548          21 :         gel(V,j) = colcomb(nf, a,gneg(b), T0,S0);
    1549             :       }
    1550             :     }
    1551         238 :     for (j=i-1; j>=1; j--)
    1552             :     {
    1553             :       GEN ri, rj;
    1554         133 :       b = gcoeff(A,j,i); if (gequal0(b)) continue;
    1555             : 
    1556          56 :       a = gcoeff(A,i,i);
    1557          56 :       d = nfbezout(nf, a,b, gel(I,i),gel(I,j), &u,&v,&w,&dinv);
    1558          56 :       ri = rowcomb(nf, u,v,       i,j, A, i);
    1559          56 :       rj = rowcomb(nf, a,gneg(b), j,i, A, i);
    1560         210 :       for (k=1; k<=i; k++) {
    1561         154 :         gcoeff(A,j,k) = gel(rj,k);
    1562         154 :         gcoeff(A,i,k) = gel(ri,k);
    1563             :       }
    1564          56 :       if (U)
    1565             :       {
    1566          21 :         ri = rowcomb(nf, u,v,       i,j, U, m);
    1567          21 :         rj = rowcomb(nf, a,gneg(b), j,i, U, m);
    1568          84 :         for (k=1; k<=m; k++) {
    1569          63 :           gcoeff(U,j,k) = gel(rj,k);
    1570          63 :           gcoeff(U,i,k) = gel(ri,k);
    1571             :         }
    1572             :       }
    1573          56 :       gel(I,i) = d; gel(I,j) = w; c = 1;
    1574             :     }
    1575         147 :     if (c) { i++; continue; }
    1576             : 
    1577          63 :     Aii = gcoeff(A,i,i); if (gequal0(Aii)) continue;
    1578          63 :     gel(J,i) = idealmul(nf, gel(J,i), Aii);
    1579          63 :     gcoeff(A,i,i) = gen_1;
    1580          63 :     if (V) gel(V,i) = nfC_nf_mul(nf, gel(V,i), nfinv(nf,Aii));
    1581          63 :     gel(z,i) = idealmul(nf,gel(J,i),gel(I,i));
    1582          63 :     b = Q_remove_denom(gel(z,i), &db);
    1583         126 :     for (k=1; k<i; k++)
    1584         168 :       for (l=1; l<i; l++)
    1585             :       {
    1586         105 :         GEN d, D, p1, p2, p3, Akl = gcoeff(A,k,l);
    1587             :         long t;
    1588         105 :         if (gequal0(Akl)) continue;
    1589             : 
    1590          91 :         p1 = idealmul(nf,Akl,gel(J,l));
    1591          91 :         p3 = idealmul(nf, p1, gel(I,k));
    1592          91 :         if (db) p3 = RgM_Rg_mul(p3, db);
    1593          91 :         if (RgM_is_ZM(p3) && hnfdivide(b, p3)) continue;
    1594             : 
    1595             :         /* find d in D = I[k]/I[i] not in J[i]/(A[k,l] J[l]) */
    1596           0 :         D = idealdiv(nf,gel(I,k),gel(I,i));
    1597           0 :         p2 = idealdiv(nf,gel(J,i), p1);
    1598           0 :         d = RgV_find_denom( RgM_solve(p2, D) );
    1599           0 :         if (!d) pari_err_BUG("nfsnf");
    1600           0 :         p1 = element_mulvecrow(nf,d,A,k,i);
    1601           0 :         for (t=1; t<=i; t++) gcoeff(A,i,t) = gadd(gcoeff(A,i,t),gel(p1,t));
    1602           0 :         if (U)
    1603             :         {
    1604           0 :           p1 = element_mulvecrow(nf,d,U,k,i);
    1605           0 :           for (t=1; t<=i; t++) gcoeff(U,i,t) = gadd(gcoeff(U,i,t),gel(p1,t));
    1606             :         }
    1607             : 
    1608           0 :         k = i; c = 1; break;
    1609             :       }
    1610          63 :     if (gc_needed(av,1))
    1611             :     {
    1612           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"nfsnf");
    1613           0 :       gerepileall(av,U?6:4, &A,&I,&J,&z,&U,&V);
    1614             :     }
    1615          63 :     if (c) i++; /* iterate on row/column i */
    1616             :   }
    1617          21 :   if (U) z = mkvec3(z,U,V);
    1618          21 :   return gerepilecopy(av, z);
    1619             : }
    1620             : GEN
    1621           0 : nfsnf(GEN nf, GEN x) { return nfsnf0(nf,x,0); }
    1622             : 
    1623             : /* Given a pseudo-basis x, outputs a multiple of its ideal determinant */
    1624             : GEN
    1625          14 : nfdetint(GEN nf, GEN x)
    1626             : {
    1627             :   GEN pass,c,v,det1,piv,pivprec,vi,p1,A,I,id,idprod;
    1628          14 :   long i, j, k, rg, n, m, m1, cm=0, N;
    1629          14 :   pari_sp av = avma, av1;
    1630             : 
    1631          14 :   nf = checknf(nf); N = nf_get_degree(nf);
    1632          14 :   check_ZKmodule(x, "nfdetint");
    1633          14 :   A = gel(x,1);
    1634          14 :   I = gel(x,2);
    1635          14 :   n = lg(A)-1; if (!n) return gen_1;
    1636             : 
    1637          14 :   m1 = lgcols(A); m = m1-1;
    1638          14 :   id = matid(N);
    1639          14 :   c = new_chunk(m1); for (k=1; k<=m; k++) c[k] = 0;
    1640          14 :   piv = pivprec = gen_1;
    1641             : 
    1642          14 :   av1 = avma;
    1643          14 :   det1 = idprod = gen_0; /* dummy for gerepileall */
    1644          14 :   pass = cgetg(m1,t_MAT);
    1645          14 :   v = cgetg(m1,t_COL);
    1646          49 :   for (j=1; j<=m; j++)
    1647             :   {
    1648          35 :     gel(pass,j) = zerocol(m);
    1649          35 :     gel(v,j) = gen_0; /* dummy */
    1650             :   }
    1651          63 :   for (rg=0,k=1; k<=n; k++)
    1652             :   {
    1653          49 :     long t = 0;
    1654         182 :     for (i=1; i<=m; i++)
    1655         133 :       if (!c[i])
    1656             :       {
    1657          77 :         vi=nfmul(nf,piv,gcoeff(A,i,k));
    1658         287 :         for (j=1; j<=m; j++)
    1659         210 :           if (c[j]) vi=gadd(vi,nfmul(nf,gcoeff(pass,i,j),gcoeff(A,j,k)));
    1660          77 :         gel(v,i) = vi; if (!t && !gequal0(vi)) t=i;
    1661             :       }
    1662          49 :     if (t)
    1663             :     {
    1664          49 :       pivprec = piv;
    1665          49 :       if (rg == m-1)
    1666             :       {
    1667          28 :         if (!cm)
    1668             :         {
    1669          14 :           cm=1; idprod = id;
    1670          49 :           for (i=1; i<=m; i++)
    1671          35 :             if (i!=t)
    1672          56 :               idprod = (idprod==id)? gel(I,c[i])
    1673          35 :                                    : idealmul(nf,idprod,gel(I,c[i]));
    1674             :         }
    1675          28 :         p1 = idealmul(nf,gel(v,t),gel(I,k)); c[t]=0;
    1676          28 :         det1 = (typ(det1)==t_INT)? p1: idealadd(nf,p1,det1);
    1677             :       }
    1678             :       else
    1679             :       {
    1680          21 :         rg++; piv=gel(v,t); c[t]=k;
    1681          77 :         for (i=1; i<=m; i++)
    1682          56 :           if (!c[i])
    1683             :           {
    1684         105 :             for (j=1; j<=m; j++)
    1685          77 :               if (c[j] && j!=t)
    1686             :               {
    1687          14 :                 p1 = gsub(nfmul(nf,piv,gcoeff(pass,i,j)),
    1688          14 :                           nfmul(nf,gel(v,i),gcoeff(pass,t,j)));
    1689          21 :                 gcoeff(pass,i,j) = rg>1? nfdiv(nf,p1,pivprec)
    1690          14 :                                        : p1;
    1691             :               }
    1692          28 :             gcoeff(pass,i,t) = gneg(gel(v,i));
    1693             :           }
    1694             :       }
    1695             :     }
    1696          49 :     if (gc_needed(av1,1))
    1697             :     {
    1698           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"nfdetint");
    1699           0 :       gerepileall(av1,6, &det1,&piv,&pivprec,&pass,&v,&idprod);
    1700             :     }
    1701             :   }
    1702          14 :   if (!cm) { avma = av; return cgetg(1,t_MAT); }
    1703          14 :   return gerepileupto(av, idealmul(nf,idprod,det1));
    1704             : }
    1705             : 
    1706             : /* reduce in place components of x[1..lim] mod D (destroy x). D in HNF */
    1707             : static void
    1708       15624 : nfcleanmod(GEN nf, GEN x, long lim, GEN D)
    1709             : {
    1710             :   long i;
    1711             :   GEN DZ, DZ2, dD;
    1712       15624 :   D = Q_remove_denom(D, &dD);
    1713       15624 :   if (dD) x = RgC_Rg_mul(x, dD);
    1714       15624 :   DZ = gcoeff(D,1,1);
    1715       15624 :   DZ2 = shifti(DZ,-1);
    1716       75180 :   for (i=1; i<=lim; i++) {
    1717       59556 :     GEN c = gel(x,i);
    1718       59556 :     c = nf_to_scalar_or_basis(nf, c);
    1719       59556 :     switch(typ(c)) /* c = centermod(c, D) */
    1720             :     {
    1721             :       case t_INT:
    1722       58436 :         if (!signe(c)) break;
    1723       32809 :         c = centermodii(c, DZ, DZ2);
    1724       32809 :         if (dD) c = Qdivii(c,dD);
    1725       32809 :         break;
    1726             :       case t_FRAC: {
    1727          21 :         GEN dc = gel(c,2), nc = gel(c,1), N = mulii(DZ, dc);
    1728          21 :         c = centermodii(nc, N, shifti(N,-1));
    1729          21 :         c = Qdivii(c, dD ? mulii(dc,dD): dc);
    1730          21 :         break;
    1731             :       }
    1732             :       case t_COL: {
    1733             :         GEN dc;
    1734        1099 :         c = Q_remove_denom(c, &dc);
    1735        1099 :         c = ZC_hnfrem(c, dc? ZM_Z_mul(D,dc): D);
    1736        1099 :         if (ZV_isscalar(c))
    1737             :         {
    1738          63 :           c = gel(c,1);
    1739          63 :           if (dD) c = Qdivii(c,dD);
    1740             :         }
    1741             :         else
    1742        1036 :           if (dD) c = RgC_Rg_div(c, dD);
    1743        1099 :         break;
    1744             :       }
    1745             :     }
    1746       59556 :     gel(x,i) = c;
    1747             :   }
    1748       15624 : }
    1749             : 
    1750             : GEN
    1751        1827 : nfhnfmod(GEN nf, GEN x, GEN detmat)
    1752             : {
    1753             :   long li, co, i, j, def, ldef;
    1754        1827 :   pari_sp av0=avma, av;
    1755             :   GEN dA, dI, d0, w, p1, d, u, v, A, I, J, di;
    1756             : 
    1757        1827 :   nf = checknf(nf);
    1758        1827 :   check_ZKmodule(x, "nfhnfmod");
    1759        1827 :   A = gel(x,1);
    1760        1827 :   I = gel(x,2);
    1761        1827 :   co = lg(A); if (co==1) return cgetg(1,t_MAT);
    1762             : 
    1763        1827 :   li = lgcols(A);
    1764        1827 :   if (typ(detmat)!=t_MAT) detmat = idealhnf_shallow(nf, detmat);
    1765        1827 :   detmat = Q_remove_denom(detmat, NULL);
    1766        1827 :   RgM_check_ZM(detmat, "nfhnfmod");
    1767             : 
    1768        1827 :   av = avma;
    1769        1827 :   A = RgM_to_nfM(nf, A);
    1770        1827 :   A = Q_remove_denom(A, &dA);
    1771        1827 :   I = Q_remove_denom(leafcopy(I), &dI);
    1772        1827 :   dA = mul_denom(dA,dI);
    1773        1827 :   if (dA) detmat = ZM_Z_mul(detmat, powiu(dA, minss(li,co)));
    1774             : 
    1775        1827 :   def = co; ldef = (li>co)? li-co+1: 1;
    1776       10647 :   for (i=li-1; i>=ldef; i--)
    1777             :   {
    1778        8820 :     def--; j=def; while (j>=1 && isintzero(gcoeff(A,i,j))) j--;
    1779        8820 :     if (!j) continue;
    1780        8820 :     if (j==def) j--;
    1781             :     else {
    1782        1386 :       swap(gel(A,j), gel(A,def));
    1783        1386 :       swap(gel(I,j), gel(I,def));
    1784             :     }
    1785       40103 :     for (  ; j; j--)
    1786             :     {
    1787       31283 :       GEN a, b, S, T, S0, T0 = gel(A,j);
    1788       31283 :       b = gel(T0,i); if (isintzero(b)) continue;
    1789             : 
    1790        8351 :       S0 = gel(A,def); a = gel(S0,i);
    1791        8351 :       d = nfbezout(nf, a,b, gel(I,def),gel(I,j), &u,&v,&w,&di);
    1792        8351 :       S = colcomb(nf, u,v, S0,T0);
    1793        8351 :       T = colcomb(nf, a,gneg(b), T0,S0);
    1794        8351 :       if (u != gen_0 && v != gen_0) /* already reduced otherwise */
    1795         280 :         nfcleanmod(nf, S, i, idealmul(nf,detmat,di));
    1796        8351 :       nfcleanmod(nf, T, i, idealdiv(nf,detmat,w));
    1797        8351 :       gel(A,def) = S; gel(A,j) = T;
    1798        8351 :       gel(I,def) = d; gel(I,j) = w;
    1799             :     }
    1800        8820 :     if (gc_needed(av,2))
    1801             :     {
    1802           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"[1]: nfhnfmod, i = %ld", i);
    1803           0 :       gerepileall(av,dA? 4: 3, &A,&I,&detmat,&dA);
    1804             :     }
    1805             :   }
    1806        1827 :   def--; d0 = detmat;
    1807        1827 :   A += def; A[0] = evaltyp(t_MAT)|evallg(li);
    1808        1827 :   I += def; I[0] = evaltyp(t_VEC)|evallg(li);
    1809        1827 :   J = cgetg(li,t_VEC);
    1810       10647 :   for (i=li-1; i>=1; i--)
    1811             :   {
    1812        8820 :     GEN b = gcoeff(A,i,i);
    1813        8820 :     d = nfbezout(nf, gen_1,b, d0,gel(I,i), &u,&v,&w,&di);
    1814        8820 :     p1 = nfC_nf_mul(nf,gel(A,i),v);
    1815        8820 :     if (i > 1)
    1816             :     {
    1817        6993 :       d0 = idealmul(nf,d0,di);
    1818        6993 :       nfcleanmod(nf, p1, i, d0);
    1819             :     }
    1820        8820 :     gel(A,i) = p1; gel(p1,i) = gen_1;
    1821        8820 :     gel(I,i) = d;
    1822        8820 :     gel(J,i) = di;
    1823             :   }
    1824        8820 :   for (i=li-2; i>=1; i--)
    1825             :   {
    1826        6993 :     d = gel(I,i);
    1827       28161 :     for (j=i+1; j<li; j++)
    1828             :     {
    1829       21168 :       GEN c = gcoeff(A,i,j); if (isintzero(c)) continue;
    1830        7294 :       c = element_close(nf, c, idealmul(nf,d,gel(J,j)));
    1831        7294 :       gel(A,j) = colcomb1(nf, gneg(c), gel(A,j),gel(A,i));
    1832             :     }
    1833        6993 :     if (gc_needed(av,2))
    1834             :     {
    1835           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"[2]: nfhnfmod, i = %ld", i);
    1836           0 :       gerepileall(av,dA? 4: 3, &A,&I,&J,&dA);
    1837             :     }
    1838             :   }
    1839        1827 :   idV_simplify(I);
    1840        1827 :   if (dA) I = gdiv(I,dA);
    1841        1827 :   return gerepilecopy(av0, mkvec2(A, I));
    1842             : }

Generated by: LCOV version 1.11