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 to exceed 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 - hnf_snf.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.12.1 lcov report (development 24988-2584e74448) Lines: 1524 1669 91.3 %
Date: 2020-01-26 05:57:03 Functions: 88 95 92.6 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2000  The PARI group.
       2             : 
       3             : This file is part of the PARI/GP package.
       4             : 
       5             : PARI/GP is free software; you can redistribute it and/or modify it under the
       6             : terms of the GNU General Public License as published by the Free Software
       7             : Foundation. It is distributed in the hope that it will be useful, but WITHOUT
       8             : ANY WARRANTY WHATSOEVER.
       9             : 
      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             : #include "pari.h"
      15             : #include "paripriv.h"
      16             : /**************************************************************/
      17             : /**                                                          **/
      18             : /**                HERMITE NORMAL FORM REDUCTION             **/
      19             : /**                                                          **/
      20             : /**************************************************************/
      21             : static GEN ZV_hnfgcdext(GEN A);
      22             : static GEN
      23          14 : hnfallgen(GEN x)
      24             : {
      25          14 :   GEN z = cgetg(3, t_VEC);
      26          14 :   gel(z,1) = RgM_hnfall(x, (GEN*)(z+2), 1);
      27          14 :   return z;
      28             : }
      29             : GEN
      30         252 : mathnf0(GEN x, long flag)
      31             : {
      32         252 :   switch(typ(x))
      33             :   {
      34             :     case t_VEC:
      35          70 :       if (RgV_is_ZV(x))
      36          42 :         switch (flag)
      37             :         {
      38             :           case 0:
      39          14 :             if (lg(x) == 1) return cgetg(1, t_MAT);
      40           7 :             retmkmat(mkcol(ZV_content(x)));
      41             :           case 1:
      42             :           case 4:
      43          21 :             return ZV_hnfgcdext(x);
      44             :         }
      45          35 :       x = gtomat(x); break;
      46         182 :     case t_MAT: break;
      47           0 :     default: pari_err_TYPE("mathnf0",x);
      48             :   }
      49             : 
      50         217 :   switch(flag)
      51             :   {
      52         168 :     case 0: case 2: return RgM_is_ZM(x)? ZM_hnf(x): RgM_hnfall(x,NULL,1);
      53          28 :     case 1: case 3: return RgM_is_ZM(x)? hnfall(x): hnfallgen(x);
      54           7 :     case 4: RgM_check_ZM(x, "mathnf0"); return hnflll(x);
      55          14 :     case 5: RgM_check_ZM(x, "mathnf0"); return hnfperm(x);
      56           0 :     default: pari_err_FLAG("mathnf");
      57             :   }
      58             :   return NULL; /* LCOV_EXCL_LINE */
      59             : }
      60             : 
      61             : /*******************************************************************/
      62             : /*                                                                 */
      63             : /*                SPECIAL HNF (FOR INTERNAL USE !!!)               */
      64             : /*                                                                 */
      65             : /*******************************************************************/
      66             : static int
      67     3777716 : count(GEN mat, long row, long len, long *firstnonzero)
      68             : {
      69     3777716 :   long j, n = 0;
      70             : 
      71   320556270 :   for (j=1; j<=len; j++)
      72             :   {
      73   317922735 :     long p = mael(mat,j,row);
      74   317922735 :     if (p)
      75             :     {
      76    10334048 :       if (labs(p)!=1) return -1;
      77     9189867 :       n++; *firstnonzero=j;
      78             :     }
      79             :   }
      80     2633535 :   return n;
      81             : }
      82             : 
      83             : static int
      84      157736 : count2(GEN mat, long row, long len)
      85             : {
      86             :   long j;
      87     4491563 :   for (j=len; j; j--)
      88     4413239 :     if (labs(mael(mat,j,row)) == 1) return j;
      89       78324 :   return 0;
      90             : }
      91             : 
      92             : static GEN
      93       49671 : hnffinal(GEN matgen,GEN perm,GEN* ptdep,GEN* ptB,GEN* ptC)
      94             : {
      95             :   GEN p1,p2,U,H,Hnew,Bnew,Cnew,diagH1;
      96       49671 :   GEN B = *ptB, C = *ptC, dep = *ptdep, depnew;
      97             :   pari_sp av;
      98             :   long i,j,k,s,i1,j1,zc;
      99       49671 :   long co = lg(C);
     100       49671 :   long col = lg(matgen)-1;
     101             :   long lnz, nlze, lig;
     102             : 
     103       49671 :   if (col == 0) return matgen;
     104       49160 :   lnz = nbrows(matgen); /* was called lnz-1 - nr in hnfspec */
     105       49160 :   nlze = nbrows(dep);
     106       49160 :   lig = nlze + lnz;
     107             :   /* H: lnz x lnz [disregarding initial 0 cols], U: col x col */
     108       49160 :   H = ZM_hnflll(matgen, &U, 0);
     109       49160 :   H += (lg(H)-1 - lnz); H[0] = evaltyp(t_MAT) | evallg(lnz+1);
     110             :   /* Only keep the part above the H (above the 0s is 0 since the dep rows
     111             :    * are dependent from the ones in matgen) */
     112       49160 :   zc = col - lnz; /* # of 0 columns, correspond to units */
     113       49160 :   if (nlze) { dep = ZM_mul(dep,U); dep += zc; }
     114             : 
     115       49160 :   diagH1 = new_chunk(lnz+1); /* diagH1[i] = 0 iff H[i,i] != 1 (set later) */
     116             : 
     117       49160 :   av = avma;
     118       49160 :   Cnew = cgetg(co, typ(C));
     119       49160 :   setlg(C, col+1); p1 = gmul(C,U);
     120       49160 :   for (j=1; j<=col; j++) gel(Cnew,j) = gel(p1,j);
     121       49160 :   for (   ; j<co ; j++)  gel(Cnew,j) = gel(C,j);
     122             : 
     123             :   /* Clean up B using new H */
     124      233199 :   for (s=0,i=lnz; i; i--)
     125             :   {
     126      184039 :     GEN Di = gel(dep,i), Hi = gel(H,i);
     127      184039 :     GEN h = gel(Hi,i); /* H[i,i] */
     128      184039 :     if ( (diagH1[i] = is_pm1(h)) ) { h = NULL; s++; }
     129    11506009 :     for (j=col+1; j<co; j++)
     130             :     {
     131    11321970 :       GEN z = gel(B,j-col);
     132    11321970 :       p1 = gel(z,i+nlze);
     133    11321970 :       if (h) p1 = truedivii(p1,h);
     134    11321970 :       if (!signe(p1)) continue;
     135     6805008 :       for (k=1; k<=nlze; k++) gel(z,k) = subii(gel(z,k), mulii(p1, gel(Di,k)));
     136     6805008 :       for (   ; k<=lig;  k++) gel(z,k) = subii(gel(z,k), mulii(p1, gel(Hi,k-nlze)));
     137     6805008 :       gel(Cnew,j) = gsub(gel(Cnew,j), gmul(p1, gel(Cnew,i+zc)));
     138             :     }
     139      184039 :     if (gc_needed(av,2))
     140             :     {
     141         375 :       if(DEBUGMEM>1) pari_warn(warnmem,"hnffinal, i = %ld",i);
     142         375 :       gerepileall(av, 2, &Cnew, &B);
     143             :     }
     144             :   }
     145       49160 :   p1 = cgetg(lnz+1,t_VEC); p2 = perm + nlze;
     146      233199 :   for (i1=0, j1=lnz-s, i=1; i<=lnz; i++) /* push the 1 rows down */
     147      184039 :     if (diagH1[i])
     148       85703 :       gel(p1,++j1) = gel(p2,i);
     149             :     else
     150       98336 :       gel(p2,++i1) = gel(p2,i);
     151       49160 :   for (i=i1+1; i<=lnz; i++) gel(p2,i) = gel(p1,i);
     152             : 
     153             :   /* s = # extra redundant generators taken from H
     154             :    *          zc  col-s  co   zc = col - lnz
     155             :    *       [ 0 |dep |     ]    i = nlze + lnz - s = lig - s
     156             :    *  nlze [--------|  B' ]
     157             :    *       [ 0 | H' |     ]    H' = H minus the s rows with a 1 on diagonal
     158             :    *     i [--------|-----] lig-s           (= "1-rows")
     159             :    *       [   0    | Id  ]
     160             :    *       [        |     ] li */
     161       49160 :   lig -= s; col -= s; lnz -= s;
     162       49160 :   Hnew = cgetg(lnz+1,t_MAT);
     163       49160 :   depnew = cgetg(lnz+1,t_MAT); /* only used if nlze > 0 */
     164       49160 :   Bnew = cgetg(co-col,t_MAT);
     165       49160 :   C = shallowcopy(Cnew);
     166      233199 :   for (j=1,i1=j1=0; j<=lnz+s; j++)
     167             :   {
     168      184039 :     GEN z = gel(H,j);
     169      184039 :     if (diagH1[j])
     170             :     { /* hit exactly s times */
     171       85703 :       i1++; C[i1+col] = Cnew[j+zc];
     172       85703 :       p1 = cgetg(lig+1,t_COL); gel(Bnew,i1) = p1;
     173       85703 :       for (i=1; i<=nlze; i++) gel(p1,i) = gcoeff(dep,i,j);
     174       85703 :       p1 += nlze;
     175             :     }
     176             :     else
     177             :     {
     178       98336 :       j1++; C[j1+zc] = Cnew[j+zc];
     179       98336 :       p1 = cgetg(lnz+1,t_COL); gel(Hnew,j1) = p1;
     180       98336 :       depnew[j1] = dep[j];
     181             :     }
     182      953524 :     for (i=k=1; k<=lnz; i++)
     183      769485 :       if (!diagH1[i]) p1[k++] = z[i];
     184             :   }
     185     1890937 :   for (j=s+1; j<co-col; j++)
     186             :   {
     187     1841777 :     GEN z = gel(B,j-s);
     188     1841777 :     p1 = cgetg(lig+1,t_COL); gel(Bnew,j) = p1;
     189     1841777 :     for (i=1; i<=nlze; i++) gel(p1,i) = gel(z,i);
     190     1841777 :     z += nlze; p1 += nlze;
     191     9355648 :     for (i=k=1; k<=lnz; i++)
     192     7513871 :       if (!diagH1[i]) gel(p1,k++) = gel(z,i);
     193             :   }
     194       49160 :   *ptdep = depnew;
     195       49160 :   *ptC = C;
     196       49160 :   *ptB = Bnew; return Hnew;
     197             : }
     198             : 
     199             : /* for debugging */
     200             : static void
     201           0 : p_mat(GEN mat, GEN perm, long k)
     202             : {
     203           0 :   pari_sp av = avma;
     204           0 :   perm = vecslice(perm, k+1, lg(perm)-1);
     205           0 :   err_printf("Permutation: %Ps\n",perm);
     206           0 :   if (DEBUGLEVEL > 6)
     207           0 :     err_printf("matgen = %Ps\n", zm_to_ZM( rowpermute(mat, perm) ));
     208           0 :   set_avma(av);
     209           0 : }
     210             : 
     211             : static GEN
     212      688240 : col_dup(long l, GEN col)
     213             : {
     214      688240 :   GEN c = new_chunk(l);
     215      688240 :   memcpy(c,col,l * sizeof(long)); return c;
     216             : }
     217             : 
     218             : /* permutation giving imagecompl(x') | image(x'), x' = transpose of x */
     219             : static GEN
     220       49160 : ZM_rowrankprofile(GEN x, long *nlze)
     221             : {
     222       49160 :   pari_sp av = avma;
     223             :   GEN d, y;
     224             :   long i, j, k, l, r;
     225             : 
     226       49160 :   x = shallowtrans(x); l = lg(x);
     227       49160 :   (void)new_chunk(l); /* HACK */
     228       49160 :   d = ZM_pivots(x,&r); set_avma(av);
     229       49160 :   *nlze = r;
     230       49160 :   if (!d) return identity_perm(l-1);
     231       47387 :   y = cgetg(l,t_VECSMALL);
     232      356416 :   for (i = j = 1, k = r+1; i<l; i++)
     233      309029 :     if (d[i]) y[k++] = i; else y[j++] = i;
     234       47387 :   return y;
     235             : }
     236             : 
     237             : /* HNF reduce a relation matrix (column operations + row permutation)
     238             : ** Input:
     239             : **   mat = (li-1) x (co-1) matrix of long
     240             : **   C   = r x (co-1) matrix of GEN
     241             : **   perm= permutation vector (length li-1), indexing the rows of mat: easier
     242             : **     to maintain perm than to copy rows. For columns we can do it directly
     243             : **     using e.g. swap(mat[i], mat[j])
     244             : **   k0 = integer. The k0 first lines of mat are dense, the others are sparse.
     245             : ** Output: cf ASCII art in the function body
     246             : **
     247             : ** row permutations applied to perm
     248             : ** column operations applied to C. IN PLACE
     249             : **/
     250             : GEN
     251       28554 : hnfspec_i(GEN mat0, GEN perm, GEN* ptdep, GEN* ptB, GEN* ptC, long k0)
     252             : {
     253             :   pari_sp av;
     254             :   long co, n, s, nlze, lnz, nr, i, j, k, lk0, col, lig, *p;
     255             :   GEN mat;
     256             :   GEN p1, p2, matb, matbnew, vmax, matt, T, extramat, B, C, H, dep, permpro;
     257       28554 :   const long li = lg(perm); /* = lgcols(mat0) */
     258       28554 :   const long CO = lg(mat0);
     259             : 
     260       28554 :   n = 0; /* -Wall */
     261             : 
     262       28554 :   C = *ptC; co = CO;
     263       28554 :   if (co > 300 && co > 1.5 * li)
     264             :   { /* treat the rest at the end */
     265           0 :     co = (long)(1.2 * li);
     266           0 :     setlg(C, co);
     267             :   }
     268             : 
     269       28554 :   if (DEBUGLEVEL>5)
     270             :   {
     271           0 :     err_printf("Entering hnfspec\n");
     272           0 :     p_mat(mat0,perm,0);
     273             :   }
     274       28554 :   matt = cgetg(co, t_MAT); /* dense part of mat (top) */
     275       28554 :   mat = cgetg(co, t_MAT);
     276      716794 :   for (j = 1; j < co; j++)
     277             :   {
     278      688240 :     GEN matj = col_dup(li, gel(mat0,j));
     279      688240 :     p1 = cgetg(k0+1,t_COL); gel(matt,j) = p1; gel(mat,j) = matj;
     280      688240 :     for (i=1; i<=k0; i++) gel(p1,i) = stoi(matj[perm[i]]);
     281             :   }
     282       28554 :   av = avma;
     283             : 
     284       28554 :   i = lig = li-1; col = co-1; lk0 = k0;
     285       28554 :   T = (k0 || (lg(C) > 1 && lgcols(C) > 1))? matid(col): NULL;
     286             :   /* Look for lines with a single non-0 entry, equal to 1 in absolute value */
     287     2830801 :   while (i > lk0 && col)
     288     2773693 :     switch( count(mat,perm[i],col,&n) )
     289             :     {
     290             :       case 0: /* move zero lines between k0+1 and lk0 */
     291       21792 :         lk0++; lswap(perm[i], perm[lk0]);
     292       21792 :         i = lig; continue;
     293             : 
     294             :       case 1: /* move trivial generator between lig+1 and li */
     295      139208 :         lswap(perm[i], perm[lig]);
     296      139208 :         if (T) swap(gel(T,n), gel(T,col));
     297      139208 :         swap(gel(mat,n), gel(mat,col)); p = gel(mat,col);
     298      139208 :         if (p[perm[lig]] < 0) /* = -1 */
     299             :         { /* convert relation -g = 0 to g = 0 */
     300      113338 :           for (i=lk0+1; i<lig; i++) p[perm[i]] = -p[perm[i]];
     301      113338 :           if (T)
     302             :           {
     303      113338 :             p1 = gel(T,col);
     304     5854201 :             for (i=1; ; i++) /* T is a permuted identity: single non-0 entry */
     305    11595064 :               if (signe(gel(p1,i))) { togglesign_safe(&gel(p1,i)); break; }
     306             :           }
     307             :         }
     308      139208 :         lig--; col--; i = lig; continue;
     309             : 
     310     2612693 :       default: i--;
     311             :     }
     312       28554 :   if (DEBUGLEVEL>5) { err_printf("    after phase1:\n"); p_mat(mat,perm,0); }
     313             : 
     314             : #define absmax(s,z) {long _z; _z = labs(z); if (_z > s) s = _z;}
     315             :   /* Get rid of all lines containing only 0 and +/- 1, keeping track of column
     316             :    * operations in T. Leave the rows 1..lk0 alone [up to k0, coefficient
     317             :    * explosion, between k0+1 and lk0, row is 0] */
     318       28554 :   s = 0;
     319      284582 :   while (lig > lk0 && col && s < (long)(HIGHBIT>>1))
     320             :   {
     321     1029125 :     for (i=lig; i>lk0; i--)
     322     1004023 :       if (count(mat,perm[i],col,&n) > 0) break;
     323      252576 :     if (i == lk0) break;
     324             : 
     325             :     /* only 0, +/- 1 entries, at least 2 of them non-zero */
     326      227474 :     lswap(perm[i], perm[lig]);
     327      227474 :     swap(gel(mat,n), gel(mat,col)); p = gel(mat,col);
     328      227474 :     if (T) swap(gel(T,n), gel(T,col));
     329      227474 :     if (p[perm[lig]] < 0)
     330             :     {
     331      157168 :       for (i=lk0+1; i<=lig; i++) p[perm[i]] = -p[perm[i]];
     332      157168 :       if (T) ZV_togglesign(gel(T,col));
     333             :     }
     334    10311692 :     for (j=1; j<col; j++)
     335             :     {
     336    10084218 :       GEN matj = gel(mat,j);
     337             :       long t;
     338    10084218 :       if (! (t = matj[perm[lig]]) ) continue;
     339      621532 :       if (t == 1) {
     340      301283 :         for (i=lk0+1; i<=lig; i++) absmax(s, matj[perm[i]] -= p[perm[i]]);
     341             :       }
     342             :       else { /* t = -1 */
     343      320249 :         for (i=lk0+1; i<=lig; i++) absmax(s, matj[perm[i]] += p[perm[i]]);
     344             :       }
     345      621532 :       if (T) ZC_lincomb1_inplace(gel(T,j), gel(T,col), stoi(-t));
     346             :     }
     347      227474 :     lig--; col--;
     348      227474 :     if (gc_needed(av,3))
     349             :     {
     350           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"hnfspec[1]");
     351           0 :       if (T) T = gerepilecopy(av, T); else set_avma(av);
     352             :     }
     353             :   }
     354             :   /* As above with lines containing a +/- 1 (no other assumption).
     355             :    * Stop when single precision becomes dangerous */
     356       28554 :   vmax = cgetg(co,t_VECSMALL);
     357      350112 :   for (j=1; j<=col; j++)
     358             :   {
     359      321558 :     GEN matj = gel(mat,j);
     360      321558 :     for (s=0, i=lk0+1; i<=lig; i++) absmax(s, matj[i]);
     361      321558 :     vmax[j] = s;
     362             :   }
     363      136515 :   while (lig > lk0 && col)
     364             :   {
     365      164706 :     for (i=lig; i>lk0; i--)
     366      157736 :       if ( (n = count2(mat,perm[i],col)) ) break;
     367       86382 :     if (i == lk0) break;
     368             : 
     369       79412 :     lswap(vmax[n], vmax[col]);
     370       79412 :     lswap(perm[i], perm[lig]);
     371       79412 :     swap(gel(mat,n), gel(mat,col)); p = gel(mat,col);
     372       79412 :     if (T) swap(gel(T,n), gel(T,col));
     373       79412 :     if (p[perm[lig]] < 0)
     374             :     {
     375       33432 :       for (i=lk0+1; i<=lig; i++) p[perm[i]] = -p[perm[i]];
     376       33432 :       if (T) ZV_togglesign(gel(T,col));
     377             :     }
     378     1405246 :     for (j=1; j<col; j++)
     379             :     {
     380     1325839 :       GEN matj = gel(mat,j);
     381             :       long t;
     382     1325839 :       if (! (t = matj[perm[lig]]) ) continue;
     383      714467 :       if (vmax[col] && (ulong)labs(t) >= (HIGHBIT-vmax[j]) / vmax[col])
     384           5 :         goto END2;
     385             : 
     386      714462 :       for (s=0, i=lk0+1; i<=lig; i++) absmax(s, matj[perm[i]] -= t*p[perm[i]]);
     387      714462 :       vmax[j] = s;
     388      714462 :       if (T) ZC_lincomb1_inplace(gel(T,j), gel(T,col), stoi(-t));
     389             :     }
     390       79407 :     lig--; col--;
     391       79407 :     if (gc_needed(av,3))
     392             :     {
     393           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"hnfspec[2]");
     394           0 :       gerepileall(av, T? 2: 1, &vmax, &T);
     395             :     }
     396             :   }
     397             : 
     398             : END2: /* clean up mat: remove everything to the right of the 1s on diagonal */
     399             :   /* go multiprecision first */
     400       28554 :   matb = cgetg(co,t_MAT); /* bottom part (complement of matt) */
     401      716794 :   for (j=1; j<co; j++)
     402             :   {
     403      688240 :     GEN matj = gel(mat,j);
     404      688240 :     p1 = cgetg(li-k0,t_COL); gel(matb,j) = p1;
     405      688240 :     p1 -= k0;
     406      688240 :     for (i=k0+1; i<li; i++) gel(p1,i) = stoi(matj[perm[i]]);
     407             :   }
     408       28554 :   if (DEBUGLEVEL>5)
     409             :   {
     410           0 :     err_printf("    after phase2:\n");
     411           0 :     p_mat(mat,perm,lk0);
     412             :   }
     413      446242 :   for (i=li-2; i>lig; i--)
     414             :   {
     415      417688 :     long h, i0 = i - k0, k = i + co-li;
     416      417688 :     GEN Bk = gel(matb,k);
     417    15616097 :     for (j=k+1; j<co; j++)
     418             :     {
     419    15198409 :       GEN Bj = gel(matb,j), v = gel(Bj,i0);
     420    15198409 :       s = signe(v); if (!s) continue;
     421             : 
     422     3112634 :       gel(Bj,i0) = gen_0;
     423     3112634 :       if (is_pm1(v))
     424             :       {
     425     1566331 :         if (s > 0) /* v = 1 */
     426      781279 :         { for (h=1; h<i0; h++) gel(Bj,h) = subii(gel(Bj,h), gel(Bk,h)); }
     427             :         else /* v = -1 */
     428      785052 :         { for (h=1; h<i0; h++) gel(Bj,h) = addii(gel(Bj,h), gel(Bk,h)); }
     429             :       }
     430             :       else {
     431     1546303 :         for (h=1; h<i0; h++) gel(Bj,h) = subii(gel(Bj,h), mulii(v,gel(Bk,h)));
     432             :       }
     433     3112634 :       if (T) ZC_lincomb1_inplace(gel(T,j), gel(T,k), negi(v));
     434             :     }
     435      417688 :     if (gc_needed(av,2))
     436             :     {
     437          42 :       if(DEBUGMEM>1) pari_warn(warnmem,"hnfspec[3], i = %ld", i);
     438          42 :       for (h=1; h<co; h++) setlg(matb[h], i0+1); /* bottom can be forgotten */
     439          42 :       gerepileall(av, T? 2: 1, &matb, &T);
     440             :     }
     441             :   }
     442       28554 :   for (j=1; j<co; j++) setlg(matb[j], lig-k0+1); /* bottom can be forgotten */
     443       28554 :   gerepileall(av, T? 2: 1, &matb, &T);
     444       28554 :   if (DEBUGLEVEL>5) err_printf("    matb cleaned up (using Id block)\n");
     445             : 
     446       28554 :   nlze = lk0 - k0;  /* # of 0 rows */
     447       28554 :   lnz = lig-nlze+1; /* 1 + # of non-0 rows (!= 0...0 1 0 ... 0) */
     448       28554 :   if (T) matt = ZM_mul(matt,T); /* update top rows */
     449       28554 :   extramat = cgetg(col+1,t_MAT); /* = new C minus the 0 rows */
     450      270705 :   for (j=1; j<=col; j++)
     451             :   {
     452      242151 :     GEN z = gel(matt,j);
     453      242151 :     GEN t = (gel(matb,j)) + nlze - k0;
     454      242151 :     p2=cgetg(lnz,t_COL); gel(extramat,j) = p2;
     455      242151 :     for (i=1; i<=k0; i++) gel(p2,i) = gel(z,i); /* top k0 rows */
     456      242151 :     for (   ; i<lnz; i++) gel(p2,i) = gel(t,i); /* other non-0 rows */
     457             :   }
     458       28554 :   if (!col) {
     459         511 :     permpro = identity_perm(lnz);
     460         511 :     nr = lnz;
     461             :   }
     462             :   else
     463       28043 :     permpro = ZM_rowrankprofile(extramat, &nr);
     464             :   /* lnz = lg(permpro) */
     465       28554 :   if (nlze)
     466             :   { /* put the nlze 0 rows (trivial generators) at the top */
     467        4012 :     p1 = new_chunk(lk0+1);
     468        4012 :     for (i=1; i<=nlze; i++) p1[i] = perm[i + k0];
     469        4012 :     for (   ; i<=lk0; i++)  p1[i] = perm[i - nlze];
     470        4012 :     for (i=1; i<=lk0; i++)  perm[i] = p1[i];
     471             :   }
     472             :   /* sort other rows according to permpro (nr redundant generators first) */
     473       28554 :   p1 = new_chunk(lnz); p2 = perm + nlze;
     474       28554 :   for (i=1; i<lnz; i++) p1[i] = p2[permpro[i]];
     475       28554 :   for (i=1; i<lnz; i++) p2[i] = p1[i];
     476             :   /* perm indexes the rows of mat
     477             :    *   |_0__|__redund__|__dense__|__too big__|_____done______|
     478             :    *   0  nlze                              lig             li
     479             :    *         \___nr___/ \___k0__/
     480             :    *         \____________lnz ______________/
     481             :    *
     482             :    *               col   co
     483             :    *       [dep     |     ]
     484             :    *    i0 [--------|  B  ] (i0 = nlze + nr)
     485             :    *       [matbnew |     ] matbnew has maximal rank = lnz-1 - nr
     486             :    * mat = [--------|-----] lig
     487             :    *       [   0    | Id  ]
     488             :    *       [        |     ] li */
     489             : 
     490       28554 :   matbnew = cgetg(col+1,t_MAT); /* dense+toobig, maximal rank. For hnffinal */
     491       28554 :   dep    = cgetg(col+1,t_MAT); /* rows dependent from the ones in matbnew */
     492      270705 :   for (j=1; j<=col; j++)
     493             :   {
     494      242151 :     GEN z = gel(extramat,j);
     495      242151 :     p1 = cgetg(nlze+nr+1,t_COL); gel(dep,j) = p1;
     496      242151 :     p2 = cgetg(lnz-nr,t_COL); gel(matbnew,j) = p2;
     497      242151 :     for (i=1; i<=nlze; i++) gel(p1,i) = gen_0;
     498      242151 :     p1 += nlze; for (i=1; i<=nr; i++) p1[i] = z[permpro[i]];
     499      242151 :     p2 -= nr;   for (   ; i<lnz; i++) p2[i] = z[permpro[i]];
     500             :   }
     501             : 
     502             :   /* redundant generators in terms of the genuine generators
     503             :    * (x_i) = - (g_i) B */
     504       28554 :   B = cgetg(co-col,t_MAT);
     505      474643 :   for (j=col+1; j<co; j++)
     506             :   {
     507      446089 :     GEN y = gel(matt,j);
     508      446089 :     GEN z = gel(matb,j);
     509      446089 :     p1=cgetg(lig+1,t_COL); gel(B,j-col) = p1;
     510      446089 :     for (i=1; i<=nlze; i++) gel(p1,i) = gel(z,i);
     511      446089 :     p1 += nlze; z += nlze-k0;
     512     4437329 :     for (k=1; k<lnz; k++)
     513             :     {
     514     3991240 :       i = permpro[k];
     515     3991240 :       gel(p1,k) = (i <= k0)? gel(y,i): gel(z,i);
     516             :     }
     517             :   }
     518       28554 :   if (T) C = typ(C)==t_MAT? RgM_mul(C,T): RgV_RgM_mul(C,T);
     519       28554 :   gerepileall(av, 4, &matbnew, &B, &dep, &C);
     520       28554 :   *ptdep = dep;
     521       28554 :   *ptB = B;
     522       28554 :   H = hnffinal(matbnew, perm, ptdep, ptB, &C);
     523       28554 :   if (CO > co)
     524             :   { /* treat the rest, N cols at a time (hnflll slow otherwise) */
     525           0 :     const long N = 300;
     526           0 :     long a, L = CO - co, l = minss(L, N); /* L columns to add */
     527           0 :     GEN CC = *ptC, m0 = mat0;
     528           0 :     setlg(CC, CO); /* restore */
     529           0 :     CC += co-1;
     530           0 :     m0 += co-1;
     531           0 :     for (a = l;;)
     532           0 :     {
     533             :       GEN MAT, emb;
     534           0 :       gerepileall(av, 4, &H,&C,ptB,ptdep);
     535           0 :       MAT = cgetg(l + 1, t_MAT);
     536           0 :       emb = cgetg(l + 1, typ(C));
     537           0 :       for (j = 1 ; j <= l; j++)
     538             :       {
     539           0 :         gel(MAT,j) = gel(m0,j);
     540           0 :         emb[j] = CC[j];
     541             :       }
     542           0 :       H = hnfadd_i(H, perm, ptdep, ptB, &C, MAT, emb);
     543           0 :       if (a == L) break;
     544           0 :       CC += l;
     545           0 :       m0 += l;
     546           0 :       a += l; if (a > L) { l = L - (a - l); a = L; }
     547             :     }
     548             :   }
     549       28554 :   *ptC = C; return H;
     550             : }
     551             : 
     552             : GEN
     553         343 : hnfspec(GEN mat, GEN perm, GEN* ptdep, GEN* ptB, GEN* ptC, long k0)
     554             : {
     555         343 :   pari_sp av = avma;
     556         343 :   GEN H = hnfspec_i(mat, perm, ptdep, ptB, ptC, k0);
     557         343 :   gerepileall(av, 4, ptC, ptdep, ptB, &H); return H;
     558             : }
     559             : 
     560             : /* HNF reduce x, apply same transforms to C */
     561             : GEN
     562         343 : mathnfspec(GEN x, GEN *ptperm, GEN *ptdep, GEN *ptB, GEN *ptC)
     563             : {
     564         343 :   long i,j,k,ly,lx = lg(x);
     565             :   GEN z, perm;
     566         343 :   if (lx == 1) return cgetg(1, t_MAT);
     567         343 :   ly = lgcols(x);
     568         343 :   *ptperm = perm = identity_perm(ly-1);
     569         343 :   z = cgetg(lx,t_MAT);
     570        2674 :   for (i=1; i<lx; i++)
     571             :   {
     572        2331 :     GEN C = cgetg(ly,t_COL), D = gel(x,i);
     573        2331 :     gel(z,i) = C;
     574       88578 :     for (j=1; j<ly; j++)
     575             :     {
     576       86247 :       GEN d = gel(D,j);
     577       86247 :       if (is_bigint(d)) goto TOOLARGE;
     578       86247 :       C[j] = itos(d);
     579             :     }
     580             :   }
     581             :   /*  [ dep |     ]
     582             :    *  [-----|  B  ]
     583             :    *  [  H  |     ]
     584             :    *  [-----|-----]
     585             :    *  [  0  | Id  ] */
     586         343 :   return hnfspec(z,perm, ptdep, ptB, ptC, 0);
     587             : 
     588             : TOOLARGE:
     589           0 :   if (lg(*ptC) > 1 && lgcols(*ptC) > 1)
     590           0 :     pari_err_IMPL("mathnfspec with large entries");
     591           0 :   x = ZM_hnf(x); lx = lg(x); j = ly; k = 0;
     592           0 :   for (i=1; i<ly; i++)
     593             :   {
     594           0 :     if (equali1(gcoeff(x,i,i + lx-ly)))
     595           0 :       perm[--j] = i;
     596             :     else
     597           0 :       perm[++k] = i;
     598             :   }
     599           0 :   setlg(perm,k+1);
     600           0 :   x = rowpermute(x, perm); /* upper part */
     601           0 :   setlg(perm,ly);
     602           0 :   *ptB = vecslice(x, j+lx-ly, lx-1);
     603           0 :   setlg(x, j);
     604           0 :   *ptdep = rowslice(x, 1, lx-ly);
     605           0 :   return rowslice(x, lx-ly+1, k); /* H */
     606             : }
     607             : 
     608             : /* add new relations to a matrix treated by hnfspec (extramat / extraC) */
     609             : GEN
     610       21117 : hnfadd_i(GEN H, GEN perm, GEN* ptdep, GEN* ptB, GEN* ptC, /* cf hnfspec */
     611             :        GEN extramat,GEN extraC)
     612             : {
     613       21117 :   GEN matb, extratop, Cnew, permpro, B = *ptB, C = *ptC, dep = *ptdep;
     614             :   long i, lH, lB, li, lig, co, col, nlze;
     615             : 
     616       21117 :   if (lg(extramat) == 1) return H;
     617       21117 :   co = lg(C)-1;
     618       21117 :   lH = lg(H)-1;
     619       21117 :   lB = lg(B)-1;
     620       21117 :   li = lg(perm)-1;
     621       21117 :   lig = li - lB;
     622       21117 :   col = co - lB;
     623       21117 :   nlze = lig - lH;
     624             : 
     625             :  /*               col    co
     626             :   *       [ 0 |dep |     ]
     627             :   *  nlze [--------|  B  ]
     628             :   *       [ 0 | H  |     ]
     629             :   *       [--------|-----] lig
     630             :   *       [   0    | Id  ]
     631             :   *       [        |     ] li */
     632       21117 :   extratop = zm_to_ZM( rowslicepermute(extramat, perm, 1, lig) );
     633       21117 :   if (li != lig)
     634             :   { /* zero out bottom part, using the Id block */
     635       21117 :     GEN A = vecslice(C, col+1, co);
     636       21117 :     GEN c = rowslicepermute(extramat, perm, lig+1, li);
     637       21117 :     extraC   = gsub(extraC, typ(A)==t_MAT? RgM_zm_mul(A, c): RgV_zm_mul(A,c));
     638       21117 :     extratop = ZM_sub(extratop, ZM_zm_mul(B, c));
     639             :   }
     640             : 
     641       21117 :   extramat = shallowconcat(extratop, vconcat(dep, H));
     642       21117 :   Cnew     = shallowconcat(extraC, vecslice(C, col-lH+1, co));
     643       21117 :   if (DEBUGLEVEL>5) err_printf("    1st phase done\n");
     644       21117 :   permpro = ZM_rowrankprofile(extramat, &nlze);
     645       21117 :   extramat = rowpermute(extramat, permpro);
     646       21117 :   *ptB     = rowpermute(B,        permpro);
     647       21117 :   permpro = vecsmallpermute(perm, permpro);
     648       21117 :   for (i=1; i<=lig; i++) perm[i] = permpro[i]; /* perm o= permpro */
     649             : 
     650       21117 :   *ptdep  = rowslice(extramat, 1, nlze);
     651       21117 :   matb    = rowslice(extramat, nlze+1, lig);
     652       21117 :   if (DEBUGLEVEL>5) err_printf("    2nd phase done\n");
     653       21117 :   H = hnffinal(matb,perm,ptdep,ptB,&Cnew);
     654       21117 :   *ptC = shallowconcat(vecslice(C, 1, col-lH), Cnew);
     655       21117 :   return H;
     656             : }
     657             : 
     658             : GEN
     659           0 : hnfadd(GEN H, GEN perm, GEN* ptdep, GEN* ptB, GEN* ptC, /* cf hnfspec */
     660             :        GEN extramat,GEN extraC)
     661             : {
     662           0 :   pari_sp av = avma;
     663           0 :   H = hnfadd_i(H, perm, ptdep, ptB, ptC, ZM_to_zm(extramat), extraC);
     664           0 :   gerepileall(av, 4, ptC, ptdep, ptB, &H); return H;
     665             : }
     666             : 
     667             : /* zero aj = Aij (!= 0)  using  ak = Aik (maybe 0), via linear combination of
     668             :  * A[j] and A[k] of determinant 1. If U != NULL, likewise update its columns */
     669             : static void
     670    24443333 : ZC_elem(GEN aj, GEN ak, GEN A, GEN U, long j, long k)
     671             : {
     672             :   GEN p1,u,v,d;
     673             : 
     674    24443333 :   if (!signe(ak)) {
     675       20245 :     swap(gel(A,j), gel(A,k));
     676       20245 :     if (U) swap(gel(U,j), gel(U,k));
     677    21391863 :     return;
     678             :   }
     679    24423088 :   d = bezout(aj,ak,&u,&v);
     680             :   /* frequent special case (u,v) = (1,0) or (0,1) */
     681    24423142 :   if (!signe(u))
     682             :   { /* ak | aj */
     683    10493563 :     p1 = diviiexact(aj,ak); togglesign(p1);
     684    10492967 :     ZC_lincomb1_inplace(gel(A,j), gel(A,k), p1);
     685    10493381 :     if (U)
     686      533991 :       ZC_lincomb1_inplace(gel(U,j), gel(U,k), p1);
     687    10493426 :     return;
     688             :   }
     689    13929579 :   if (!signe(v))
     690             :   { /* aj | ak */
     691    10858044 :     p1 = diviiexact(ak,aj); togglesign(p1);
     692    10857951 :     ZC_lincomb1_inplace(gel(A,k), gel(A,j), p1);
     693    10857947 :     swap(gel(A,j), gel(A,k));
     694    10857947 :     if (U) {
     695       62466 :       ZC_lincomb1_inplace(gel(U,k), gel(U,j), p1);
     696       62466 :       swap(gel(U,j), gel(U,k));
     697             :     }
     698    10857947 :     return;
     699             :   }
     700             : 
     701     3071535 :   if (!is_pm1(d)) { aj = diviiexact(aj, d); ak = diviiexact(ak, d); }
     702     3071535 :   p1 = gel(A,k); aj = negi(aj); /* NOT togglesign */
     703     3071574 :   gel(A,k) = ZC_lincomb(u,v, gel(A,j),p1);
     704     3071529 :   gel(A,j) = ZC_lincomb(aj,ak, p1,gel(A,j));
     705     3071551 :   if (U)
     706             :   {
     707      184705 :     p1 = gel(U,k);
     708      184705 :     gel(U,k) = ZC_lincomb(u,v, gel(U,j),p1);
     709      184695 :     gel(U,j) = ZC_lincomb(aj,ak, p1,gel(U,j));
     710             :   }
     711             : }
     712             : 
     713             : INLINE int
     714        1001 : is_RgX(GEN a, long v) { return typ(a) == t_POL && varn(a)==v; }
     715             : /* set u,v such that au + bv = gcd(a,b), divide a,b by the gcd */
     716             : static GEN
     717         266 : gbezout_step(GEN *pa, GEN *pb, GEN *pu, GEN *pv, long vx)
     718             : {
     719         266 :   GEN a = *pa, b = *pb, d;
     720         266 :   if (gequal0(a))
     721             :   {
     722           0 :     *pa = gen_0; *pu = gen_0;
     723           0 :     *pb = gen_1; *pv = gen_1; return b;
     724             :   }
     725         266 :   a = is_RgX(a,vx)? RgX_renormalize(a): scalarpol(a, vx);
     726         266 :   b = is_RgX(b,vx)? RgX_renormalize(b): scalarpol(b, vx);
     727         266 :   d = RgX_extgcd(a,b, pu,pv);
     728         266 :   if (degpol(d)) { a = RgX_div(a, d); b = RgX_div(b, d); }
     729         154 :   else if (typ(gel(d,2)) == t_REAL && lg(gel(d,2)) <= 3)
     730             : #if 1
     731             :   { /* possible accuracy problem */
     732           0 :     GEN D = RgX_gcd_simple(a,b);
     733           0 :     if (degpol(D)) {
     734           0 :       D = RgX_normalize(D);
     735           0 :       a = RgX_div(a, D);
     736           0 :       b = RgX_div(b, D);
     737           0 :       d = RgX_extgcd(a,b, pu,pv); /* retry now */
     738           0 :       d = RgX_mul(d, D);
     739             :     }
     740             :   }
     741             : #else
     742             :   { /* less stable */
     743             :     d = RgX_extgcd_simple(a,b, pu,pv);
     744             :     if (degpol(d)) { a = RgX_div(a, d); b = RgX_div(b, d); }
     745             :   }
     746             : #endif
     747         266 :   *pa = a;
     748         266 :   *pb = b; return d;
     749             : }
     750             : static GEN
     751      646216 : col_mul(GEN x, GEN c)
     752             : {
     753      646216 :   if (typ(x) == t_INT)
     754             :   {
     755      646216 :     long s = signe(x);
     756      646216 :     if (!s) return NULL;
     757      488843 :     if (is_pm1(x)) return (s > 0)? c: RgC_neg(c);
     758             :   }
     759       89150 :   return RgC_Rg_mul(c, x);
     760             : }
     761             : static void
     762           0 : do_zero(GEN x)
     763             : {
     764           0 :   long i, lx = lg(x);
     765           0 :   for (i=1; i<lx; i++) gel(x,i) = gen_0;
     766           0 : }
     767             : 
     768             : /* (c1, c2) *= [u,-b; v,a] */
     769             : static void
     770      161554 : update(GEN u, GEN v, GEN a, GEN b, GEN *c1, GEN *c2)
     771             : {
     772             :   GEN p1,p2;
     773             : 
     774      161554 :   u = col_mul(u,*c1);
     775      161554 :   v = col_mul(v,*c2);
     776      161554 :   if (u) p1 = v? gadd(u,v): u;
     777        1528 :   else   p1 = v? v: NULL;
     778             : 
     779      161554 :   a = col_mul(a,*c2);
     780      161554 :   b = col_mul(gneg_i(b),*c1);
     781      161554 :   if (a) p2 = b? RgC_add(a,b): a;
     782           0 :   else   p2 = b? b: NULL;
     783             : 
     784      161554 :   if (!p1) do_zero(*c1); else *c1 = p1;
     785      161554 :   if (!p2) do_zero(*c2); else *c2 = p2;
     786      161554 : }
     787             : 
     788             : /* zero aj = Aij (!= 0)  using  ak = Aik (maybe 0), via linear combination of
     789             :  * A[j] and A[k] of determinant 1. If U != NULL, likewise update its columns */
     790             : static void
     791           7 : RgC_elem(GEN aj, GEN ak, GEN A, GEN V, long j, long k, long li, long vx)
     792             : {
     793           7 :   GEN u,v, d = gbezout_step(&aj, &ak, &u, &v, vx);
     794             :   long l;
     795             :   /* (A[,k], A[,j]) *= [v, -aj; u, ak ] */
     796          14 :   for (l = 1; l < li; l++)
     797             :   {
     798           7 :     GEN t = gadd(gmul(u,gcoeff(A,l,j)), gmul(v,gcoeff(A,l,k)));
     799           7 :     gcoeff(A,l,j) = gsub(gmul(ak,gcoeff(A,l,j)), gmul(aj,gcoeff(A,l,k)));
     800           7 :     gcoeff(A,l,k) = t;
     801             :   }
     802           7 :   gcoeff(A,li,j) = gen_0;
     803           7 :   gcoeff(A,li,k) = d;
     804           7 :   if (V) update(v,u,ak,aj,(GEN*)(V+k),(GEN*)(V+j));
     805           7 : }
     806             : 
     807             : /* reduce A[i,j] mod A[i,j0] for j=j0+1... via column operations */
     808             : static void
     809     1460362 : ZM_reduce(GEN A, GEN U, long i, long j0)
     810             : {
     811     1460362 :   long j, lA = lg(A);
     812     1460362 :   GEN d = gcoeff(A,i,j0);
     813     1460362 :   if (signe(d) < 0)
     814             :   {
     815        5601 :     ZV_neg_inplace(gel(A,j0));
     816        5601 :     if (U) ZV_togglesign(gel(U,j0));
     817        5601 :     d = gcoeff(A,i,j0);
     818             :   }
     819     3586971 :   for (j=j0+1; j<lA; j++)
     820             :   {
     821     2126609 :     GEN q = truedivii(gcoeff(A,i,j), d);
     822     2126609 :     if (!signe(q)) continue;
     823             : 
     824      124049 :     togglesign(q);
     825      124049 :     ZC_lincomb1_inplace(gel(A,j), gel(A,j0), q);
     826      124049 :     if (U) ZC_lincomb1_inplace(gel(U,j), gel(U,j0), q);
     827             :   }
     828     1460362 : }
     829             : 
     830             : /* normalize T as if it were a t_POL in variable v */
     831             : static GEN
     832         210 : normalize_as_RgX(GEN T, long v, GEN *pd)
     833             : {
     834             :   GEN d;
     835         210 :   if (!is_RgX(T,v)) { *pd = T; return gen_1; }
     836         175 :   d = leading_coeff(T);
     837         350 :   while (gequal0(d) || (typ(d) == t_REAL && lg(d) == 3
     838           0 :                         && gexpo(T) - expo(d) > (long)BITS_IN_LONG)) {
     839          14 :     T = normalizepol_lg(T, lg(T)-1);
     840          14 :     if (!signe(T)) { *pd = gen_1; return T; }
     841           0 :     d = leading_coeff(T);
     842             :   }
     843         161 :   if (degpol(T)) T = RgX_Rg_div(T,d); else { d = gel(T,2); T = gen_1; }
     844         161 :   *pd = d; return T;
     845             : }
     846             : /* reduce A[i,j] mod A[i,j0] for j=j0+1... via column operations */
     847             : static void
     848          42 : RgM_reduce(GEN A, GEN U, long i, long j0, long vx)
     849             : {
     850          42 :   long j, lA = lg(A);
     851          42 :   GEN d, T = normalize_as_RgX(gcoeff(A,i,j0), vx, &d);
     852          42 :   if (U && !gequal1(d)) gel(U,j0) = RgC_Rg_div(gel(U,j0), d);
     853          42 :   gcoeff(A,i,j0) = T;
     854             : 
     855          49 :   for (j=j0+1; j<lA; j++)
     856             :   {
     857           7 :     GEN t = gcoeff(A,i,j), q;
     858           7 :     if (gequal0(t)) continue;
     859           7 :     if (T == gen_1)
     860           0 :       q = t;
     861           7 :     else if (is_RgX(t,vx))
     862           7 :       q = RgX_div(t, T);
     863           0 :     else continue;
     864             : 
     865           7 :     if (gequal0(q)) continue;
     866           0 :     gel(A,j) = RgC_sub(gel(A,j), RgC_Rg_mul(gel(A,j0), q));
     867           0 :     if (U) gel(U,j) = RgC_sub(gel(U,j), RgC_Rg_mul(gel(U,j0), q));
     868             :   }
     869          42 : }
     870             : 
     871             : /* A,B square integral in upper HNF, of the same dimension > 0. Return Au
     872             :  * in Z^n (v in Z^n not computed), such that Au + Bv = [1, 0, ..., 0] */
     873             : GEN
     874      169469 : hnfmerge_get_1(GEN A, GEN B)
     875             : {
     876      169469 :   pari_sp av = avma;
     877      169469 :   long j, k, l = lg(A), lb;
     878      169469 :   GEN b, U = cgetg(l + 1, t_MAT), C = cgetg(l + 1, t_VEC);
     879             : 
     880      169434 :   b = gcoeff(B,1,1); lb = lgefint(b);
     881      326332 :   for (j = 1; j < l; j++)
     882             :   {
     883             :     GEN t;
     884      326331 :     long c = j+1;
     885      326331 :     gel(U,j) = col_ei(l-1, j);
     886      326410 :     gel(U,c) = zerocol(l-1); /* dummy */
     887      326426 :     gel(C,j) = vecslice(gel(A,j), 1,j);
     888      326413 :     gel(C,c) = vecslice(gel(B,j), 1,j);
     889      920556 :     for (k = j; k > 0; k--)
     890             :     {
     891      594223 :       t = gcoeff(C,k,c);
     892      594223 :       if (gequal0(t)) continue;
     893      522959 :       setlg(C[c], k+1);
     894      522937 :       ZC_elem(t, gcoeff(C,k,k), C, U, c, k);
     895      522612 :       if (lgefint(gcoeff(C,k,k)) > lb) gel(C,k) = FpC_red(gel(C,k), b);
     896      522612 :       if (j > 4)
     897             :       {
     898       52583 :         GEN u = gel(U,k);
     899             :         long h;
     900      703957 :         for (h=1; h<l; h++)
     901      651374 :           if (lgefint(gel(u,h)) > lb) gel(u,h) = remii(gel(u,h), b);
     902             :       }
     903             :     }
     904      326333 :     if (j == 1)
     905      169430 :       t = gcoeff(C,1,1);
     906             :     else
     907             :     {
     908             :       GEN u;
     909      156903 :       t = bezout(gcoeff(C,1,1), b, &u, NULL); /* >= 0 */
     910      156939 :       if (signe(u) && !equali1(u)) gel(U,1) = ZC_Z_mul(gel(U,1), u);
     911      156940 :       gcoeff(C,1,1) = t;
     912             :     }
     913      326370 :     if (equali1(t)) break;
     914             :   }
     915      169467 :   if (j >= l) return NULL;
     916      169467 :   b = lcmii(gcoeff(A,1,1),b);
     917      169447 :   A = FpC_red(ZM_ZC_mul(A,gel(U,1)), b);
     918      169452 :   return gerepileupto(av, FpC_center(A, b, shifti(b,-1)));
     919             : }
     920             : 
     921             : /* remove the first r columns */
     922             : static void
     923       75988 : remove_0cols(long r, GEN *pA, GEN *pB, long remove)
     924             : {
     925       75988 :   GEN A = *pA, B = *pB;
     926       75988 :   long l = lg(A);
     927       75988 :   A += r; A[0] = evaltyp(t_MAT) | evallg(l-r);
     928       75988 :   if (B && remove == 2) { B += r; B[0] = A[0]; }
     929       75988 :   *pA = A; *pB = B;
     930       75988 : }
     931             : 
     932             : /* Inefficient compared to hnfall. 'remove' = throw away lin.dep columns */
     933             : static GEN
     934       19432 : hnf_i(GEN A, int remove)
     935             : {
     936       19432 :   pari_sp av0 = avma, av;
     937             :   long s, n, m, j, k, li, def, ldef;
     938             : 
     939       19432 :   RgM_dimensions(A, &m, &n);
     940       19432 :   if (!n) return cgetg(1,t_MAT);
     941       19376 :   av = avma;
     942       19376 :   A = RgM_shallowcopy(A);
     943       19376 :   def = n; ldef = (m>n)? m-n: 0;
     944       67123 :   for (li=m; li>ldef; li--)
     945             :   {
     946      114660 :     for (j=def-1; j; j--)
     947             :     {
     948       66913 :       GEN a = gcoeff(A,li,j);
     949       66913 :       if (!signe(a)) continue;
     950             : 
     951             :       /* zero a = Aij  using  b = Aik */
     952       36260 :       k = (j==1)? def: j-1;
     953       36260 :       ZC_elem(a,gcoeff(A,li,k), A,NULL, j,k);
     954       36260 :       if (gc_needed(av,1))
     955             :       {
     956           0 :         if (DEBUGMEM>1) pari_warn(warnmem,"ZM_hnf[1]. li=%ld",li);
     957           0 :         A = gerepilecopy(av, A);
     958             :       }
     959             :     }
     960       47747 :     s = signe(gcoeff(A,li,def));
     961       47747 :     if (s)
     962             :     {
     963       47684 :       if (s < 0) ZV_neg_inplace(gel(A,def));
     964       47684 :       ZM_reduce(A, NULL, li,def);
     965       47684 :       def--;
     966             :     }
     967             :     else
     968          63 :       if (ldef) ldef--;
     969       47747 :     if (gc_needed(av,1))
     970             :     {
     971           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"ZM_hnf[2]. li=%ld",li);
     972           0 :       A = gerepilecopy(av, A);
     973             :     }
     974             :   }
     975             :   /* rank A = n - def */
     976       19376 :   if (remove) { GEN B = NULL; remove_0cols(def, &A, &B, remove); }
     977       19376 :   return gerepileupto(av0, ZM_copy(A));
     978             : }
     979             : 
     980             : GEN
     981       23215 : ZM_hnf(GEN x) { return lg(x) > 8? ZM_hnfall(x, NULL, 1): hnf_i(x, 1); }
     982             : 
     983             : /* u*z[1..k] mod p, in place */
     984             : static void
     985     2948184 : FpV_Fp_mul_part_ip(GEN z, GEN u, GEN p, long k)
     986             : {
     987             :   long i;
     988     2948184 :   if (is_pm1(u)) {
     989       75817 :     if (signe(u) > 0) {
     990      226183 :       for (i = 1; i <= k; i++)
     991      181142 :         if (signe(gel(z,i))) gel(z,i) = modii(gel(z,i), p);
     992             :     } else {
     993       77337 :       for (i = 1; i <= k; i++)
     994       46561 :         if (signe(gel(z,i))) gel(z,i) = modii(negi(gel(z,i)), p);
     995             :     }
     996             :   }
     997             :   else {
     998    10704695 :     for (i = 1; i <= k; i++)
     999     7832518 :       if (signe(gel(z,i))) gel(z,i) = Fp_mul(u,gel(z,i), p);
    1000             :   }
    1001     2947994 : }
    1002             : static void
    1003    21881294 : FpV_red_part_ipvec(GEN z, GEN p, long k)
    1004             : {
    1005             :   long i;
    1006    21881294 :   for (i = 1; i <= k; i++) gel(z,i) = modii(gel(z,i), gel(p,i));
    1007    21881282 : }
    1008             : 
    1009             : /* return x * U, in echelon form (mod p^m), where (det(U),p) = 1.
    1010             :  * If early_abort is set, return NULL as soon as one pivot is 0 (mod p^m) */
    1011             : GEN
    1012       55461 : ZpM_echelon(GEN x, long early_abort, GEN p, GEN pm)
    1013             : {
    1014       55461 :   pari_sp av0 = avma, av;
    1015             :   long m, li, co, i, j, k, def, ldef;
    1016             : 
    1017       55461 :   co = lg(x); if (co == 1) return cgetg(1,t_MAT);
    1018       55461 :   li = lgcols(x);
    1019       55461 :   av = avma;
    1020       55461 :   x = RgM_shallowcopy(x);
    1021       55461 :   m = Z_pval(pm, p);
    1022             : 
    1023       55461 :   ldef = (li > co)? li - co: 0;
    1024      494217 :   for (def = co-1,i = li-1; i > ldef; i--)
    1025             :   {
    1026      439336 :     long vmin = LONG_MAX, kmin = 0;
    1027      439336 :     GEN umin = gen_0, pvmin, q;
    1028     3915248 :     for (k = 1; k <= def; k++)
    1029             :     {
    1030     3562487 :       GEN u = gcoeff(x,i,k);
    1031             :       long v;
    1032     3562487 :       if (!signe(u)) continue;
    1033     1681609 :       v = Z_pvalrem(u, p, &u);
    1034     1681609 :       if (v >= m) gcoeff(x,i,k) = gen_0;
    1035     1098395 :       else if (v < vmin) {
    1036      417884 :         vmin = v; kmin = k; umin = u;
    1037      417884 :         if (!vmin) break;
    1038             :       }
    1039             :     }
    1040      439336 :     if (!kmin)
    1041             :     {
    1042      180859 :       if (early_abort) return NULL;
    1043      180279 :       gcoeff(x,i,def) = gen_0;
    1044      180279 :       ldef--;
    1045      180279 :       if (ldef < 0) ldef = 0;
    1046      180279 :       continue;
    1047             :     }
    1048      258477 :     if (kmin != def) swap(gel(x,def), gel(x,kmin));
    1049      258477 :     q = vmin? powiu(p, m-vmin): pm;
    1050             :     /* pivot has valuation vmin */
    1051      258477 :     umin = modii(umin, q);
    1052      258477 :     if (!equali1(umin))
    1053      197165 :       FpV_Fp_mul_part_ip(gel(x,def), Fp_inv(umin,q), pm, i-1);
    1054      258477 :     gcoeff(x, i, def) = pvmin = powiu(p, vmin);
    1055     2860871 :     for (j = def-1; j; j--)
    1056             :     { /* zero x[i, 1..def-1] using x[i,def] = pvmin */
    1057     2602394 :       GEN t, a = gcoeff(x,i,j) = modii(gcoeff(x,i,j), pm);
    1058     2602394 :       if (!signe(a)) continue;
    1059             : 
    1060     1271970 :       t = diviiexact(a, pvmin); togglesign(t);
    1061     1271970 :       ZC_lincomb1_inplace(gel(x,j), gel(x,def), t);
    1062     1271970 :       if (gc_needed(av,1))
    1063             :       {
    1064          14 :         if (DEBUGMEM>1) pari_warn(warnmem,"ZpM_echelon. i=%ld",i);
    1065          14 :         x = gerepilecopy(av, x); pvmin = gcoeff(x,i,def);
    1066             :       }
    1067             :     }
    1068      258477 :     def--;
    1069             :   }
    1070       54881 :   if (co > li)
    1071             :   {
    1072           0 :     x += co - li;
    1073           0 :     x[0] = evaltyp(t_MAT) | evallg(li);
    1074             :   }
    1075       54881 :   return gerepilecopy(av0, x);
    1076             : }
    1077             : GEN
    1078      974025 : zlm_echelon(GEN x, long early_abort, ulong p, ulong pm)
    1079             : {
    1080      974025 :   pari_sp av0 = avma;
    1081             :   long li, co, i, j, k, def, ldef;
    1082             :   ulong m;
    1083             : 
    1084      974025 :   co = lg(x); if (co == 1) return cgetg(1,t_MAT);
    1085      974025 :   li = lgcols(x);
    1086      974026 :   x = Flm_copy(x);
    1087      974032 :   m = u_lval(pm, p);
    1088             : 
    1089      974032 :   ldef = (li > co)? li - co: 0;
    1090     5132530 :   for (def = co-1,i = li-1; i > ldef; i--)
    1091             :   {
    1092     4230588 :     long vmin = LONG_MAX, kmin = 0;
    1093     4230588 :     ulong umin = 0, pvmin, q;
    1094    16328140 :     for (k = 1; k <= def; k++)
    1095             :     {
    1096    14572981 :       ulong u = ucoeff(x,i,k);
    1097             :       long v;
    1098    14572981 :       if (!u) continue;
    1099     7119191 :       v = u_lvalrem(u, p, &u);
    1100     7119212 :       if (v >= (long) m) ucoeff(x,i,k) = 0;
    1101     7119212 :       else if (v < vmin) {
    1102     4329919 :         vmin = v; kmin = k; umin = u;
    1103     4329919 :         if (!vmin) break;
    1104             :       }
    1105             :     }
    1106     4230609 :     if (!kmin)
    1107             :     {
    1108      283663 :       if (early_abort) return NULL;
    1109      211560 :       ucoeff(x,i,def) = 0;
    1110      211560 :       ldef--;
    1111      211560 :       if (ldef < 0) ldef = 0;
    1112      211560 :       continue;
    1113             :     }
    1114     3946946 :     if (kmin != def) swap(gel(x,def), gel(x,kmin));
    1115     3946946 :     q = vmin? upowuu(p, m-vmin): pm;
    1116             :     /* pivot has valuation vmin */
    1117     3946953 :     umin %= q;
    1118     3946953 :     if (umin != 1)
    1119     2078991 :       Flv_Fl_mul_part_inplace(gel(x,def), Fl_inv(umin,q), pm, i-1);
    1120     3946944 :     ucoeff(x, i, def) = pvmin = upowuu(p, vmin);
    1121    23601720 :     for (j = def-1; j; j--)
    1122             :     { /* zero x[i, 1..def-1] using x[i,def] = pvmin */
    1123    19654782 :       ulong t, a = ucoeff(x,i,j);
    1124    19654782 :       if (!a) continue;
    1125             : 
    1126    11494199 :       t = Fl_neg(a / pvmin, q);
    1127    11494204 :       Flc_lincomb1_inplace(gel(x,j), gel(x,def), t, pm);
    1128             :     }
    1129     3946938 :     def--;
    1130             :   }
    1131      901942 :   if (co > li)
    1132             :   {
    1133           0 :     x += co - li;
    1134           0 :     x[0] = evaltyp(t_MAT) | evallg(li);
    1135             :   }
    1136      901942 :   return gerepilecopy(av0, x);
    1137             : }
    1138             : 
    1139             : static int
    1140       18494 : ZV_allequal(GEN v)
    1141             : {
    1142       18494 :   long i, l = lg(v);
    1143       18494 :   if (l > 1)
    1144             :   {
    1145       18494 :     GEN x = gel(v,1);
    1146       18494 :     for (i = 2; i < l; i++) if (!equalii(x,gel(v,i))) return 0;
    1147             :   }
    1148       10262 :   return 1;
    1149             : }
    1150             : /* compute optimal D for hnfmod: x upper triangular */
    1151             : static GEN
    1152     3189431 : optimal_D(GEN x, GEN D)
    1153             : {
    1154     3189431 :   long i, n = nbrows(x);
    1155     3189431 :   GEN C = shallowcopy(D);
    1156     3189448 :   gel(C,1) = gcoeff(x,1,1);
    1157     3552380 :   for (i = 2; i < n; i++)
    1158             :   {
    1159     1214286 :     GEN c = mulii(gel(C,i-1), gcoeff(x,i,i));
    1160     1214271 :     if (signe(c) < 0) togglesign(c);
    1161     1214272 :     if (cmpii(c, gel(D,i)) >= 0) break;
    1162      362932 :     gel(C,i) = c;
    1163             :   }
    1164     3189445 :   return C;
    1165             : }
    1166             : 
    1167             : /* D = multiple of det x (usually detint(x)) or vector of positive moduli
    1168             :  * (compute hnf(x | D))
    1169             :  * flag & hnf_MODID: reduce mod D * matid [ otherwise as above ].
    1170             :  * flag & hnf_PART: don't reduce once diagonal is known
    1171             :  * flag & hnf_CENTER: centermod HNF (2|x[i,j]| <] x[i,i]) */
    1172             : GEN
    1173     3189525 : ZM_hnfmodall_i(GEN x, GEN D, long flag)
    1174             : {
    1175             :   pari_sp av;
    1176     3189525 :   const long center = (flag & hnf_CENTER);
    1177             :   long moddiag, modsame, nli, li, co, i, j, k, def, ldef;
    1178             :   GEN u, LDM;
    1179             : 
    1180     3189525 :   co = lg(x);
    1181     3189525 :   if (co == 1)
    1182             :   {
    1183          56 :     if (typ(D) == t_INT || lg(D) == 1) return cgetg(1,t_MAT);
    1184          14 :     x = diagonal_shallow(D); /* handle flags properly */
    1185          14 :     co = lg(x);
    1186             :   }
    1187     3189483 :   li = lgcols(x);
    1188     3189488 :   if (li == 1)
    1189             :   {
    1190          21 :     if (typ(D) != t_INT && lg(D) != li) pari_err_DIM("ZM_hnfmod");
    1191          21 :     return cgetg(1,t_MAT);
    1192             :   }
    1193     3189467 :   nli = li - 1;
    1194     3189467 :   modsame = typ(D)==t_INT;
    1195     3189467 :   if (!modsame)
    1196             :   {
    1197       18501 :     if (lg(D) != li) pari_err_DIM("ZM_hnfmod");
    1198       18494 :     if (ZV_allequal(D)) { modsame = 1; D = gel(D,1); }
    1199             :   }
    1200     3189460 :   moddiag = (flag & hnf_MODID) || !modsame;
    1201             :   /* modsame: triangularize mod fixed d*Id;
    1202             :    * moddiag: modulo diagonal matrix, else modulo multiple of determinant */
    1203             : 
    1204     3189460 :   if (modsame)
    1205             :   {
    1206     3181228 :     LDM = const_vecsmall(nli, 2*lgefint(D)-2);
    1207     3181238 :     D = const_vec(nli,D);
    1208             :   }
    1209             :   else
    1210             :   {
    1211        8232 :     LDM = cgetg(li, t_VECSMALL);
    1212        8232 :     for (i=1; i<li; i++) LDM[i] = lgefint(gel(D,i));
    1213             :   }
    1214     3189475 :   av = avma;
    1215     3189475 :   x = RgM_shallowcopy(x);
    1216             : 
    1217     3189479 :   ldef = 0;
    1218     3189479 :   if (li > co)
    1219             :   {
    1220       48531 :     ldef = li - co;
    1221       48531 :     if (!moddiag)
    1222           7 :       pari_err_DOMAIN("ZM_hnfmod","nb lines",">", strtoGENstr("nb columns"), x);
    1223             :   }
    1224    11980268 :   for (def = co-1,i = nli; i > ldef; i--,def--)
    1225             :   {
    1226     8790829 :     GEN d = gel(D,i);
    1227     8790829 :     long add_N = modsame;
    1228    43246242 :     for (j = 1; j < def; j++)
    1229             :     {
    1230    34455417 :       GEN p1, p2, b, a = gcoeff(x,i,j) = remii(gcoeff(x,i,j), d);
    1231    54537275 :       if (!signe(a)) continue;
    1232             : 
    1233    23515331 :       k = j+1;
    1234    23515331 :       b = gcoeff(x,i,k) = remii(gcoeff(x,i,k), d);
    1235    23515333 :       if (!signe(b)) { swap(gel(x,j), gel(x,k)); continue; }
    1236    14373529 :       if (add_N)
    1237             :       { /* ensure the moving pivot on row i divides d from now on */
    1238     4566851 :         add_N = 0;
    1239     4566851 :         if (!equali1(a))
    1240             :         { /* x[j] *= u; after this, a = x[i,j] | d */
    1241     3868863 :           GEN u = Fp_invgen(a, d, &a);
    1242             :           long t;
    1243     3868863 :           p1 = gel(x,j);
    1244     3868863 :           for (t = 1; t < i; t++) gel(p1,t) = mulii(gel(p1,t), u);
    1245     3868863 :           FpV_red_part_ipvec(p1, D, i-1);
    1246     3868863 :           gel(p1,i) = a;
    1247     3868863 :           if (2*lg(a) < lg(b))
    1248             :           { /* reduce x[i,k] mod x[i,j]: helps ZC_elem */
    1249        9577 :             GEN r, q = dvmdii(b, a, &r);
    1250        9577 :             togglesign(q);
    1251        9577 :             ZC_lincomb1_inplace_i(gel(x,k), gel(x,j), q, i-1);
    1252        9577 :             FpV_red_part_ipvec(gel(x,k), D, i-1);
    1253        9577 :             gcoeff(x,i,k) = b = r;
    1254             :           }
    1255             :         }
    1256             :       }
    1257    14373527 :       ZC_elem(a,b, x, NULL, j,k);
    1258    14373537 :       p1 = gel(x,j);
    1259    14373537 :       p2 = gel(x,k);
    1260             :       /* prevent coeffs explosion */
    1261    67336379 :       for (k = 1; k < i; k++)
    1262             :       {
    1263    52962840 :         if (lgefint(gel(p1,k)) > LDM[k]) gel(p1,k) = remii(gel(p1,k),gel(D,k));
    1264    52962842 :         if (lgefint(gel(p2,k)) > LDM[k]) gel(p2,k) = remii(gel(p2,k),gel(D,k));
    1265             :       }
    1266             :     }
    1267     8790825 :     if (gc_needed(av,2))
    1268             :     {
    1269           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"ZM_hnfmod[1]. i=%ld",i);
    1270           0 :       x = gerepilecopy(av, x);
    1271             :     }
    1272     8790825 :     if (moddiag && !signe(gcoeff(x,i,def)))
    1273             :     { /* missing pivot on line i, insert column */
    1274      392841 :       GEN a = cgetg(co + 1, t_MAT);
    1275      392841 :       for (k = 1; k <= def; k++) gel(a,k) = gel(x,k);
    1276      392841 :       gel(a,k++) = Rg_col_ei(gel(D,i), nli, i);
    1277      392841 :       for (     ; k <= co;  k++) gel(a,k) = gel(x,k-1);
    1278      392841 :       ldef--; if (ldef < 0) ldef = 0;
    1279      392841 :       co++; def++; x = a;
    1280             :     }
    1281             :   }
    1282     3189439 :   if (co < li)
    1283             :   { /* implies moddiag, add missing diag(D) components */
    1284       29225 :     GEN a = cgetg(li+1, t_MAT);
    1285       29225 :     for (k = 1; k <= li-co; k++) gel(a,k) = Rg_col_ei(gel(D,k), nli, k);
    1286       29225 :     for (i = 1; i < co; i++) gel(a,k-1+i) = gel(x,i);
    1287       29225 :     gel(a,li) = zerocol(nli); x = a;
    1288             :   }
    1289             :   else
    1290             :   {
    1291     3160214 :     x += co - li;
    1292     3160214 :     x[0] = evaltyp(t_MAT) | evallg(li); /* kill 0 columns */
    1293     3160210 :     if (moddiag) x = shallowconcat(x, zerocol(nli));
    1294             :   }
    1295     3189459 :   if (moddiag)
    1296             :   { /* x[li]: extra column, an accumulator discarded at the end */
    1297             :     GEN D2;
    1298     3126144 :     gcoeff(x,1,1) = gcdii(gcoeff(x,1,1), gel(D,1));
    1299     3126117 :     D2 = optimal_D(x,D);
    1300             :     /* add up missing diag(D) components */
    1301    11805803 :     for (i = nli; i > 0; i--)
    1302             :     {
    1303     8679665 :       gcoeff(x, i, li) = gel(D,i);
    1304    31956802 :       for (j = i; j > 0; j--)
    1305             :       {
    1306    23277143 :         GEN a = gcoeff(x, j, li);
    1307    23277143 :         if (!signe(a)) continue;
    1308     9001500 :         ZC_elem(a, gcoeff(x,j,j), x, NULL, li,j);
    1309     9001513 :         FpV_red_part_ipvec(gel(x,li), D, j-1);
    1310     9001472 :         FpV_red_part_ipvec(gel(x,j),  D, j-1);
    1311             :       }
    1312     8679659 :       if (gc_needed(av,1))
    1313             :       {
    1314           6 :         if (DEBUGMEM>1) pari_warn(warnmem,"ZM_hnfmod[2]. i=%ld", i);
    1315           6 :         gerepileall(av, 2, &x, &D2);
    1316             :       }
    1317             :     }
    1318     3126138 :     D = D2;
    1319             :   }
    1320             :   else
    1321             :   {
    1322       63315 :     GEN b = gel(D,1);
    1323      237244 :     for (i = nli; i > 0; i--)
    1324             :     {
    1325      173929 :       GEN d = bezout(gcoeff(x,i,i),b, &u,NULL);
    1326      173929 :       gcoeff(x,i,i) = d;
    1327      173929 :       FpV_Fp_mul_part_ip(gel(x,i), u, b, i-1);
    1328      173929 :       if (i > 1) b = diviiexact(b,d);
    1329             :     }
    1330       63315 :     D = optimal_D(x,D);
    1331             :   }
    1332     3189453 :   x[0] = evaltyp(t_MAT) | evallg(li); /* kill 0 columns / discard accumulator */
    1333     3189500 :   if (flag & hnf_PART) return x;
    1334             : 
    1335    12043057 :   for (i = nli; i > 0; i--)
    1336             :   {
    1337     8853631 :     GEN diag = gcoeff(x,i,i);
    1338     8853631 :     if (signe(diag) < 0) { gel(x,i) = ZC_neg(gel(x,i)); diag = gcoeff(x,i,i); }
    1339     8853608 :     if (i != nli)
    1340     5664195 :       for (j = 1; j < i; j++) gcoeff(x,j,i) = remii(gcoeff(x,j,i), gel(D,j));
    1341    23745284 :     for (j = i+1; j < li; j++)
    1342             :     {
    1343    14891731 :       GEN b = gcoeff(x,i,j) = remii(gcoeff(x,i,j), gel(D,i));
    1344    14891702 :       GEN r, q = truedvmdii(b, diag, &r);
    1345             :       /* ensure -diag/2 <= r < diag/2 */
    1346    14891679 :       if (center && signe(r) && abscmpii(shifti(r,1),diag) >= 0)
    1347      294560 :       { r = subii(r,diag); q = addiu(q,1); }
    1348    14891678 :       if (!signe(q)) continue;
    1349     6148698 :       togglesign(q);
    1350     6148698 :       ZC_lincomb1_inplace_i(gel(x,j), gel(x,i), q, i-1);
    1351     6148698 :       gcoeff(x,i,j) = r;
    1352             :     }
    1353     8853553 :     if (gc_needed(av,1))
    1354             :     {
    1355           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"ZM_hnfmod[3]. i=%ld", i);
    1356           0 :       gerepileall(av, 2, &x, &D);
    1357             :     }
    1358             :   }
    1359     3189426 :   return x;
    1360             : }
    1361             : GEN
    1362     3181677 : ZM_hnfmodall(GEN x, GEN dm, long flag)
    1363             : {
    1364     3181677 :   pari_sp av = avma;
    1365     3181677 :   return gerepilecopy(av, ZM_hnfmodall_i(x, dm, flag));
    1366             : }
    1367             : GEN
    1368       63308 : ZM_hnfmod(GEN x, GEN d) { return ZM_hnfmodall(x,d,0); }
    1369             : GEN
    1370     2863524 : ZM_hnfmodid(GEN x, GEN d)
    1371     2863524 : { return ZM_hnfmodall(x,d,hnf_MODID); }
    1372             : /* return the column echelon form of x with 1's as pivots,
    1373             :  * P contains the row indices containing the pivots in increasing order */
    1374             : static GEN
    1375     1774137 : FpM_echelon(GEN x, GEN *pP, GEN p)
    1376             : {
    1377             :   pari_sp av;
    1378             :   long iP, li, co, i, j, k, def, ldef;
    1379             :   GEN P;
    1380             : 
    1381     1774137 :   co = lg(x); if (co == 1) { *pP = cgetg(1,t_VECSMALL); return cgetg(1,t_MAT); }
    1382     1774137 :   li = lgcols(x);
    1383     1774135 :   iP = 1;
    1384     1774135 :   *pP = P = cgetg(li, t_VECSMALL);
    1385     1774136 :   av = avma;
    1386     1774136 :   x = FpM_red(x, p);
    1387             : 
    1388     1773945 :   ldef = (li > co)? li - co: 0;
    1389     6539412 :   for (def = co-1,i = li-1; i > ldef; i--)
    1390             :   {
    1391     4765274 :     GEN u = NULL;
    1392     7297549 :     for (k = def; k; k--)
    1393             :     {
    1394     5407660 :       u = gcoeff(x,i,k);
    1395     5407660 :       if (signe(u)) break;
    1396             :     }
    1397     4765274 :     if (!k)
    1398             :     {
    1399     1889858 :       if (--ldef < 0) ldef = 0;
    1400     1889858 :       continue;
    1401             :     }
    1402     2875416 :     P[iP++] = i;
    1403     2875416 :     if (k != def) swap(gel(x,def), gel(x,k));
    1404     2875416 :     if (!equali1(u))
    1405     2577239 :       FpV_Fp_mul_part_ip(gel(x,def), Fp_inv(u,p), p, i-1);
    1406     2875204 :     gcoeff(x, i, def) = gen_1;
    1407     9484376 :     for (j = def-1; j; j--)
    1408             :     { /* zero x[i, 1..def-1] using x[i,def] = 1*/
    1409     6608774 :       GEN xj = gel(x,j), u = gel(xj,i);
    1410     6608774 :       if (!signe(u)) continue;
    1411             : 
    1412     4194626 :       ZC_lincomb1_inplace(xj, gel(x,def), negi(u));
    1413     4195025 :       for (k = 1; k < i; k++) gel(xj,k) = modii(gel(xj,k), p);
    1414             :     }
    1415     2875602 :     if (gc_needed(av,2))
    1416             :     {
    1417           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"FpM_echelon. i=%ld",i);
    1418           0 :       x = gerepilecopy(av, x);
    1419             :     }
    1420     2875609 :     def--;
    1421             :   }
    1422             :   /* rank = iP - 1 */
    1423     1774138 :   setlg(P, iP); vecsmall_sort(P);
    1424     1774124 :   if (co > iP) x += co - iP;
    1425     1774124 :   x[0] = evaltyp(t_MAT) | evallg(iP);
    1426     1774120 :   return x;
    1427             : }
    1428             : /* given x square of maximal rank with 1 or p on diagonal from hnfmodid
    1429             :  * (=> a column containing p has its other entries at 0 ), return the HNF */
    1430             : static GEN
    1431     1774157 : FpM_hnfend(pari_sp av, GEN x, GEN p)
    1432             : {
    1433     1774157 :   long i, l = lgcols(x);
    1434     6539556 :   for (i = l-1; i > 0; i--)
    1435             :   {
    1436     4765547 :     GEN diag = gcoeff(x,i,i);
    1437             :     long j;
    1438     4765547 :     if (is_pm1(diag))
    1439     6192882 :       for (j = i+1; j < l; j++)
    1440             :       {
    1441     3317204 :         GEN xj = gel(x,j), b = gel(xj,i);
    1442             :         long k;
    1443     3317204 :         if (!signe(b)) continue;
    1444     2488116 :         ZC_lincomb1_inplace(xj, gel(x,i), negi(b));
    1445    11297700 :         for (k=1; k<i; k++)
    1446     8809642 :           if (lgefint(gel(xj,k)) > 3) gel(xj,k) = remii(gel(xj,k), p);
    1447             :       }
    1448             :     else
    1449     1889897 :       for (j = i+1; j < l; j++) gcoeff(x,i,j) = modii(gcoeff(x,i,j), p);
    1450     4765437 :     if (gc_needed(av,2))
    1451             :     {
    1452           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"FpM_hnfend. i=%ld",i);
    1453           0 :       x = gerepilecopy(av, x);
    1454             :     }
    1455             :   }
    1456     1774009 :   return x;
    1457             : }
    1458             : GEN
    1459     1774142 : ZM_hnfmodprime(GEN x, GEN p)
    1460             : {
    1461     1774142 :   pari_sp av = avma;
    1462             :   GEN P, y;
    1463             :   long l, lP, i;
    1464     1774142 :   if (lg(x) == 1) return cgetg(1, t_MAT);
    1465     1774142 :   l = lgcols(x);
    1466     1774137 :   x = FpM_echelon(x, &P, p);
    1467     1774113 :   lP = lg(P); /* rank = lP-1 */
    1468     1774113 :   if (lP == l) { set_avma(av); return matid(l-1); }
    1469     1774113 :   y = scalarmat_shallow(p, l-1);
    1470     1774155 :   for (i = 1; i < lP; i++) gel(y,P[i]) = gel(x,i);
    1471     1774155 :   return gerepilecopy(av, FpM_hnfend(av,y,p));
    1472             : }
    1473             : 
    1474             : static GEN
    1475      253789 : allhnfmod(GEN x, GEN dm, int flag)
    1476             : {
    1477      253789 :   if (typ(x)!=t_MAT) pari_err_TYPE("allhnfmod",x);
    1478      253789 :   RgM_check_ZM(x, "allhnfmod");
    1479      253782 :   if (isintzero(dm)) return ZM_hnf(x);
    1480      253782 :   return ZM_hnfmodall(x, dm, flag);
    1481             : }
    1482             : GEN
    1483          14 : hnfmod(GEN x, GEN d)
    1484             : {
    1485          14 :   if (typ(d) != t_INT) pari_err_TYPE("mathnfmod",d);
    1486          14 :   return allhnfmod(x, d, 0);
    1487             : }
    1488             : GEN
    1489      253775 : hnfmodid(GEN x, GEN d)
    1490             : {
    1491      253775 :   switch(typ(d))
    1492             :   {
    1493      252459 :     case t_INT: break;
    1494             :     case t_VEC: case t_COL:
    1495        1316 :       if (RgV_is_ZV(d)) break;
    1496           0 :     default: pari_err_TYPE("mathnfmodid",d);
    1497             :   }
    1498      253775 :   return allhnfmod(x, d, hnf_MODID);
    1499             : }
    1500             : 
    1501             : /* M a ZM in HNF. Normalize with *centered* residues */
    1502             : GEN
    1503        2092 : ZM_hnfcenter(GEN M)
    1504             : {
    1505        2092 :   long i, j, k, N = lg(M)-1;
    1506        2092 :   pari_sp av = avma;
    1507             : 
    1508        8017 :   for (j=N-1; j>0; j--) /* skip last line */
    1509             :   {
    1510        5925 :     GEN Mj = gel(M,j), a = gel(Mj,j);
    1511       31073 :     for (k = j+1; k <= N; k++)
    1512             :     {
    1513       25148 :       GEN Mk = gel(M,k), q = diviiround(gel(Mk,j), a);
    1514       25148 :       long s = signe(q);
    1515       25148 :       if (!s) continue;
    1516       16240 :       if (is_pm1(q))
    1517             :       {
    1518        3458 :         if (s < 0)
    1519         478 :           for (i = 1; i <= j; i++) gel(Mk,i) = addii(gel(Mk,i), gel(Mj,i));
    1520             :         else
    1521        2980 :           for (i = 1; i <= j; i++) gel(Mk,i) = subii(gel(Mk,i), gel(Mj,i));
    1522             :       }
    1523             :       else
    1524       12782 :         for (i = 1; i <= j; i++) gel(Mk,i) = subii(gel(Mk,i), mulii(q,gel(Mj,i)));
    1525       16240 :       if (gc_needed(av,1))
    1526             :       {
    1527           0 :         if (DEBUGMEM>1) pari_warn(warnmem,"ZM_hnfcenter, j = %ld",j);
    1528           0 :         M = gerepilecopy(av, M);
    1529             :       }
    1530             :     }
    1531             :   }
    1532        2092 :   return M;
    1533             : }
    1534             : 
    1535             : /***********************************************************************/
    1536             : /*                                                                     */
    1537             : /*                 HNFLLL (Havas, Majewski, Mathews)                   */
    1538             : /*                                                                     */
    1539             : /***********************************************************************/
    1540             : 
    1541             : static void
    1542      573198 : Minus(long j, GEN lambda)
    1543             : {
    1544      573198 :   long k, n = lg(lambda);
    1545             : 
    1546      573198 :   for (k=1  ; k<j; k++) togglesign_safe(&gcoeff(lambda,k,j));
    1547      573198 :   for (k=j+1; k<n; k++) togglesign_safe(&gcoeff(lambda,j,k));
    1548      573198 : }
    1549             : 
    1550             : /* index of first non-zero entry */
    1551             : static long
    1552    46568870 : findi(GEN M)
    1553             : {
    1554    46568870 :   long i, n = lg(M);
    1555   439874716 :   for (i=1; i<n; i++)
    1556   414907066 :     if (signe(gel(M,i))) return i;
    1557    24967650 :   return 0;
    1558             : }
    1559             : 
    1560             : static long
    1561    46568870 : findi_normalize(GEN Aj, GEN B, long j, GEN lambda)
    1562             : {
    1563    46568870 :   long r = findi(Aj);
    1564    46568870 :   if (r && signe(gel(Aj,r)) < 0)
    1565             :   {
    1566      573198 :     ZV_togglesign(Aj); if (B) ZV_togglesign(gel(B,j));
    1567      573198 :     Minus(j,lambda);
    1568             :   }
    1569    46568870 :   return r;
    1570             : }
    1571             : 
    1572             : static void
    1573    23283784 : reduce2(GEN A, GEN B, long k, long j, long *row0, long *row1, GEN lambda, GEN D)
    1574             : {
    1575             :   GEN q;
    1576             :   long i;
    1577             : 
    1578    23283784 :   *row0 = findi_normalize(gel(A,j), B,j,lambda);
    1579    23283784 :   *row1 = findi_normalize(gel(A,k), B,k,lambda);
    1580    23283784 :   if (*row0)
    1581     7886100 :     q = truedivii(gcoeff(A,*row0,k), gcoeff(A,*row0,j));
    1582    15397684 :   else if (abscmpii(shifti(gcoeff(lambda,j,k), 1), gel(D,j)) > 0)
    1583     3364089 :     q = diviiround(gcoeff(lambda,j,k), gel(D,j));
    1584             :   else
    1585    12033595 :     return;
    1586             : 
    1587    11250189 :   if (signe(q))
    1588             :   {
    1589     6319840 :     GEN Lk = gel(lambda,k), Lj = gel(lambda,j);
    1590     6319840 :     togglesign_safe(&q);
    1591     6319840 :     if (*row0) ZC_lincomb1_inplace(gel(A,k),gel(A,j),q);
    1592     6319840 :     if (B) ZC_lincomb1_inplace(gel(B,k),gel(B,j),q);
    1593     6319840 :     gel(Lk,j) = addmulii(gel(Lk,j), q, gel(D,j));
    1594     6319840 :     if (is_pm1(q))
    1595             :     {
    1596     2057334 :       if (signe(q) > 0)
    1597             :       {
    1598     9070280 :         for (i=1; i<j; i++)
    1599     8298896 :           if (signe(gel(Lj,i))) gel(Lk,i) = addii(gel(Lk,i), gel(Lj,i));
    1600             :       }
    1601             :       else
    1602             :       {
    1603    10576531 :         for (i=1; i<j; i++)
    1604     9290581 :           if (signe(gel(Lj,i))) gel(Lk,i) = subii(gel(Lk,i), gel(Lj,i));
    1605             :       }
    1606             :     }
    1607             :     else
    1608             :     {
    1609    29864638 :       for (i=1; i<j; i++)
    1610    25602132 :         if (signe(gel(Lj,i))) gel(Lk,i) = addmulii(gel(Lk,i), q, gel(Lj,i));
    1611             :     }
    1612             :   }
    1613             : }
    1614             : 
    1615             : static void
    1616     3106272 : hnfswap(GEN A, GEN B, long k, GEN lambda, GEN D)
    1617             : {
    1618     3106272 :   GEN t, p1, p2, Lk = gel(lambda,k);
    1619     3106272 :   long i,j,n = lg(A);
    1620             : 
    1621     3106272 :   swap(gel(A,k), gel(A,k-1));
    1622     3106272 :   if (B) swap(gel(B,k), gel(B,k-1));
    1623     3106272 :   for (j=k-2; j; j--) swap(gcoeff(lambda,j,k-1), gel(Lk,j));
    1624    41535288 :   for (i=k+1; i<n; i++)
    1625             :   {
    1626    38429016 :     GEN Li = gel(lambda,i);
    1627    38429016 :     p1 = mulii(gel(Li,k-1), gel(D,k));
    1628    38429016 :     p2 = mulii(gel(Li,k), gel(Lk,k-1));
    1629    38429016 :     t = subii(p1,p2);
    1630             : 
    1631    38429016 :     p1 = mulii(gel(Li,k), gel(D,k-2));
    1632    38429016 :     p2 = mulii(gel(Li,k-1), gel(Lk,k-1));
    1633    38429016 :     gel(Li,k-1) = diviiexact(addii(p1,p2), gel(D,k-1));
    1634    38429016 :     gel(Li,k) = diviiexact(t, gel(D,k-1));
    1635             :   }
    1636     3106272 :   p1 = mulii(gel(D,k-2), gel(D,k));
    1637     3106272 :   p2 = sqri(gel(Lk,k-1));
    1638     3106272 :   gel(D,k-1) = diviiexact(addii(p1,p2), gel(D,k-1));
    1639     3106272 : }
    1640             : 
    1641             : /* reverse row order in matrix A, IN PLACE */
    1642             : static GEN
    1643      124220 : reverse_rows(GEN A)
    1644             : {
    1645      124220 :   long i, j, h, n = lg(A);
    1646      124220 :   if (n == 1) return A;
    1647      124220 :   h = lgcols(A);
    1648     1035594 :   for (j=1; j<n; j++)
    1649             :   {
    1650      911374 :     GEN c = gel(A,j);
    1651             :     /* start at (h-1) >>1 : if h = 2i even, no need to swap c[i] and itself */
    1652      911374 :     for (i=(h-1)>>1; i; i--) swap(gel(c,i), gel(c,h-i));
    1653             :   }
    1654      124220 :   return A;
    1655             : }
    1656             : /* decide whether to swap */
    1657             : static int
    1658     1569010 : must_swap(long k, GEN lambda, GEN D)
    1659             : {
    1660     1569010 :   pari_sp av = avma;
    1661     1569010 :   GEN z = addii(mulii(gel(D,k-2),gel(D,k)), sqri(gcoeff(lambda,k-1,k)));
    1662     1569010 :   return gc_bool(av, cmpii(z, sqri(gel(D,k-1))) < 0);
    1663             : }
    1664             : 
    1665             : GEN
    1666       62110 : ZM_hnflll(GEN A, GEN *ptB, int remove)
    1667             : {
    1668       62110 :   pari_sp av = avma;
    1669             :   long n, k, kmax;
    1670             :   GEN B, lambda, D;
    1671             : 
    1672       62110 :   n = lg(A);
    1673       62110 :   A = reverse_rows(ZM_copy(A)); /* ZM_copy for in place findi_normalize() */
    1674       62110 :   B = ptB? matid(n-1): NULL;
    1675       62110 :   D = const_vec(n, gen_1) + 1;
    1676       62110 :   lambda = zeromatcopy(n-1,n-1);
    1677       62110 :   k = kmax = 2;
    1678     5362461 :   while (k < n)
    1679             :   {
    1680             :     long row0, row1;
    1681             :     int do_swap;
    1682     5238241 :     reduce2(A,B,k,k-1,&row0,&row1,lambda,D);
    1683     5238241 :     if      (row0) do_swap = (!row1 || row0 <= row1);
    1684     1775863 :     else if (row1) do_swap = 0;
    1685     1497620 :     else do_swap = must_swap(k,lambda,D);
    1686     5238241 :     if (do_swap)
    1687             :     {
    1688     3006107 :       hnfswap(A,B,k,lambda,D);
    1689     3006107 :       if (k > 2) k--;
    1690             :     }
    1691             :     else
    1692             :     {
    1693             :       long i;
    1694    20277677 :       for (i=k-2; i; i--)
    1695             :       {
    1696             :         long row0, row1;
    1697    18045543 :         reduce2(A,B,k,i,&row0,&row1,lambda,D);
    1698    18045543 :         if (gc_needed(av,3))
    1699             :         {
    1700          14 :           GEN b = D-1;
    1701          14 :           if (DEBUGMEM>1) pari_warn(warnmem,"hnflll (reducing), kmax = %ld",kmax);
    1702          14 :           gerepileall(av, B? 4: 3, &A, &lambda, &b, &B);
    1703          14 :           D = b+1;
    1704             :         }
    1705             :       }
    1706     2232134 :       if (++k > kmax) kmax = k;
    1707             :     }
    1708     5238241 :     if (gc_needed(av,3))
    1709             :     {
    1710         140 :       GEN b = D-1;
    1711         140 :       if (DEBUGMEM>1) pari_warn(warnmem,"hnflll, kmax = %ld / %ld",kmax,n-1);
    1712         140 :       gerepileall(av, B? 4: 3, &A, &lambda, &b, &B);
    1713         140 :       D = b+1;
    1714             :     }
    1715             :   }
    1716             :   /* handle trivial case: return negative diag coefficient otherwise */
    1717       62110 :   if (n == 2) (void)findi_normalize(gel(A,1), B,1,lambda);
    1718       62110 :   A = reverse_rows(A);
    1719       62110 :   if (remove)
    1720             :   {
    1721             :     long i;
    1722       82164 :     for (i = 1; i < n; i++)
    1723       74450 :       if (!ZV_equal0(gel(A,i))) break;
    1724       12873 :     remove_0cols(i-1, &A, &B, remove);
    1725             :   }
    1726       62110 :   gerepileall(av, B? 2: 1, &A, &B);
    1727       62110 :   if (B) *ptB = B;
    1728       62110 :   return A;
    1729             : }
    1730             : 
    1731             : GEN
    1732           7 : hnflll(GEN x)
    1733             : {
    1734           7 :   GEN z = cgetg(3, t_VEC);
    1735           7 :   gel(z,1) = ZM_hnflll(x, &gel(z,2), 1);
    1736           7 :   return z;
    1737             : }
    1738             : 
    1739             : /* Variation on HNFLLL: Extended GCD */
    1740             : 
    1741             : static void
    1742      469886 : reduce1(GEN A, GEN B, long k, long j, GEN lambda, GEN D)
    1743             : {
    1744             :   GEN q;
    1745             :   long i;
    1746             : 
    1747      469886 :   if (signe(gel(A,j)))
    1748       92842 :     q = diviiround(gel(A,k),gel(A,j));
    1749      377044 :   else if (abscmpii(shifti(gcoeff(lambda,j,k), 1), gel(D,j)) > 0)
    1750       13871 :     q = diviiround(gcoeff(lambda,j,k), gel(D,j));
    1751             :   else
    1752      363173 :     return;
    1753             : 
    1754      106713 :   if (signe(q))
    1755             :   {
    1756       89233 :     GEN Lk = gel(lambda,k), Lj = gel(lambda,j);
    1757       89233 :     togglesign_safe(&q);
    1758       89233 :     gel(A,k) = addmulii(gel(A,k), q, gel(A,j));
    1759       89233 :     ZC_lincomb1_inplace(gel(B,k),gel(B,j),q);
    1760       89233 :     gel(Lk,j) = addmulii(gel(Lk,j), q, gel(D,j));
    1761      122013 :     for (i=1; i<j; i++)
    1762       32780 :       if (signe(gel(Lj,i))) gel(Lk,i) = addmulii(gel(Lk,i), q, gel(Lj,i));
    1763             :   }
    1764             : }
    1765             : 
    1766             : static GEN
    1767       87595 : ZV_gcdext_i(GEN A)
    1768             : {
    1769       87595 :   long k, n = lg(A);
    1770             :   GEN B, lambda, D;
    1771             : 
    1772       87595 :   if (n == 1) retmkvec2(gen_1, cgetg(1,t_MAT));
    1773       87595 :   A = leafcopy(A);
    1774       87595 :   B = matid(n-1);
    1775       87595 :   lambda = zeromatcopy(n-1,n-1);
    1776       87595 :   D = const_vec(n, gen_1) + 1;
    1777       87595 :   k = 2;
    1778      457363 :   while (k < n)
    1779             :   {
    1780             :     int do_swap;
    1781             : 
    1782      282173 :     reduce1(A,B,k,k-1,lambda,D);
    1783      282173 :     if    (signe(gel(A,k-1))) do_swap = 1;
    1784      189331 :     else if (signe(gel(A,k))) do_swap = 0;
    1785       71390 :     else do_swap = must_swap(k,lambda,D);
    1786      282173 :     if (do_swap)
    1787             :     {
    1788      100165 :       hnfswap(A,B,k,lambda,D);
    1789      100165 :       if (k > 2) k--;
    1790             :     }
    1791             :     else
    1792             :     {
    1793             :       long i;
    1794      182008 :       for (i=k-2; i; i--) reduce1(A,B,k,i,lambda,D);
    1795      182008 :       k++;
    1796             :     }
    1797             :   }
    1798       87595 :   if (signe(gel(A,n-1)) < 0)
    1799             :   {
    1800        7705 :     gel(A,n-1) = negi(gel(A,n-1));
    1801        7705 :     ZV_togglesign(gel(B,n-1));
    1802             :   }
    1803       87595 :   return mkvec2(gel(A,n-1), B);
    1804             : }
    1805             : GEN
    1806       87581 : ZV_extgcd(GEN A)
    1807             : {
    1808       87581 :   pari_sp av = avma;
    1809       87581 :   return gerepilecopy(av, ZV_gcdext_i(A));
    1810             : }
    1811             : /* as ZV_extgcd, transforming the gcd into a t_MAT, for mathnf0 */
    1812             : static GEN
    1813          21 : ZV_hnfgcdext(GEN A)
    1814             : {
    1815          21 :   pari_sp av = avma;
    1816             :   GEN z;
    1817          21 :   if (lg(A) == 1) retmkvec2(cgetg(1,t_MAT),cgetg(1,t_MAT));
    1818          14 :   z = ZV_gcdext_i(A);
    1819          14 :   gel(z,1) = mkmat(mkcol(gel(z,1)));
    1820          14 :   return gerepilecopy(av, z);
    1821             : }
    1822             : 
    1823             : /* HNF with permutation. */
    1824             : GEN
    1825         301 : ZM_hnfperm(GEN A, GEN *ptU, GEN *ptperm)
    1826             : {
    1827             :   GEN U, c, l, perm, d, p, q, b;
    1828         301 :   pari_sp av = avma, av1;
    1829             :   long r, t, i, j, j1, k, m, n;
    1830             : 
    1831         301 :   n = lg(A)-1;
    1832         301 :   if (!n)
    1833             :   {
    1834           7 :     if (ptU) *ptU = cgetg(1,t_MAT);
    1835           7 :     if (ptperm) *ptperm = cgetg(1,t_VEC);
    1836           7 :     return cgetg(1, t_MAT);
    1837             :   }
    1838         294 :   m = nbrows(A);
    1839         294 :   c = zero_zv(m);
    1840         294 :   l = zero_zv(n);
    1841         294 :   perm = cgetg(m+1, t_VECSMALL);
    1842         294 :   av1 = avma;
    1843         294 :   A = RgM_shallowcopy(A);
    1844         294 :   U = ptU? matid(n): NULL;
    1845             :   /* U base change matrix : A0*U = A all along */
    1846        1470 :   for (r=0, k=1; k <= n; k++)
    1847             :   {
    1848        3647 :     for (j=1; j<k; j++)
    1849             :     {
    1850        2471 :       if (!l[j]) continue;
    1851        2338 :       t=l[j]; b=gcoeff(A,t,k);
    1852        2338 :       if (!signe(b)) continue;
    1853             : 
    1854         553 :       ZC_elem(b,gcoeff(A,t,j), A,U,k,j);
    1855         553 :       d = gcoeff(A,t,j);
    1856         553 :       if (signe(d) < 0)
    1857             :       {
    1858           0 :         ZV_neg_inplace(gel(A,j));
    1859           0 :         if (U) ZV_togglesign(gel(U,j));
    1860           0 :         d = gcoeff(A,t,j);
    1861             :       }
    1862        1428 :       for (j1=1; j1<j; j1++)
    1863             :       {
    1864         875 :         if (!l[j1]) continue;
    1865         861 :         q = truedivii(gcoeff(A,t,j1),d);
    1866         861 :         if (!signe(q)) continue;
    1867             : 
    1868         329 :         togglesign(q);
    1869         329 :         ZC_lincomb1_inplace(gel(A,j1), gel(A,j), q);
    1870         329 :         if (U) ZC_lincomb1_inplace(gel(U,j1), gel(U,j), q);
    1871             :       }
    1872             :     }
    1873        1176 :     t = m; while (t && (c[t] || !signe(gcoeff(A,t,k)))) t--;
    1874        1176 :     if (t)
    1875             :     {
    1876        1078 :       p = gcoeff(A,t,k);
    1877       19313 :       for (i=t-1; i; i--)
    1878             :       {
    1879       18235 :         q = gcoeff(A,i,k);
    1880       18235 :         if (signe(q) && abscmpii(p,q) > 0) { p = q; t = i; }
    1881             :       }
    1882        1078 :       perm[++r] = l[k] = t; c[t] = k;
    1883        1078 :       if (signe(p) < 0)
    1884             :       {
    1885         174 :         ZV_neg_inplace(gel(A,k));
    1886         174 :         if (U) ZV_togglesign(gel(U,k));
    1887         174 :         p = gcoeff(A,t,k);
    1888             :       }
    1889             :       /* p > 0 */
    1890        3206 :       for (j=1; j<k; j++)
    1891             :       {
    1892        2128 :         if (!l[j]) continue;
    1893        2107 :         q = truedivii(gcoeff(A,t,j),p);
    1894        2107 :         if (!signe(q)) continue;
    1895             : 
    1896         231 :         togglesign(q);
    1897         231 :         ZC_lincomb1_inplace(gel(A,j), gel(A,k), q);
    1898         231 :         if (U) ZC_lincomb1_inplace(gel(U,j), gel(U,k), q);
    1899             :       }
    1900             :     }
    1901        1176 :     if (gc_needed(av1,1))
    1902             :     {
    1903           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"hnfperm, k=%ld",k);
    1904           0 :       gerepileall(av1, U? 2: 1, &A, &U);
    1905             :     }
    1906             :   }
    1907         294 :   if (r < m)
    1908             :   {
    1909        4732 :     for (i=1,k=r; i<=m; i++)
    1910        4494 :       if (!c[i]) perm[++k] = i;
    1911             :   }
    1912             : 
    1913             : /* We have A0*U=A, U in Gl(n,Z)
    1914             :  * basis for Im(A):  columns of A s.t l[j]>0 (r   cols)
    1915             :  * basis for Ker(A): columns of U s.t l[j]=0 (n-r cols) */
    1916         294 :   p = cgetg(r+1,t_MAT);
    1917         294 :   for (i=1; i<=m/2; i++) lswap(perm[i], perm[m+1-i]);
    1918         294 :   if (U)
    1919             :   {
    1920          70 :     GEN u = cgetg(n+1,t_MAT);
    1921         322 :     for (t=1,k=r,j=1; j<=n; j++)
    1922         252 :       if (l[j])
    1923             :       {
    1924         154 :         u[k + n-r] = U[j];
    1925         154 :         gel(p,k--) = vecpermute(gel(A,j), perm);
    1926             :       }
    1927             :       else
    1928          98 :         u[t++] = U[j];
    1929          70 :     *ptU = u;
    1930          70 :     if (ptperm) *ptperm = perm;
    1931          70 :     gerepileall(av, ptperm? 3: 2, &p, ptU, ptperm);
    1932             :   }
    1933             :   else
    1934             :   {
    1935        1148 :     for (k=r,j=1; j<=n; j++)
    1936         924 :       if (l[j]) gel(p,k--) = vecpermute(gel(A,j), perm);
    1937         224 :     if (ptperm) *ptperm = perm;
    1938         224 :     gerepileall(av, ptperm? 2: 1, &p, ptperm);
    1939             :   }
    1940         294 :   return p;
    1941             : }
    1942             : 
    1943             : GEN
    1944         217 : ZM_hnf_knapsack(GEN x)
    1945             : {
    1946         217 :   GEN t, perm, H = ZM_hnfperm(x,NULL,&perm);
    1947         217 :   long i,j, l = lg(H), h = lgcols(H);
    1948        3136 :   for (i=1; i<h; i++)
    1949             :   {
    1950        3003 :     int fl = 0;
    1951       16681 :     for (j=1; j<l; j++)
    1952             :     {
    1953       13762 :       t = gcoeff(H,i,j);
    1954       13762 :       if (signe(t))
    1955             :       {
    1956        3038 :         if (!is_pm1(t) || fl) return NULL;
    1957        2954 :         fl = 1;
    1958             :       }
    1959             :     }
    1960             :   }
    1961         133 :   return rowpermute(H, perm_inv(perm));
    1962             : }
    1963             : 
    1964             : GEN
    1965          14 : hnfperm(GEN A)
    1966             : {
    1967          14 :   GEN y = cgetg(4, t_VEC);
    1968          14 :   gel(y,1) = ZM_hnfperm(A, &gel(y,2), &gel(y,3));
    1969          14 :   return y;
    1970             : }
    1971             : 
    1972             : /* Hermite Normal Form, with base change matrix if ptB != NULL.
    1973             :  * If 'remove' = 1, remove 0 columns (do NOT update *ptB accordingly)
    1974             :  * If 'remove' = 2, remove 0 columns and update *ptB accordingly */
    1975             : GEN
    1976      135488 : ZM_hnfall_i(GEN A, GEN *ptB, long remove)
    1977             : {
    1978             :   pari_sp av;
    1979             :   long m, n, r, i, j, k, li;
    1980             :   GEN B, c, h, a;
    1981             : 
    1982      135488 :   RgM_dimensions(A, &m,&n);
    1983      135488 :   if (!n)
    1984             :   {
    1985           7 :     if (ptB) *ptB = cgetg(1,t_MAT);
    1986           7 :     return cgetg(1,t_MAT);
    1987             :   }
    1988      135481 :   c = zero_zv(m);
    1989      135481 :   h = const_vecsmall(n, m);
    1990      135481 :   av = avma;
    1991      135481 :   A = RgM_shallowcopy(A);
    1992      135481 :   B = ptB? matid(n): NULL;
    1993      135481 :   r = n+1;
    1994      448362 :   for (li=m; li; li--)
    1995             :   {
    1996      941357 :     for (j=1; j<r; j++)
    1997             :     {
    1998     1188429 :       for (i=h[j]; i>li; i--)
    1999             :       {
    2000      247254 :         a = gcoeff(A,i,j);
    2001      247254 :         k = c[i];
    2002             :         /* zero a = Aij  using  Aik */
    2003      247254 :         if (signe(a)) ZC_elem(a,gcoeff(A,i,k), A,B,j,k);
    2004      247254 :         ZM_reduce(A,B, i,k); /* ensure reduced entries even if a = 0 */
    2005      247254 :         if (gc_needed(av,1))
    2006             :         {
    2007           0 :           if (DEBUGMEM>1)
    2008           0 :             pari_warn(warnmem,"ZM_hnfall[1], li = %ld, j = %ld", li, j);
    2009           0 :           gerepileall(av, B? 2: 1, &A, &B);
    2010             :         }
    2011             :       }
    2012      941175 :       if (signe( gcoeff(A,li,j) )) break;
    2013      628476 :       h[j] = li-1;
    2014             :     }
    2015      312881 :     if (j == r) continue;
    2016      312699 :     r--;
    2017      312699 :     if (j < r) /* A[j] != 0 */
    2018             :     {
    2019      230901 :       swap(gel(A,j), gel(A,r));
    2020      230901 :       if (B) swap(gel(B,j), gel(B,r));
    2021      230901 :       h[j] = h[r]; h[r] = li; c[li] = r;
    2022             :     }
    2023      312699 :     if (signe(gcoeff(A,li,r)) < 0)
    2024             :     {
    2025       76051 :       ZV_neg_inplace(gel(A,r));
    2026       76051 :       if (B) ZV_togglesign(gel(B,r));
    2027             :     }
    2028      312699 :     ZM_reduce(A,B, li,r);
    2029      312699 :     if (gc_needed(av,1))
    2030             :     {
    2031           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"ZM_hnfall[2], li = %ld", li);
    2032           0 :       gerepileall(av, B? 2: 1, &A, &B);
    2033             :     }
    2034             :   }
    2035             : 
    2036      135481 :   if (DEBUGLEVEL>5) err_printf("\nhnfall, final phase: ");
    2037      135481 :   r--; /* first r cols are in the image the n-r (independent) last ones */
    2038      760566 :   for (j=1; j<=r; j++)
    2039     1477810 :     for (i=h[j]; i; i--)
    2040             :     {
    2041      852725 :       a = gcoeff(A,i,j);
    2042      852725 :       k = c[i];
    2043      852725 :       if (signe(a)) ZC_elem(a,gcoeff(A,i,k), A,B, j,k);
    2044      852725 :       ZM_reduce(A,B, i,k); /* ensure reduced entries, even if a = 0 */
    2045      852725 :       if (gc_needed(av,1))
    2046             :       {
    2047           0 :         if (DEBUGMEM>1) pari_warn(warnmem,"ZM_hnfall[3], j = %ld", j);
    2048           0 :         gerepileall(av, B? 2: 1, &A, &B);
    2049             :       }
    2050             :     }
    2051      135481 :   if (DEBUGLEVEL>5) err_printf("\n");
    2052      135481 :   if (remove) remove_0cols(r, &A, &B, remove);
    2053      135481 :   if (ptB) *ptB = B;
    2054      135481 :   return A;
    2055             : }
    2056             : GEN
    2057        4385 : ZM_hnfall(GEN A, GEN *ptB, long remove)
    2058             : {
    2059        4385 :   pari_sp av = avma;
    2060        4385 :   A = ZM_hnfall_i(A, ptB, remove);
    2061        4385 :   gerepileall(av, ptB? 2: 1, &A, ptB);
    2062        4385 :   return A;
    2063             : }
    2064             : 
    2065             : GEN
    2066          14 : hnfall(GEN x)
    2067             : {
    2068          14 :   GEN z = cgetg(3, t_VEC);
    2069          14 :   gel(z,1) = ZM_hnfall(x, (GEN*)(z+2), 1);
    2070          14 :   return z;
    2071             : }
    2072             : GEN
    2073          14 : hnf(GEN x) { return mathnf0(x,0); }
    2074             : 
    2075             : /* C = A^(-1)t where A and C are integral, A is upper triangular, t t_INT */
    2076             : GEN
    2077       29897 : hnf_invscale(GEN A, GEN t)
    2078             : {
    2079       29897 :   long n = lg(A)-1, i,j,k;
    2080       29897 :   GEN m, c = cgetg(n+1,t_MAT);
    2081             : 
    2082       29897 :   if (!n) return c;
    2083      148022 :   for (k=1; k<=n; k++)
    2084             :   { /* cf hnf_divscale with B = id, thus b = e_k */
    2085      118125 :     GEN u = cgetg(n+1, t_COL);
    2086      118125 :     pari_sp av = avma;
    2087      118125 :     gel(c,k) = u;
    2088      118125 :     gel(u,n) = k == n? gerepileuptoint(av, diviiexact(t, gcoeff(A,n,n))): gen_0;
    2089      717577 :     for (i=n-1; i>0; i--)
    2090             :     {
    2091      599452 :       av = avma; m = i == k? t: gen_0;
    2092      599452 :       for (j=i+1; j<=n; j++) m = subii(m, mulii(gcoeff(A,i,j),gel(u,j)));
    2093      599452 :       gel(u,i) = gerepileuptoint(av, diviiexact(m, gcoeff(A,i,i)));
    2094             :     }
    2095             :   }
    2096       29897 :   return c;
    2097             : }
    2098             : 
    2099             : /* C = A^(-1)(tB) where A, B, C are integral, A is upper triangular, t t_INT */
    2100             : GEN
    2101      138561 : hnf_divscale(GEN A, GEN B, GEN t)
    2102             : {
    2103      138561 :   long n = lg(A)-1, i,j,k;
    2104      138561 :   GEN m, c = cgetg(n+1,t_MAT);
    2105             : 
    2106      138561 :   if (!n) return c;
    2107      558028 :   for (k=1; k<=n; k++)
    2108             :   {
    2109      419467 :     GEN u = cgetg(n+1, t_COL), b = gel(B,k);
    2110      419467 :     pari_sp av = avma;
    2111      419467 :     gel(c,k) = u; m = mulii(gel(b,n),t);
    2112      419467 :     gel(u,n) = gerepileuptoint(av, diviiexact(m, gcoeff(A,n,n)));
    2113     1947953 :     for (i=n-1; i>0; i--)
    2114             :     {
    2115     1528486 :       av = avma; m = mulii(gel(b,i),t);
    2116     1528486 :       for (j=i+1; j<=n; j++) m = subii(m, mulii(gcoeff(A,i,j),gel(u,j)));
    2117     1528486 :       gel(u,i) = gerepileuptoint(av, diviiexact(m, gcoeff(A,i,i)));
    2118             :     }
    2119             :   }
    2120      138561 :   return c;
    2121             : }
    2122             : 
    2123             : /* A, B integral upper HNF. A^(-1) B integral ? */
    2124             : int
    2125          91 : hnfdivide(GEN A, GEN B)
    2126             : {
    2127          91 :   pari_sp av = avma;
    2128          91 :   long n = lg(A)-1, i,j,k;
    2129             :   GEN u, b, m, r;
    2130             : 
    2131          91 :   if (!n) return 1;
    2132          91 :   if (lg(B)-1 != n) pari_err_DIM("hnfdivide");
    2133          91 :   u = cgetg(n+1, t_COL);
    2134         413 :   for (k=1; k<=n; k++)
    2135             :   {
    2136         322 :     b = gel(B,k);
    2137         322 :     m = gel(b,k);
    2138         322 :     gel(u,k) = dvmdii(m, gcoeff(A,k,k), &r);
    2139         322 :     if (r != gen_0) return gc_long(av, 0);
    2140         798 :     for (i=k-1; i>0; i--)
    2141             :     {
    2142         476 :       m = gel(b,i);
    2143         476 :       for (j=i+1; j<=k; j++) m = subii(m, mulii(gcoeff(A,i,j),gel(u,j)));
    2144         476 :       m = dvmdii(m, gcoeff(A,i,i), &r);
    2145         476 :       if (r != gen_0) return gc_long(av, 0);
    2146         476 :       gel(u,i) = m;
    2147             :     }
    2148             :   }
    2149          91 :   return gc_long(av, 1);
    2150             : }
    2151             : 
    2152             : /* A upper HNF, b integral vector. Return A^(-1) b if integral,
    2153             :  * NULL otherwise. Assume #A[,1] = #b. */
    2154             : GEN
    2155      115445 : hnf_invimage(GEN A, GEN b)
    2156             : {
    2157      115445 :   pari_sp av = avma;
    2158      115445 :   long n = lg(A)-1, m, i, k;
    2159             :   GEN u, r;
    2160             : 
    2161      115445 :   if (!n) return lg(b)==1? cgetg(1,t_COL):NULL;
    2162      115403 :   m = nbrows(A); /* m >= n */
    2163      115403 :   u = cgetg(n+1, t_COL);
    2164      324005 :   for (i = n, k = m; k > 0; k--)
    2165             :   {
    2166      324005 :     pari_sp av2 = avma;
    2167             :     long j;
    2168      324005 :     GEN t = gel(b,k), Aki = gcoeff(A,k,i);
    2169      324005 :     if (typ(t) != t_INT) pari_err_TYPE("hnf_invimage",t);
    2170      324005 :     for (j=i+1; j<=n; j++) t = subii(t, mulii(gcoeff(A,k,j),gel(u,j)));
    2171      324005 :     if (!signe(Aki))
    2172             :     {
    2173           0 :       if (signe(t)) return gc_NULL(av);
    2174           0 :       set_avma(av2); gel(u,i) = gen_0; continue;
    2175             :     }
    2176      324005 :     t = dvmdii(t, Aki, &r);
    2177      324005 :     if (r != gen_0) return gc_NULL(av);
    2178      269866 :     gel(u,i) = gerepileuptoint(av2, t);
    2179      269866 :     if (--i == 0) break;
    2180             :   }
    2181             :   /* If there is a solution, it must be u. Check remaining equations */
    2182      122528 :   for (; k > 0; k--)
    2183             :   {
    2184       61264 :     pari_sp av2 = avma;
    2185             :     long j;
    2186       61264 :     GEN t = gel(b,k);
    2187       61264 :     if (typ(t) != t_INT) pari_err_TYPE("hnf_invimage",t);
    2188       61264 :     for (j=1; j<=n; j++) t = subii(t, mulii(gcoeff(A,k,j),gel(u,j)));
    2189       61264 :     if (signe(t)) return gc_NULL(av);
    2190       61264 :     set_avma(av2);
    2191             :   }
    2192       61264 :   return u;
    2193             : }
    2194             : 
    2195             : /* A upper HNF, B integral matrix or column. Return A^(-1) B if integral,
    2196             :  * NULL otherwise */
    2197             : GEN
    2198      102074 : hnf_solve(GEN A, GEN B)
    2199             : {
    2200             :   pari_sp av;
    2201             :   long i, l;
    2202             :   GEN C;
    2203             : 
    2204      102074 :   if (typ(B) == t_COL) return hnf_invimage(A, B);
    2205       73122 :   av = avma; C = cgetg_copy(B, &l);
    2206      114037 :   for (i = 1; i < l; i++)
    2207             :   {
    2208       75824 :     GEN c = hnf_invimage(A, gel(B,i));
    2209       75824 :     if (!c) return gc_NULL(av);
    2210       40915 :     gel(C,i) = c;
    2211             :   }
    2212       38213 :   return C;
    2213             : }
    2214             : 
    2215             : /***************************************************************/
    2216             : /**                                                           **/
    2217             : /**               SMITH NORMAL FORM REDUCTION                 **/
    2218             : /**                                                           **/
    2219             : /***************************************************************/
    2220             : 
    2221             : static GEN
    2222           0 : trivsmith(long all)
    2223             : {
    2224             :   GEN z;
    2225           0 :   if (!all) return cgetg(1,t_VEC);
    2226           0 :   z=cgetg(4,t_VEC);
    2227           0 :   gel(z,1) = cgetg(1,t_MAT);
    2228           0 :   gel(z,2) = cgetg(1,t_MAT);
    2229           0 :   gel(z,3) = cgetg(1,t_MAT); return z;
    2230             : }
    2231             : 
    2232             : static void
    2233          38 : snf_pile1(pari_sp av, GEN *x, GEN *U)
    2234             : {
    2235             :   GEN *gptr[2];
    2236          38 :   int c = 1; gptr[0]=x;
    2237          38 :   if (*U) gptr[c++] = U;
    2238          38 :   gerepilemany(av,gptr,c);
    2239          38 : }
    2240             : static void
    2241      203887 : snf_pile(pari_sp av, GEN *x, GEN *U, GEN *V)
    2242             : {
    2243             :   GEN *gptr[3];
    2244      203887 :   int c = 1; gptr[0]=x;
    2245      203887 :   if (*U) gptr[c++] = U;
    2246      203887 :   if (*V) gptr[c++] = V;
    2247      203887 :   gerepilemany(av,gptr,c);
    2248      203887 : }
    2249             : 
    2250             : static GEN
    2251      183569 : bezout_step(GEN *pa, GEN *pb, GEN *pu, GEN *pv)
    2252             : {
    2253      183569 :   GEN a = *pa, b = *pb, d;
    2254      183569 :   if (absequalii(a,b))
    2255             :   {
    2256       82791 :     long sa = signe(a), sb = signe(b);
    2257       82791 :     *pv = gen_0;
    2258       82791 :     if (sb == sa) {
    2259       80273 :       *pa = *pb = gen_1;
    2260       80273 :       if (sa > 0) {*pu=gen_1; return a;} else {*pu=gen_m1; return absi(a);}
    2261             :     }
    2262        2518 :     if (sa > 0) { *pa = *pu = gen_1; *pb = gen_m1; return a; }
    2263          35 :     *pa = *pu = gen_m1; *pb = gen_1; return b;
    2264             :   }
    2265      100778 :   d = bezout(a,b, pu,pv);
    2266      100778 :   *pa = diviiexact(a, d);
    2267      100778 :   *pb = diviiexact(b, d); return d;
    2268             : }
    2269             : 
    2270             : static int
    2271      188447 : negcmpii(void *E, GEN x, GEN y) { (void)E; return -cmpii(x,y); }
    2272             : 
    2273             : /* x square of maximal rank; does b = x[i,i] divide all entries in
    2274             :  * x[1..i-1, 1..i-1] ? If so, return 0; else the index of a problematic row */
    2275             : static long
    2276      232277 : ZM_snf_no_divide(GEN x, long i)
    2277             : {
    2278      232277 :   GEN b = gcoeff(x,i,i);
    2279             :   long j, k;
    2280             : 
    2281      232277 :   if (is_pm1(b)) return 0;
    2282      400083 :   for (k = 1; k < i; k++)
    2283     1877259 :     for (j = 1; j < i; j++)
    2284     1599166 :       if (!dvdii(gcoeff(x,k,j),b)) return k;
    2285       84854 :   return 0;
    2286             : }
    2287             : 
    2288             : static void
    2289      284972 : ZM_redpart(GEN x, GEN p, long I)
    2290             : {
    2291      284972 :   long l = lgefint(p), i, j;
    2292     1864818 :   for (i = 1; i <= I; i++)
    2293    15745674 :     for (j = 1; j <= I; j++)
    2294             :     {
    2295    14165828 :       GEN c = gcoeff(x,i,j);
    2296    14165828 :       if (lgefint(c) > l) gcoeff(x,i,j) = remii(c, p);
    2297             :     }
    2298      284972 : }
    2299             : static void
    2300      175341 : ZMrow_divexact_inplace(GEN M, long i, GEN c)
    2301             : {
    2302      175341 :   long j, l = lg(M);
    2303      175341 :   for (j = 1; j < l; j++) gcoeff(M,i,j) = diviiexact(gcoeff(M,i,j), c);
    2304      175341 : }
    2305             : 
    2306             : /* Return the SNF D of matrix X. If ptU/ptV non-NULL set them to U/V
    2307             :  * to that D = UXV */
    2308             : GEN
    2309      107764 : ZM_snfall_i(GEN x, GEN *ptU, GEN *ptV, long flag)
    2310             : {
    2311      107764 :   pari_sp av0 = avma, av;
    2312      107764 :   const long return_vec = flag & 1;
    2313             :   long i, j, k, m0, m, n0, n;
    2314      107764 :   GEN u, v, U, V, V0, mdet, A = NULL, perm = NULL;
    2315             : 
    2316      107764 :   n0 = n = lg(x)-1;
    2317      107764 :   if (!n) {
    2318        4607 :     if (ptU) *ptU = cgetg(1,t_MAT);
    2319        4607 :     if (ptV) *ptV = cgetg(1,t_MAT);
    2320        4607 :     return cgetg(1, return_vec? t_VEC: t_MAT);
    2321             :   }
    2322      103157 :   m0 = m = nbrows(x);
    2323             : 
    2324      103157 :   U = V = V0 = NULL; /* U = TRANSPOSE of row transform matrix [act on columns]*/
    2325      103157 :   if (m == n && ZM_ishnf(x))
    2326             :   {
    2327       93588 :     mdet = ZM_det_triangular(x); /* != 0 */
    2328       93588 :     if (ptV) *ptV = matid(n);
    2329             :   }
    2330             :   else
    2331             :   {
    2332        9569 :     mdet = ZM_detmult(x);
    2333        9569 :     if (!signe(mdet))
    2334          21 :       x = ZM_hnfperm(x, ptV, ptU? &perm: NULL);
    2335             :     else
    2336             :     { /* m <= n */
    2337        9548 :       if (!ptV)
    2338         287 :         x = ZM_hnfmod(x,mdet);
    2339        9261 :       else if (m == n)
    2340             :       {
    2341        9240 :         GEN H = ZM_hnfmod(x,mdet);
    2342        9240 :         *ptV = ZM_gauss(x,H);
    2343        9240 :         x = H;
    2344             :       }
    2345             :       else
    2346          21 :         x = ZM_hnfperm(x, ptV, ptU? &perm: NULL);
    2347        9548 :       mdet = ZM_det_triangular(x); /* > 0 */
    2348             :     }
    2349        9569 :     n = lg(x)-1; /* n independent columns */
    2350        9569 :     if (ptV)
    2351             :     {
    2352        9275 :       V = *ptV;
    2353        9275 :       if (n != n0)
    2354             :       {
    2355          21 :         V0 = vecslice(V, 1, n0 - n); /* kernel */
    2356          21 :         V  = vecslice(V, n0-n+1, n0);
    2357             :       }
    2358             :     }
    2359        9569 :     if (!signe(mdet))
    2360             :     {
    2361          21 :       if (n)
    2362             :       {
    2363          21 :         x = ZM_snfall_i(shallowtrans(x), ptV, ptU, return_vec); /* swap V,U */
    2364          21 :         if (!return_vec && n != m) x = shallowtrans(x);
    2365          21 :         if (ptV) V = ZM_mul(V, shallowtrans(*ptV));
    2366          21 :         if (ptU) U = *ptU; /* TRANSPOSE */
    2367             :       }
    2368             :       else /* 0 matrix */
    2369             :       {
    2370           0 :         x = cgetg(1,t_MAT);
    2371           0 :         if (ptV) V = cgetg(1, t_MAT);
    2372           0 :         if (ptU) U = matid(m);
    2373             :       }
    2374          21 :       goto THEEND;
    2375             :     }
    2376             :   }
    2377      103136 :   if (ptV || ptU) U = matid(n); /* we will compute V in terms of U */
    2378      103136 :   if (DEBUGLEVEL>7) err_printf("starting SNF loop");
    2379             : 
    2380             :   /* square, maximal rank n */
    2381      103136 :   A = x; x = shallowcopy(x); av = avma;
    2382      596554 :   for (i = n; i > 1; i--)
    2383             :   {
    2384      195141 :     if (DEBUGLEVEL>7) err_printf("\ni = %ld: ",i);
    2385             :     for(;;)
    2386      145550 :     {
    2387      340691 :       int c = 0;
    2388             :       GEN a, b;
    2389     1567803 :       for (j = i-1; j >= 1; j--)
    2390             :       {
    2391     1227112 :         b = gcoeff(x,i,j); if (!signe(b)) continue;
    2392       85371 :         a = gcoeff(x,i,i);
    2393       85371 :         ZC_elem(b, a, x,NULL, j,i);
    2394       85371 :         if (gc_needed(av,1))
    2395             :         {
    2396           1 :           if (DEBUGMEM>1) pari_warn(warnmem,"[1]: ZM_snfall i = %ld", i);
    2397           1 :           snf_pile1(av, &x,&U);
    2398             :         }
    2399             :       }
    2400      340691 :       if (DEBUGLEVEL>7) err_printf("; ");
    2401     1567803 :       for (j=i-1; j>=1; j--)
    2402             :       {
    2403             :         GEN d;
    2404     1227112 :         b = gcoeff(x,j,i); if (!signe(b)) continue;
    2405      183569 :         a = gcoeff(x,i,i);
    2406      183569 :         d = bezout_step(&a, &b, &u, &v);
    2407     2785860 :         for (k = 1; k < i; k++)
    2408             :         {
    2409     2602291 :           GEN t = addii(mulii(u,gcoeff(x,i,k)),mulii(v,gcoeff(x,j,k)));
    2410     5204582 :           gcoeff(x,j,k) = subii(mulii(a,gcoeff(x,j,k)),
    2411     2602291 :                                 mulii(b,gcoeff(x,i,k)));
    2412     2602291 :           gcoeff(x,i,k) = t;
    2413             :         }
    2414      183569 :         gcoeff(x,j,i) = gen_0;
    2415      183569 :         gcoeff(x,i,i) = d;
    2416      183569 :         if (U) update(u,v,a,b,(GEN*)(U+i),(GEN*)(U+j));
    2417      183569 :         if (gc_needed(av,1))
    2418             :         {
    2419          34 :           if (DEBUGMEM>1) pari_warn(warnmem,"[2]: ZM_snfall, i = %ld", i);
    2420          34 :           snf_pile1(av, &x,&U);
    2421             :         }
    2422      183569 :         c = 1;
    2423             :       }
    2424      340691 :       if (!c)
    2425             :       {
    2426      232277 :         k = ZM_snf_no_divide(x, i);
    2427      232277 :         if (!k) break;
    2428             : 
    2429             :         /* x[k,j] != 0 mod b */
    2430      233921 :         for (j = 1; j <= i; j++)
    2431      196785 :           gcoeff(x,i,j) = addii(gcoeff(x,i,j),gcoeff(x,k,j));
    2432       37136 :         if (U) gel(U,i) = gadd(gel(U,i),gel(U,k));
    2433             :       }
    2434      145550 :       ZM_redpart(x, mdet, i);
    2435      145550 :       if (U && (flag & 2)) ZM_redpart(U, mdet, n);
    2436      145550 :       if (gc_needed(av,1))
    2437             :       {
    2438           3 :         if (DEBUGMEM>1) pari_warn(warnmem,"[3]: ZM_snfall");
    2439           3 :         snf_pile1(av, &x,&U);
    2440             :       }
    2441             :     }
    2442             :   }
    2443      103136 :   if (DEBUGLEVEL>7) err_printf("\n");
    2444      401413 :   for (k = n; k; k--)
    2445             :   {
    2446      298277 :     GEN d = gcdii(gcoeff(x,k,k), mdet);
    2447      298277 :     gcoeff(x,k,k) = d;
    2448      298277 :     if (!is_pm1(d)) mdet = diviiexact(mdet,d);
    2449             :   }
    2450             : THEEND:
    2451      103157 :   if (U) U = shallowtrans(U);
    2452      103157 :   if (ptV && A)
    2453             :   { /* U A V = D => D^(-1) U A = V^(-1) */
    2454       99048 :     long l = lg(x);
    2455       99048 :     GEN W = ZM_mul(U, A);
    2456      274389 :     for (i = 1; i < l; i++)
    2457             :     {
    2458      239651 :       GEN c = gcoeff(x,i,i);
    2459      239651 :       if (is_pm1(c)) break; /* only 1 from now on */
    2460      175341 :       ZMrow_divexact_inplace(W, i, c);
    2461             :     }
    2462       99048 :     if (flag & 2)
    2463             :     {
    2464       91257 :       W = FpM_red(W, gcoeff(x,1,1));
    2465       91257 :       W = matinvmod(W, gcoeff(x,1,1));
    2466             :     }
    2467             :     else
    2468        7791 :       W = ZM_inv(W, NULL);
    2469       99048 :     V = V? ZM_mul(V, W): W;
    2470             :   }
    2471      103157 :   if (return_vec)
    2472             :   {
    2473       95184 :     long l = lg(x)-1;
    2474       95184 :     if (typ(x) == t_MAT) x = RgM_diagonal_shallow(x);
    2475       95184 :     if (m0 > l) x = shallowconcat(zerovec(m0-l), x);
    2476             :   }
    2477             : 
    2478      103157 :   if (V0)
    2479             :   { /* add kernel */
    2480          21 :     if (!return_vec) x = shallowconcat(zeromat(m,n0-n), x);
    2481          21 :     if (ptV) V = shallowconcat(V0, V);
    2482             :   }
    2483      103157 :   if (perm && U) U = vecpermute(U, perm_inv(perm));
    2484      103157 :   snf_pile(av0, &x,&U,&V);
    2485      103157 :   if (ptU) *ptU = U;
    2486      103157 :   if (ptV) *ptV = V;
    2487      103157 :   return x;
    2488             : }
    2489             : GEN
    2490       11991 : ZM_snfall(GEN x, GEN *U, GEN *V) { return ZM_snfall_i(x, U, V, 0); }
    2491             : GEN
    2492         427 : ZM_snf(GEN x) { return ZM_snfall_i(x, NULL,NULL, 1); }
    2493             : GEN
    2494         329 : smith(GEN x) { return ZM_snfall_i(x, NULL,NULL, 1); }
    2495             : GEN
    2496          28 : smithall(GEN x)
    2497             : {
    2498          28 :   GEN z = cgetg(4, t_VEC);
    2499          28 :   gel(z,3) = ZM_snfall_i(x, (GEN*)(z+1),(GEN*)(z+2), 0);
    2500          28 :   return z;
    2501             : }
    2502             : 
    2503             : void
    2504      192555 : ZM_snfclean(GEN d, GEN u, GEN v)
    2505             : {
    2506      192555 :   long i, c, l = lg(d);
    2507             : 
    2508      192555 :   if (typ(d) == t_VEC)
    2509      192555 :     for (c=1; c<l; c++) { GEN t = gel(d,c); if (is_pm1(t)) break; }
    2510             :   else
    2511             :   {
    2512           0 :     for (c=1; c<l; c++) { GEN t = gcoeff(d,c,c); if (is_pm1(t)) break; }
    2513           0 :     if (c < l) for (i = 1; i < c; i++) setlg(gel(d,i), c);
    2514             :   }
    2515      192555 :   setlg(d, c);
    2516      192555 :   if (u) for (i=1; i<l; i++) setlg(gel(u,i), c);
    2517      192555 :   if (v) setlg(v, c);
    2518      192555 : }
    2519             : 
    2520             : /* Assume z was computed by [g]smithall(). Remove the 1s on the diagonal */
    2521             : GEN
    2522         903 : smithclean(GEN z)
    2523             : {
    2524             :   long i, j, h, l, c, d;
    2525             :   GEN U, V, y, D, t;
    2526             : 
    2527         903 :   if (typ(z) != t_VEC) pari_err_TYPE("smithclean",z);
    2528         903 :   l = lg(z); if (l == 1) return cgetg(1,t_VEC);
    2529         896 :   U = gel(z,1);
    2530         896 :   if (l != 4 || typ(U) != t_MAT)
    2531             :   { /* assume z = vector of elementary divisors */
    2532        2373 :     for (c=1; c<l; c++)
    2533        1883 :       if (gequal1(gel(z,c))) break;
    2534         875 :     return gcopy_lg(z, c);
    2535             :   }
    2536          21 :   V = gel(z,2);
    2537          21 :   D = gel(z,3);
    2538          21 :   l = lg(D);
    2539          21 :   if (l == 1) return gcopy(z);
    2540          21 :   h = lgcols(D);
    2541          21 :   if (h > l)
    2542             :   { /* D = vconcat(zero matrix, diagonal matrix) */
    2543          21 :     for (c=1+h-l, d=1; c<h; c++,d++)
    2544          21 :       if (gequal1(gcoeff(D,c,d))) break;
    2545             :   }
    2546           7 :   else if (h < l)
    2547             :   { /* D = concat(zero matrix, diagonal matrix) */
    2548           7 :     for (c=1, d=1+l-h; d<l; c++,d++)
    2549           7 :       if (gequal1(gcoeff(D,c,d))) break;
    2550             :   }
    2551             :   else
    2552             :   { /* D diagonal */
    2553           0 :     for (c=1; c<l; c++)
    2554           0 :       if (gequal1(gcoeff(D,c,c))) break;
    2555           0 :     d = c;
    2556             :   }
    2557             :   /* U was (h-1)x(h-1), V was (l-1)x(l-1), D was (h-1)x(l-1) */
    2558          21 :   y = cgetg(4,t_VEC);
    2559             :   /* truncate U to (c-1) x (h-1) */
    2560          21 :   gel(y,1) = t = cgetg(h,t_MAT);
    2561          21 :   for (j=1; j<h; j++) gel(t,j) = gcopy_lg(gel(U,j), c);
    2562             :   /* truncate V to (l-1) x (d-1) */
    2563          21 :   gel(y,2) = gcopy_lg(V, d);
    2564          21 :   gel(y,3) = t = zeromatcopy(c-1, d-1);
    2565             :   /* truncate D to a (c-1) x (d-1) matrix */
    2566          21 :   if (d > 1)
    2567             :   {
    2568          14 :     if (h > l)
    2569             :     {
    2570          14 :       for (i=1+h-l, j=1; i<c; i++,j++)
    2571           7 :         gcoeff(t,i,j) = gcopy(gcoeff(D,i,j));
    2572             :     }
    2573           7 :     else if (h < l)
    2574             :     {
    2575           7 :       for (i=1, j=1+l-h; j<d; i++,j++)
    2576           0 :         gcoeff(t,i,j) = gcopy(gcoeff(D,i,j));
    2577             :     }
    2578             :     else
    2579             :     {
    2580           0 :       for (j=1; j<d; j++)
    2581           0 :         gcoeff(t,j,j) = gcopy(gcoeff(D,j,j));
    2582             :     }
    2583             :   }
    2584          21 :   return y;
    2585             : }
    2586             : 
    2587             : /* does b = x[i,i] divide all entries in x[1..i-1,1..i-1] ? If so, return 0;
    2588             :  * else return the index of a problematic row */
    2589             : static long
    2590         112 : gsnf_no_divide(GEN x, long i, long vx)
    2591             : {
    2592         112 :   GEN b = gcoeff(x,i,i);
    2593             :   long j, k;
    2594             : 
    2595         112 :   if (gequal0(b))
    2596             :   {
    2597          14 :     for (k = 1; k < i; k++)
    2598          14 :       for (j = 1; j < i; j++)
    2599          14 :         if (!gequal0(gcoeff(x,k,j))) return k;
    2600           0 :     return 0;
    2601             :   }
    2602             : 
    2603          98 :   if (!is_RgX(b,vx) || degpol(b)<=0) return 0;
    2604         133 :   for (k = 1; k < i; k++)
    2605         231 :     for (j = 1; j < i; j++)
    2606             :     {
    2607         154 :       GEN z = gcoeff(x,k,j), r;
    2608         154 :       if (!is_RgX(z,vx)) z = scalarpol(z, vx);
    2609         154 :       r = RgX_rem(z, b);
    2610         154 :       if (signe(r) && (! isinexactreal(r) ||
    2611           0 :              gexpo(r) > 16 + gexpo(b) - prec2nbits(gprecision(r)))
    2612           7 :          ) return k;
    2613             :     }
    2614          49 :   return 0;
    2615             : }
    2616             : 
    2617             : /* Hermite Normal Form, with base change matrix if ptB != NULL.
    2618             :  * If 'remove' = 1, remove 0 columns (do NOT update *ptB accordingly)
    2619             :  * If 'remove' = 2, remove 0 columns and update *ptB accordingly */
    2620             : GEN
    2621          35 : RgM_hnfall(GEN A, GEN *pB, long remove)
    2622             : {
    2623             :   pari_sp av;
    2624             :   long li, j, k, m, n, def, ldef;
    2625             :   GEN B;
    2626          35 :   long vx = gvar(A);
    2627             : 
    2628          35 :   n = lg(A)-1;
    2629          35 :   if (vx==NO_VARIABLE || !n)
    2630             :   {
    2631           0 :     RgM_check_ZM(A, "mathnf0");
    2632           0 :     return ZM_hnfall(A, pB, remove);
    2633             :   }
    2634          35 :   m = nbrows(A);
    2635          35 :   av = avma;
    2636          35 :   A = RgM_shallowcopy(A);
    2637          35 :   B = pB? matid(n): NULL;
    2638          35 :   def = n; ldef = (m>n)? m-n: 0;
    2639          77 :   for (li=m; li>ldef; li--)
    2640             :   {
    2641             :     GEN d, T;
    2642          63 :     for (j=def-1; j; j--)
    2643             :     {
    2644          21 :       GEN a = gcoeff(A,li,j);
    2645          21 :       if (gequal0(a)) continue;
    2646             : 
    2647           7 :       k = (j==1)? def: j-1;
    2648           7 :       RgC_elem(a,gcoeff(A,li,k), A,B, j,k, li, vx);
    2649             :     }
    2650          42 :     T = normalize_as_RgX(gcoeff(A,li,def), vx, &d);
    2651          42 :     if (gequal0(T))
    2652           0 :     { if (ldef) ldef--; }
    2653             :     else
    2654             :     {
    2655          42 :       gcoeff(A,li,def) = T;
    2656          42 :       if (B && !gequal1(d)) gel(B, def) = RgC_Rg_div(gel(B, def), d);
    2657          42 :       RgM_reduce(A, B, li, def, vx);
    2658          42 :       def--;
    2659             :     }
    2660          42 :     if (gc_needed(av,1))
    2661             :     {
    2662           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"ghnfall");
    2663           0 :       gerepileall(av, B? 2: 1, &A, &B);
    2664             :     }
    2665             :   }
    2666             :   /* rank A = n - def */
    2667          35 :   if (remove) remove_0cols(def, &A, &B, remove);
    2668          35 :   gerepileall(av, B? 2: 1, &A, &B);
    2669          35 :   if (B) *pB = B;
    2670          35 :   return A;
    2671             : }
    2672             : 
    2673             : static GEN
    2674          35 : RgXM_snf(GEN x,long all)
    2675             : {
    2676             :   pari_sp av;
    2677             :   long i, j, k, n;
    2678             :   GEN z, u, v, U, V;
    2679          35 :   long vx = gvar(x);
    2680          35 :   n = lg(x)-1; if (!n) return trivsmith(all);
    2681          35 :   if (vx==NO_VARIABLE) pari_err_TYPE("RgXM_snf",x);
    2682          35 :   if (lgcols(x) != n+1) pari_err_DIM("gsmithall");
    2683          35 :   av = avma;
    2684          35 :   x = RgM_shallowcopy(x);
    2685          35 :   if (all) { U = matid(n); V = matid(n); }
    2686         252 :   for (i=n; i>=2; i--)
    2687             :   {
    2688             :     for(;;)
    2689          98 :     {
    2690             :       GEN a, b, d;
    2691         189 :       int c = 0;
    2692         546 :       for (j=i-1; j>=1; j--)
    2693             :       {
    2694         357 :         b = gcoeff(x,i,j); if (gequal0(b)) continue;
    2695         140 :         a = gcoeff(x,i,i);
    2696         140 :         d = gbezout_step(&b, &a, &v, &u, vx);
    2697         469 :         for (k = 1; k < i; k++)
    2698             :         {
    2699         329 :           GEN t = gadd(gmul(u,gcoeff(x,k,i)),gmul(v,gcoeff(x,k,j)));
    2700         329 :           gcoeff(x,k,j) = gsub(gmul(a,gcoeff(x,k,j)),gmul(b,gcoeff(x,k,i)));
    2701         329 :           gcoeff(x,k,i) = t;
    2702             :         }
    2703         140 :         gcoeff(x,i,j) = gen_0;
    2704         140 :         gcoeff(x,i,i) = d;
    2705         140 :         if (all) update(u,v,a,b,(GEN*)(V+i),(GEN*)(V+j));
    2706             :       }
    2707         546 :       for (j=i-1; j>=1; j--)
    2708             :       {
    2709         357 :         b = gcoeff(x,j,i); if (gequal0(b)) continue;
    2710         119 :         a = gcoeff(x,i,i);
    2711         119 :         d = gbezout_step(&b, &a, &v, &u, vx);
    2712         406 :         for (k = 1; k < i; k++)
    2713             :         {
    2714         287 :           GEN t = gadd(gmul(u,gcoeff(x,i,k)),gmul(v,gcoeff(x,j,k)));
    2715         287 :           gcoeff(x,j,k) = gsub(gmul(a,gcoeff(x,j,k)),gmul(b,gcoeff(x,i,k)));
    2716         287 :           gcoeff(x,i,k) = t;
    2717             :         }
    2718         119 :         gcoeff(x,j,i) = gen_0;
    2719         119 :         gcoeff(x,i,i) = d;
    2720         119 :         if (all) update(u,v,a,b,(GEN*)(U+i),(GEN*)(U+j));
    2721         119 :         c = 1;
    2722             :       }
    2723         189 :       if (!c)
    2724             :       {
    2725         112 :         k = gsnf_no_divide(x, i, vx);
    2726         112 :         if (!k) break;
    2727             : 
    2728          70 :         for (j=1; j<=i; j++)
    2729          49 :           gcoeff(x,i,j) = gadd(gcoeff(x,i,j),gcoeff(x,k,j));
    2730          21 :         if (all) gel(U,i) = gadd(gel(U,i),gel(U,k));
    2731             :       }
    2732          98 :       if (gc_needed(av,1))
    2733             :       {
    2734           0 :         if (DEBUGMEM>1) pari_warn(warnmem,"gsmithall");
    2735           0 :         gerepileall(av, all? 3: 1, &x, &U, &V);
    2736             :       }
    2737             :     }
    2738             :   }
    2739         161 :   for (k=1; k<=n; k++)
    2740             :   {
    2741         126 :     GEN d, T = normalize_as_RgX(gcoeff(x,k,k), vx, &d);
    2742         126 :     if (gequal0(T)) continue;
    2743         112 :     if (all && !gequal1(d)) gel(V,k) = RgC_Rg_div(gel(V,k), d);
    2744         112 :     gcoeff(x,k,k) = T;
    2745             :   }
    2746          35 :   z = all? mkvec3(shallowtrans(U), V, x): RgM_diagonal_shallow(x);
    2747          35 :   return gerepilecopy(av, z);
    2748             : }
    2749             : 
    2750             : GEN
    2751         392 : matsnf0(GEN x,long flag)
    2752             : {
    2753         392 :   pari_sp av = avma;
    2754         392 :   if (flag > 7) pari_err_FLAG("matsnf");
    2755         392 :   if (typ(x) == t_VEC && flag & 4) return smithclean(x);
    2756         392 :   if (typ(x)!=t_MAT) pari_err_TYPE("matsnf",x);
    2757         392 :   if (RgM_is_ZM(x)) x = flag&1 ? smithall(x): smith(x);
    2758          35 :   else              x = RgXM_snf(x, flag&1);
    2759         392 :   if (flag & 4) x = gerepileupto(av, smithclean(x));
    2760         392 :   return x;
    2761             : }
    2762             : GEN
    2763           0 : gsmith(GEN x) { return RgXM_snf(x,0); }
    2764             : GEN
    2765           0 : gsmithall(GEN x) { return RgXM_snf(x,1); }
    2766             : 
    2767             : /* H is a relation matrix, either in HNF or a t_VEC (diagonal HNF) */
    2768             : static GEN
    2769      192555 : snf_group(GEN H, GEN D, GEN *newU, GEN *newUi)
    2770             : {
    2771             :   long i, j, l;
    2772             : 
    2773      192555 :   ZM_snfclean(D, newU? *newU: NULL, newUi? *newUi: NULL);
    2774      192555 :   l = lg(D);
    2775      192555 :   if (newU) {
    2776      175966 :     GEN U = *newU;
    2777      540610 :     for (i = 1; i < l; i++)
    2778             :     {
    2779      364644 :       GEN d = gel(D,i), d2 = shifti(d, 1);
    2780     1773219 :       for (j = 1; j < lg(U); j++)
    2781     1408575 :         gcoeff(U,i,j) = centermodii(gcoeff(U,i,j), d, d2);
    2782             :     }
    2783      175966 :     *newU = U;
    2784             :   }
    2785      192555 :   if (newUi && l > 1)
    2786             :   { /* UHV=D -> U^-1 = (HV)D^-1 -> U^-1 = H(VD^-1 mod 1) mod H */
    2787             :     /* Ui = ZM_inv(U, NULL); setlg(Ui, l); */
    2788      184938 :     GEN V = *newUi, Ui;
    2789      184938 :     int Hvec = (typ(H) == t_VEC);
    2790      184938 :     for (i = 1; i < l; i++) gel(V,i) = FpC_red(gel(V,i), gel(D,i));
    2791      184938 :     if (!Hvec)
    2792             :     {
    2793       84397 :       if (ZM_isdiagonal(H)) { H = RgM_diagonal_shallow(H); Hvec = 1; }
    2794             :     }
    2795      184938 :     Ui = Hvec? ZM_diag_mul(H, V): ZM_mul(H, V);
    2796      184938 :     for (i = 1; i < l; i++) gel(Ui,i) = ZC_Z_divexact(gel(Ui,i), gel(D,i));
    2797      184938 :     if (Hvec)
    2798      135715 :     { for (i = 1; i < l; i++) gel(Ui,i) = vecmodii(gel(Ui,i), H); }
    2799             :     else
    2800       49223 :       Ui = ZM_hnfrem(Ui, H);
    2801      184938 :     *newUi = Ui;
    2802             :   }
    2803      192555 :   return D;
    2804             : }
    2805             : /* H relation matrix among row of generators g in HNF.  Let URV = D its SNF,
    2806             :  * newU R newV = newD its clean SNF (no 1 in Dnew). Return the diagonal of
    2807             :  * newD, newU and newUi such that  1/U = (newUi, ?).
    2808             :  * Rationale: let (G,0) = g Ui be the new generators then
    2809             :  * 0 = G U R --> G D = 0,  g = G newU  and  G = g newUi */
    2810             : GEN
    2811       91825 : ZM_snf_group(GEN H, GEN *newU, GEN *newUi)
    2812             : {
    2813       91825 :   GEN D = ZM_snfall_i(H, newU, newUi, 1 + 2);
    2814       91825 :   return snf_group(H, D, newU, newUi);
    2815             : }
    2816             : 
    2817             : /* D a ZV: SNF for matdiagonal(D). Faster because we only ensure elementary
    2818             :  * divisors condition: d[n] | ... | d[1] and need not convert D to matrix form*/
    2819             : GEN
    2820      100730 : ZV_snfall(GEN D, GEN *pU, GEN *pV)
    2821             : {
    2822      100730 :   pari_sp av = avma;
    2823      100730 :   long j, n = lg(D)-1;
    2824      100730 :   GEN U = pU? matid(n): NULL;
    2825      100730 :   GEN V = pV? matid(n): NULL;
    2826             :   GEN p;
    2827             : 
    2828      100730 :   D = leafcopy(D);
    2829      329567 :   for (j = n; j > 0; j--)
    2830             :   {
    2831      228837 :     GEN b = gel(D,j);
    2832      228837 :     if (signe(b) < 0)
    2833             :     {
    2834           0 :       gel(D,j) = negi(b);
    2835           0 :       if (V) ZV_togglesign(gel(V,j));
    2836             :     }
    2837             :   }
    2838             :   /* entries are non-negative integers */
    2839      100730 :   p = gen_indexsort(D, NULL, &negcmpii);
    2840      100730 :   D = vecpermute(D, p);
    2841      100730 :   if (U) U = vecpermute(U, p);
    2842      100730 :   if (V) V = vecpermute(V, p);
    2843             :   /* entries are sorted by decreasing value */
    2844      329567 :   for (j = n; j > 0; j--)
    2845             :   {
    2846      228837 :     GEN b = gel(D,j);
    2847             :     long i;
    2848      421372 :     for (i = j-1; i > 0; i--)
    2849             :     { /* invariant: a >= b. If au+bv = d is a Bezout relation, A=a/d and B=b/d
    2850             :        * we have [B,-A;u,v]*diag(a,b)*[1-u*A,1; -u*A,1]] = diag(Ab, d) */
    2851      193739 :       GEN a = gel(D,i), u,v, d = bezout(a,b, &u,&v), A, Wi, Wj;
    2852      193739 :       if (equalii(d,b)) continue;
    2853       25438 :       A = diviiexact(a,d);
    2854       25438 :       if (V)
    2855             :       {
    2856       25396 :         GEN t = mulii(u,A);
    2857       25396 :         Wi = ZC_lincomb(subui(1,t), negi(t), gel(V,i), gel(V,j));
    2858       25396 :         Wj = ZC_add(gel(V,i), gel(V,j));
    2859       25396 :         gel(V,i) = Wi;
    2860       25396 :         gel(V,j) = Wj;
    2861             :       }
    2862       25438 :       if (U)
    2863             :       {
    2864       25438 :         GEN B = diviiexact(b,d);
    2865       25438 :         Wi = ZC_lincomb(B, negi(A), gel(U,i), gel(U,j));
    2866       25438 :         Wj = ZC_lincomb(u, v, gel(U,i), gel(U,j));
    2867       25438 :         gel(U,i) = Wi;
    2868       25438 :         gel(U,j) = Wj;
    2869             :       }
    2870       25438 :       gel(D,i) = mulii(A,b); /* lcm(a,b) */
    2871       25438 :       gel(D,j) = d; /* gcd(a,b) */
    2872       25438 :       b = gel(D,j); if (equali1(b)) break;
    2873             :     }
    2874             :   }
    2875      100730 :   snf_pile(av, &D,&U,&V);
    2876      100730 :   if (U) *pU = shallowtrans(U);
    2877      100730 :   if (V) *pV = V;
    2878      100730 :   return D;
    2879             : }
    2880             : GEN
    2881      100730 : ZV_snf_group(GEN d, GEN *newU, GEN *newUi)
    2882             : {
    2883      100730 :   GEN D = ZV_snfall(d, newU, newUi);
    2884      100730 :   return snf_group(d, D, newU, newUi);
    2885             : }
    2886             : 
    2887             : /* D a vector of elementary divisors. Truncate (setlg) to leave out trivial
    2888             :  * entries (= 1) */
    2889             : void
    2890           0 : ZV_snf_trunc(GEN D)
    2891             : {
    2892           0 :   long i, l = lg(D);
    2893           0 :   for (i = 1; i < l; i++)
    2894           0 :     if (is_pm1(gel(D,i))) { setlg(D,i); break; }
    2895           0 : }

Generated by: LCOV version 1.13