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.18.0 lcov report (development 29792-e7296d7d08) Lines: 1612 1788 90.2 %
Date: 2024-12-05 09:09:40 Functions: 97 107 90.7 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2000  The PARI group.
       2             : 
       3             : This file is part of the PARI/GP package.
       4             : 
       5             : PARI/GP is free software; you can redistribute it and/or modify it under the
       6             : terms of the GNU General Public License as published by the Free Software
       7             : Foundation; either version 2 of the License, or (at your option) any later
       8             : version. It is distributed in the hope that it will be useful, but WITHOUT
       9             : ANY WARRANTY WHATSOEVER.
      10             : 
      11             : Check the License for details. You should have received a copy of it, along
      12             : with the package; see the file 'COPYING'. If not, write to the Free Software
      13             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      14             : 
      15             : #include "pari.h"
      16             : #include "paripriv.h"
      17             : 
      18             : #define DEBUGLEVEL DEBUGLEVEL_mathnf
      19             : 
      20             : /**************************************************************/
      21             : /**                                                          **/
      22             : /**                HERMITE NORMAL FORM REDUCTION             **/
      23             : /**                                                          **/
      24             : /**************************************************************/
      25             : static GEN ZV_hnfgcdext(GEN A);
      26             : static GEN
      27          21 : hnfallgen(GEN x)
      28             : {
      29          21 :   GEN z = cgetg(3, t_VEC);
      30          21 :   gel(z,1) = RgM_hnfall(x, (GEN*)(z+2), 1);
      31          21 :   return z;
      32             : }
      33             : GEN
      34         287 : mathnf0(GEN x, long flag)
      35             : {
      36         287 :   switch(typ(x))
      37             :   {
      38          70 :     case t_VEC:
      39          70 :       if (RgV_is_ZV(x))
      40             :         switch (flag)
      41             :         {
      42          14 :           case 0:
      43          14 :             if (lg(x) == 1) return cgetg(1, t_MAT);
      44           7 :             retmkmat(mkcol(ZV_content(x)));
      45          21 :           case 1:
      46             :           case 4:
      47          21 :             return ZV_hnfgcdext(x);
      48             :         }
      49          35 :       x = gtomat(x); break;
      50         217 :     case t_MAT: break;
      51           0 :     default: pari_err_TYPE("mathnf0",x);
      52             :   }
      53             : 
      54         252 :   switch(flag)
      55             :   {
      56         196 :     case 0: case 2: return RgM_is_ZM(x)? ZM_hnf(x): RgM_hnfall(x,NULL,1);
      57          35 :     case 1: case 3: return RgM_is_ZM(x)? hnfall(x): hnfallgen(x);
      58           7 :     case 4: RgM_check_ZM(x, "mathnf0"); return hnflll(x);
      59          14 :     case 5: RgM_check_ZM(x, "mathnf0"); return hnfperm(x);
      60           0 :     default: pari_err_FLAG("mathnf");
      61             :   }
      62             :   return NULL; /* LCOV_EXCL_LINE */
      63             : }
      64             : 
      65             : /*******************************************************************/
      66             : /*                                                                 */
      67             : /*                SPECIAL HNF (FOR INTERNAL USE !!!)               */
      68             : /*                                                                 */
      69             : /*******************************************************************/
      70             : static int
      71     8101449 : count(GEN mat, long row, long len, long *firstnonzero)
      72             : {
      73     8101449 :   long j, n = 0;
      74             : 
      75   731180669 :   for (j=1; j<=len; j++)
      76             :   {
      77   724860146 :     long p = mael(mat,j,row);
      78   724860146 :     if (p)
      79             :     {
      80    23025347 :       if (labs(p)!=1) return -1;
      81    21244421 :       n++; *firstnonzero=j;
      82             :     }
      83             :   }
      84     6320523 :   return n;
      85             : }
      86             : 
      87             : static int
      88      420495 : count2(GEN mat, long row, long len)
      89             : {
      90             :   long j;
      91     4444696 :   for (j=len; j; j--)
      92     4315657 :     if (labs(mael(mat,j,row)) == 1) return j;
      93      129039 :   return 0;
      94             : }
      95             : 
      96             : static GEN
      97      225509 : hnffinal(GEN matgen,GEN perm,GEN* ptdep,GEN* ptB,GEN* ptC)
      98             : {
      99             :   GEN p1,p2,U,H,Hnew,Bnew,Cnew,diagH1;
     100      225509 :   GEN B = *ptB, C = *ptC, dep = *ptdep, depnew;
     101             :   pari_sp av;
     102             :   long i,j,k,s,i1,j1,zc;
     103      225509 :   long co = lg(C);
     104      225509 :   long col = lg(matgen)-1;
     105             :   long lnz, nlze, lig;
     106             : 
     107      225509 :   if (col == 0) return matgen;
     108      225509 :   lnz = nbrows(matgen); /* was called lnz-1 - nr in hnfspec */
     109      225509 :   nlze = nbrows(dep);
     110      225509 :   lig = nlze + lnz;
     111             :   /* H: lnz x lnz [disregarding initial 0 cols], U: col x col */
     112      225509 :   H = ZM_hnflll(matgen, &U, 0);
     113      225509 :   H += lg(H)-1 - lnz; H[0] = evaltyp(t_MAT) | _evallg(lnz+1);
     114             :   /* Only keep the part above the H (above the 0s is 0 since the dep rows
     115             :    * are dependent from the ones in matgen) */
     116      225509 :   zc = col - lnz; /* # of 0 columns, correspond to units */
     117      225509 :   if (nlze) { dep = ZM_mul(dep,U); dep += zc; }
     118             : 
     119      225509 :   diagH1 = new_chunk(lnz+1); /* diagH1[i] = 0 iff H[i,i] != 1 (set later) */
     120             : 
     121      225509 :   av = avma;
     122      225509 :   Cnew = cgetg(co, typ(C));
     123      225509 :   setlg(C, col+1); p1 = gmul(C,U);
     124     1899917 :   for (j=1; j<=col; j++) gel(Cnew,j) = gel(p1,j);
     125     3194092 :   for (   ; j<co ; j++)  gel(Cnew,j) = gel(C,j);
     126             : 
     127             :   /* Clean up B using new H */
     128      895810 :   for (s=0,i=lnz; i; i--)
     129             :   {
     130      670327 :     GEN Di = gel(dep,i), Hi = gel(H,i);
     131      670327 :     GEN h = gel(Hi,i); /* H[i,i] */
     132      670327 :     if ( (diagH1[i] = is_pm1(h)) ) { h = NULL; s++; }
     133    16383232 :     for (j=col+1; j<co; j++)
     134             :     {
     135    15713293 :       GEN z = gel(B,j-col);
     136    15713293 :       p1 = gel(z,i+nlze);
     137    15713293 :       if (h) p1 = truedivii(p1,h);
     138    15713940 :       if (!signe(p1)) continue;
     139    27612746 :       for (k=1; k<=nlze; k++) gel(z,k) = subii(gel(z,k), mulii(p1, gel(Di,k)));
     140   134359062 :       for (   ; k<=lig;  k++) gel(z,k) = subii(gel(z,k), mulii(p1, gel(Hi,k-nlze)));
     141     9311168 :       gel(Cnew,j) = gsub(gel(Cnew,j), gmul(p1, gel(Cnew,i+zc)));
     142             :     }
     143      669939 :     if (gc_needed(av,2))
     144             :     {
     145         746 :       if(DEBUGMEM>1) pari_warn(warnmem,"hnffinal, i = %ld",i);
     146         746 :       gerepileall(av, 2, &Cnew, &B);
     147             :     }
     148             :   }
     149      225483 :   p1 = cgetg(lnz+1,t_VEC); p2 = perm + nlze;
     150      895842 :   for (i1=0, j1=lnz-s, i=1; i<=lnz; i++) /* push the 1 rows down */
     151      670340 :     if (diagH1[i])
     152      305215 :       gel(p1,++j1) = gel(p2,i);
     153             :     else
     154      365125 :       gel(p2,++i1) = gel(p2,i);
     155      530718 :   for (i=i1+1; i<=lnz; i++) gel(p2,i) = gel(p1,i);
     156             : 
     157             :   /* s = # extra redundant generators taken from H
     158             :    *          zc  col-s  co   zc = col - lnz
     159             :    *       [ 0 |dep |     ]    i = nlze + lnz - s = lig - s
     160             :    *  nlze [--------|  B' ]
     161             :    *       [ 0 | H' |     ]    H' = H minus the s rows with a 1 on diagonal
     162             :    *     i [--------|-----] lig-s           (= "1-rows")
     163             :    *       [   0    | Id  ]
     164             :    *       [        |     ] li */
     165      225502 :   lig -= s; col -= s; lnz -= s;
     166      225502 :   Hnew = cgetg(lnz+1,t_MAT);
     167      225503 :   depnew = cgetg(lnz+1,t_MAT); /* only used if nlze > 0 */
     168      225503 :   Bnew = cgetg(co-col,t_MAT);
     169      225505 :   C = shallowcopy(Cnew);
     170      895865 :   for (j=1,i1=j1=0; j<=lnz+s; j++)
     171             :   {
     172      670358 :     GEN z = gel(H,j);
     173      670358 :     if (diagH1[j])
     174             :     { /* hit exactly s times */
     175      305219 :       i1++; C[i1+col] = Cnew[j+zc];
     176      305219 :       p1 = cgetg(lig+1,t_COL); gel(Bnew,i1) = p1;
     177      535913 :       for (i=1; i<=nlze; i++) gel(p1,i) = gcoeff(dep,i,j);
     178      305218 :       p1 += nlze;
     179             :     }
     180             :     else
     181             :     {
     182      365139 :       j1++; C[j1+zc] = Cnew[j+zc];
     183      365139 :       p1 = cgetg(lnz+1,t_COL); gel(Hnew,j1) = p1;
     184      365136 :       depnew[j1] = dep[j];
     185             :     }
     186     2563977 :     for (i=k=1; k<=lnz; i++)
     187     1893623 :       if (!diagH1[i]) p1[k++] = z[i];
     188             :   }
     189     3194002 :   for (j=s+1; j<co-col; j++)
     190             :   {
     191     2968505 :     GEN z = gel(B,j-s);
     192     2968505 :     p1 = cgetg(lig+1,t_COL); gel(Bnew,j) = p1;
     193     5209430 :     for (i=1; i<=nlze; i++) gel(p1,i) = gel(z,i);
     194     2968495 :     z += nlze; p1 += nlze;
     195    12979892 :     for (i=k=1; k<=lnz; i++)
     196    10011397 :       if (!diagH1[i]) gel(p1,k++) = gel(z,i);
     197             :   }
     198      225497 :   *ptdep = depnew;
     199      225497 :   *ptC = C;
     200      225497 :   *ptB = Bnew; return Hnew;
     201             : }
     202             : 
     203             : /* for debugging */
     204             : static void
     205           0 : p_mat(GEN mat, GEN perm, long k)
     206             : {
     207           0 :   pari_sp av = avma;
     208           0 :   perm = vecslice(perm, k+1, lg(perm)-1);
     209           0 :   err_printf("Permutation: %Ps\n",perm);
     210           0 :   if (DEBUGLEVEL > 6)
     211           0 :     err_printf("matgen = %Ps\n", zm_to_ZM( rowpermute(mat, perm) ));
     212           0 :   set_avma(av);
     213           0 : }
     214             : 
     215             : static GEN
     216     2609161 : col_dup(long l, GEN col)
     217             : {
     218     2609161 :   GEN c = new_chunk(l);
     219     2609140 :   memcpy(c,col,l * sizeof(long)); return c;
     220             : }
     221             : 
     222             : /* permutation giving imagecompl(x') | image(x'), x' = transpose of x */
     223             : static GEN
     224      225507 : ZM_rowrankprofile(GEN x, long *nlze)
     225             : {
     226      225507 :   pari_sp av = avma;
     227             :   GEN d, y;
     228             :   long i, j, k, l, r;
     229             : 
     230      225507 :   x = shallowtrans(x); l = lg(x);
     231      225509 :   (void)new_chunk(l); /* HACK */
     232      225509 :   d = ZM_pivots(x,&r); set_avma(av);
     233      225507 :   *nlze = r;
     234      225507 :   if (!d) return identity_perm(l-1);
     235      215384 :   y = cgetg(l,t_VECSMALL);
     236      893441 :   for (i = j = 1, k = r+1; i<l; i++)
     237      678057 :     if (d[i]) y[k++] = i; else y[j++] = i;
     238      215384 :   return y;
     239             : }
     240             : 
     241             : /* HNF reduce a relation matrix (column operations + row permutation)
     242             : ** Input:
     243             : **   mat = (li-1) x (co-1) matrix of long
     244             : **   C   = r x (co-1) matrix of GEN
     245             : **   perm= permutation vector (length li-1), indexing the rows of mat: easier
     246             : **     to maintain perm than to copy rows. For columns we can do it directly
     247             : **     using e.g. swap(mat[i], mat[j])
     248             : **   k0 = integer. The k0 first lines of mat are dense, the others are sparse.
     249             : ** Output: cf ASCII art in the function body
     250             : **
     251             : ** row permutations applied to perm
     252             : ** column operations applied to C. IN PLACE
     253             : **/
     254             : GEN
     255      134279 : hnfspec_i(GEN mat0, GEN perm, GEN* ptdep, GEN* ptB, GEN* ptC, long k0)
     256             : {
     257             :   pari_sp av;
     258             :   long co, n, s, nlze, lnz, nr, i, j, k, lk0, col, lig, *p;
     259             :   GEN mat;
     260             :   GEN p1, p2, matb, matbnew, vmax, matt, T, extramat, B, C, H, dep, permpro;
     261      134279 :   const long li = lg(perm); /* = lgcols(mat0) */
     262      134279 :   const long CO = lg(mat0);
     263             : 
     264      134279 :   n = 0; /* -Wall */
     265             : 
     266      134279 :   C = *ptC; co = CO;
     267      134279 :   if (co > 300 && co > 1.5 * li)
     268             :   { /* treat the rest at the end */
     269           0 :     co = (long)(1.2 * li);
     270           0 :     setlg(C, co);
     271             :   }
     272             : 
     273      134279 :   if (DEBUGLEVEL>5)
     274             :   {
     275           0 :     err_printf("Entering hnfspec\n");
     276           0 :     p_mat(mat0,perm,0);
     277             :   }
     278      134279 :   matt = cgetg(co, t_MAT); /* dense part of mat (top) */
     279      134279 :   mat = cgetg(co, t_MAT);
     280     2743430 :   for (j = 1; j < co; j++)
     281             :   {
     282     2609150 :     GEN matj = col_dup(li, gel(mat0,j));
     283     2609148 :     p1 = cgetg(k0+1,t_COL); gel(matt,j) = p1; gel(mat,j) = matj;
     284    12286367 :     for (i=1; i<=k0; i++) gel(p1,i) = stoi(matj[perm[i]]);
     285             :   }
     286      134280 :   av = avma;
     287             : 
     288      134280 :   i = lig = li-1; col = co-1; lk0 = k0;
     289      134280 :   T = (k0 || (lg(C) > 1 && lgcols(C) > 1))? matid(col): NULL;
     290             :   /* Look for lines with a single nonzero entry, equal to 1 in absolute value */
     291     6386287 :   while (i > lk0 && col)
     292     6251988 :     switch( count(mat,perm[i],col,&n) )
     293             :     {
     294       26247 :       case 0: /* move zero lines between k0+1 and lk0 */
     295       26247 :         lk0++; lswap(perm[i], perm[lk0]);
     296       26247 :         i = lig; continue;
     297             : 
     298      484424 :       case 1: /* move trivial generator between lig+1 and li */
     299      484424 :         lswap(perm[i], perm[lig]);
     300      484424 :         if (T) swap(gel(T,n), gel(T,col));
     301      484424 :         swap(gel(mat,n), gel(mat,col)); p = gel(mat,col);
     302      484424 :         if (p[perm[lig]] < 0) /* = -1 */
     303             :         { /* convert relation -g = 0 to g = 0 */
     304    10572022 :           for (i=lk0+1; i<lig; i++) p[perm[i]] = -p[perm[i]];
     305      436995 :           if (T)
     306             :           {
     307      436995 :             p1 = gel(T,col);
     308    11271973 :             for (i=1; ; i++) /* T = permuted identity: single nonzero entry */
     309    11271973 :               if (signe(gel(p1,i))) { togglesign_safe(&gel(p1,i)); break; }
     310             :           }
     311             :         }
     312      484443 :         lig--; col--; i = lig; continue;
     313             : 
     314     5741317 :       default: i--;
     315             :     }
     316      134299 :   if (DEBUGLEVEL>5) { err_printf("    after phase1:\n"); p_mat(mat,perm,0); }
     317             : 
     318             : #define absmax(s,z) {long _z; _z = labs(z); if (_z > s) s = _z;}
     319             :   /* Get rid of all lines containing only 0 and +/- 1, keeping track of column
     320             :    * operations in T. Leave the rows 1..lk0 alone [up to k0, coefficient
     321             :    * explosion, between k0+1 and lk0, row is 0] */
     322      134278 :   s = 0;
     323      775796 :   while (lig > lk0 && col && s < (long)(HIGHBIT>>1))
     324             :   {
     325     1974430 :     for (i=lig; i>lk0; i--)
     326     1849516 :       if (count(mat,perm[i],col,&n) > 0) break;
     327      766470 :     if (i == lk0) break;
     328             : 
     329             :     /* only 0, +/- 1 entries, at least 2 of them nonzero */
     330      641548 :     lswap(perm[i], perm[lig]);
     331      641548 :     swap(gel(mat,n), gel(mat,col)); p = gel(mat,col);
     332      641548 :     if (T) swap(gel(T,n), gel(T,col));
     333      641548 :     if (p[perm[lig]] < 0)
     334             :     {
     335     7802351 :       for (i=lk0+1; i<=lig; i++) p[perm[i]] = -p[perm[i]];
     336      416978 :       if (T) ZV_togglesign(gel(T,col));
     337             :     }
     338    21443518 :     for (j=1; j<col; j++)
     339             :     {
     340    20802000 :       GEN matj = gel(mat,j);
     341             :       long t;
     342    20802000 :       if (! (t = matj[perm[lig]]) ) continue;
     343     1512801 :       if (t == 1) {
     344    26307690 :         for (i=lk0+1; i<=lig; i++) absmax(s, matj[perm[i]] -= p[perm[i]]);
     345             :       }
     346             :       else { /* t = -1 */
     347    19414952 :         for (i=lk0+1; i<=lig; i++) absmax(s, matj[perm[i]] += p[perm[i]]);
     348             :       }
     349     1512801 :       if (T) ZC_lincomb1_inplace(gel(T,j), gel(T,col), stoi(-t));
     350             :     }
     351      641518 :     lig--; col--;
     352      641518 :     if (gc_needed(av,3))
     353             :     {
     354           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"hnfspec[1]");
     355           0 :       if (T) T = gerepilecopy(av, T); else set_avma(av);
     356             :     }
     357             :   }
     358             :   /* As above with lines containing a +/- 1 (no other assumption).
     359             :    * Stop when single precision becomes dangerous */
     360      134251 :   vmax = cgetg(co,t_VECSMALL);
     361     1617538 :   for (j=1; j<=col; j++)
     362             :   {
     363     1483258 :     GEN matj = gel(mat,j);
     364     7080312 :     for (s=0, i=lk0+1; i<=lig; i++) absmax(s, matj[i]);
     365     1483258 :     vmax[j] = s;
     366             :   }
     367      425736 :   while (lig > lk0 && col)
     368             :   {
     369      448047 :     for (i=lig; i>lk0; i--)
     370      420495 :       if ( (n = count2(mat,perm[i],col)) ) break;
     371      319007 :     if (i == lk0) break;
     372             : 
     373      291454 :     lswap(vmax[n], vmax[col]);
     374      291454 :     lswap(perm[i], perm[lig]);
     375      291454 :     swap(gel(mat,n), gel(mat,col)); p = gel(mat,col);
     376      291454 :     if (T) swap(gel(T,n), gel(T,col));
     377      291454 :     if (p[perm[lig]] < 0)
     378             :     {
     379      593489 :       for (i=lk0+1; i<=lig; i++) p[perm[i]] = -p[perm[i]];
     380       95657 :       if (T) ZV_togglesign(gel(T,col));
     381             :     }
     382     3893122 :     for (j=1; j<col; j++)
     383             :     {
     384     3601711 :       GEN matj = gel(mat,j);
     385             :       long t;
     386     3601711 :       if (! (t = matj[perm[lig]]) ) continue;
     387     1723445 :       if (vmax[col] && (ulong)labs(t) >= (HIGHBIT-vmax[j]) / vmax[col])
     388           0 :         goto END2;
     389             : 
     390    20122541 :       for (s=0, i=lk0+1; i<=lig; i++) absmax(s, matj[perm[i]] -= t*p[perm[i]]);
     391     1723445 :       vmax[j] = s;
     392     1723445 :       if (T) ZC_lincomb1_inplace(gel(T,j), gel(T,col), stoi(-t));
     393             :     }
     394      291411 :     lig--; col--;
     395      291411 :     if (gc_needed(av,3))
     396             :     {
     397           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"hnfspec[2]");
     398           0 :       gerepileall(av, T? 2: 1, &vmax, &T);
     399             :     }
     400             :   }
     401             : 
     402      106728 : END2: /* clean up mat: remove everything to the right of the 1s on diagonal */
     403             :   /* go multiprecision first */
     404      134281 :   matb = cgetg(co,t_MAT); /* bottom part (complement of matt) */
     405     2743382 :   for (j=1; j<co; j++)
     406             :   {
     407     2609136 :     GEN matj = gel(mat,j);
     408     2609136 :     p1 = cgetg(li-k0,t_COL); gel(matb,j) = p1;
     409     2609108 :     p1 -= k0;
     410    82473928 :     for (i=k0+1; i<li; i++) gel(p1,i) = stoi(matj[perm[i]]);
     411             :   }
     412      134246 :   if (DEBUGLEVEL>5)
     413             :   {
     414           0 :     err_printf("    after phase2:\n");
     415           0 :     p_mat(mat,perm,lk0);
     416             :   }
     417     1417874 :   for (i=li-2; i>lig; i--)
     418             :   {
     419     1283600 :     long h, i0 = i - k0, k = i + co-li;
     420     1283600 :     GEN Bk = gel(matb,k);
     421    29087295 :     for (j=k+1; j<co; j++)
     422             :     {
     423    27803947 :       GEN Bj = gel(matb,j), v = gel(Bj,i0);
     424    27803947 :       s = signe(v); if (!s) continue;
     425             : 
     426     5391041 :       gel(Bj,i0) = gen_0;
     427     5391041 :       if (is_pm1(v))
     428             :       {
     429     3081972 :         if (s > 0) /* v = 1 */
     430    46869584 :         { for (h=1; h<i0; h++) gel(Bj,h) = subii(gel(Bj,h), gel(Bk,h)); }
     431             :         else /* v = -1 */
     432    40313091 :         { for (h=1; h<i0; h++) gel(Bj,h) = addii(gel(Bj,h), gel(Bk,h)); }
     433             :       }
     434             :       else {
     435    55615827 :         for (h=1; h<i0; h++) gel(Bj,h) = subii(gel(Bj,h), mulii(v,gel(Bk,h)));
     436             :       }
     437     5390763 :       if (T) ZC_lincomb1_inplace(gel(T,j), gel(T,k), negi(v));
     438             :     }
     439     1283348 :     if (gc_needed(av,2))
     440             :     {
     441           6 :       if(DEBUGMEM>1) pari_warn(warnmem,"hnfspec[3], i = %ld", i);
     442        1590 :       for (h=1; h<co; h++) setlg(matb[h], i0+1); /* bottom can be forgotten */
     443           6 :       gerepileall(av, T? 2: 1, &matb, &T);
     444             :     }
     445             :   }
     446     2743499 :   for (j=1; j<co; j++) setlg(matb[j], lig-k0+1); /* bottom can be forgotten */
     447      134281 :   gerepileall(av, T? 2: 1, &matb, &T);
     448      134280 :   if (DEBUGLEVEL>5) err_printf("    matb cleaned up (using Id block)\n");
     449             : 
     450      134280 :   nlze = lk0 - k0;  /* # of 0 rows */
     451      134280 :   lnz = lig-nlze+1; /* 1 + # of nonzero rows (!= 0...0 1 0 ... 0) */
     452      134280 :   if (T) matt = ZM_mul(matt,T); /* update top rows */
     453      134276 :   extramat = cgetg(col+1,t_MAT); /* = new C minus the 0 rows */
     454     1326036 :   for (j=1; j<=col; j++)
     455             :   {
     456     1191762 :     GEN z = gel(matt,j);
     457     1191762 :     GEN t = (gel(matb,j)) + nlze - k0;
     458     1191762 :     p2=cgetg(lnz,t_COL); gel(extramat,j) = p2;
     459     5094584 :     for (i=1; i<=k0; i++) gel(p2,i) = gel(z,i); /* top k0 rows */
     460     1989582 :     for (   ; i<lnz; i++) gel(p2,i) = gel(t,i); /* other nonzero rows */
     461             :   }
     462      134274 :   if (!col) {
     463           0 :     permpro = identity_perm(lnz);
     464           0 :     nr = lnz;
     465             :   }
     466             :   else
     467      134274 :     permpro = ZM_rowrankprofile(extramat, &nr);
     468             :   /* lnz = lg(permpro) */
     469      134280 :   if (nlze)
     470             :   { /* put the nlze 0 rows (trivial generators) at the top */
     471       11594 :     p1 = new_chunk(lk0+1);
     472       37841 :     for (i=1; i<=nlze; i++) p1[i] = perm[i + k0];
     473       60249 :     for (   ; i<=lk0; i++)  p1[i] = perm[i - nlze];
     474       86496 :     for (i=1; i<=lk0; i++)  perm[i] = p1[i];
     475             :   }
     476             :   /* sort other rows according to permpro (nr redundant generators first) */
     477      134280 :   p1 = new_chunk(lnz); p2 = perm + nlze;
     478      589984 :   for (i=1; i<lnz; i++) p1[i] = p2[permpro[i]];
     479      589984 :   for (i=1; i<lnz; i++) p2[i] = p1[i];
     480             :   /* perm indexes the rows of mat
     481             :    *   |_0__|__redund__|__dense__|__too big__|_____done______|
     482             :    *   0  nlze                              lig             li
     483             :    *         \___nr___/ \___k0__/
     484             :    *         \____________lnz ______________/
     485             :    *
     486             :    *               col   co
     487             :    *       [dep     |     ]
     488             :    *    i0 [--------|  B  ] (i0 = nlze + nr)
     489             :    *       [matbnew |     ] matbnew has maximal rank = lnz-1 - nr
     490             :    * mat = [--------|-----] lig
     491             :    *       [   0    | Id  ]
     492             :    *       [        |     ] li */
     493             : 
     494      134280 :   matbnew = cgetg(col+1,t_MAT); /* dense+toobig, maximal rank. For hnffinal */
     495      134280 :   dep    = cgetg(col+1,t_MAT); /* rows dependent from the ones in matbnew */
     496     1326081 :   for (j=1; j<=col; j++)
     497             :   {
     498     1191802 :     GEN z = gel(extramat,j);
     499     1191802 :     p1 = cgetg(nlze+nr+1,t_COL); gel(dep,j) = p1;
     500     1191802 :     p2 = cgetg(lnz-nr,t_COL); gel(matbnew,j) = p2;
     501     1550979 :     for (i=1; i<=nlze; i++) gel(p1,i) = gen_0;
     502     1260576 :     p1 += nlze; for (i=1; i<=nr; i++) p1[i] = z[permpro[i]];
     503     5823802 :     p2 -= nr;   for (   ; i<lnz; i++) p2[i] = z[permpro[i]];
     504             :   }
     505             : 
     506             :   /* redundant generators in terms of the genuine generators
     507             :    * (x_i) = - (g_i) B */
     508      134279 :   B = cgetg(co-col,t_MAT);
     509     1551724 :   for (j=col+1; j<co; j++)
     510             :   {
     511     1417445 :     GEN y = gel(matt,j);
     512     1417445 :     GEN z = gel(matb,j);
     513     1417445 :     p1=cgetg(lig+1,t_COL); gel(B,j-col) = p1;
     514     2979922 :     for (i=1; i<=nlze; i++) gel(p1,i) = gel(z,i);
     515     1417444 :     p1 += nlze; z += nlze-k0;
     516    10119324 :     for (k=1; k<lnz; k++)
     517             :     {
     518     8701880 :       i = permpro[k];
     519     8701880 :       gel(p1,k) = (i <= k0)? gel(y,i): gel(z,i);
     520             :     }
     521             :   }
     522      134279 :   if (T) C = typ(C)==t_MAT? RgM_ZM_mul(C,T): RgV_RgM_mul(C,T);
     523      134278 :   gerepileall(av, 4, &matbnew, &B, &dep, &C);
     524      134280 :   *ptdep = dep;
     525      134280 :   *ptB = B;
     526      134280 :   H = hnffinal(matbnew, perm, ptdep, ptB, &C);
     527      134279 :   if (CO > co)
     528             :   { /* treat the rest, N cols at a time (hnflll slow otherwise) */
     529           0 :     const long N = 300;
     530           0 :     long a, L = CO - co, l = minss(L, N); /* L columns to add */
     531           0 :     GEN CC = *ptC, m0 = mat0;
     532           0 :     setlg(CC, CO); /* restore */
     533           0 :     CC += co-1;
     534           0 :     m0 += co-1;
     535           0 :     for (a = l;;)
     536           0 :     {
     537             :       GEN MAT, emb;
     538           0 :       gerepileall(av, 4, &H,&C,ptB,ptdep);
     539           0 :       MAT = cgetg(l + 1, t_MAT);
     540           0 :       emb = cgetg(l + 1, typ(C));
     541           0 :       for (j = 1 ; j <= l; j++)
     542             :       {
     543           0 :         gel(MAT,j) = gel(m0,j);
     544           0 :         emb[j] = CC[j];
     545             :       }
     546           0 :       H = hnfadd_i(H, perm, ptdep, ptB, &C, MAT, emb);
     547           0 :       if (a == L) break;
     548           0 :       CC += l;
     549           0 :       m0 += l;
     550           0 :       a += l; if (a > L) { l = L - (a - l); a = L; }
     551             :     }
     552             :   }
     553      134279 :   *ptC = C; return H;
     554             : }
     555             : 
     556             : GEN
     557           0 : hnfspec(GEN mat, GEN perm, GEN* ptdep, GEN* ptB, GEN* ptC, long k0)
     558             : {
     559           0 :   pari_sp av = avma;
     560           0 :   GEN H = hnfspec_i(mat, perm, ptdep, ptB, ptC, k0);
     561           0 :   return gc_all(av, 4, &H, ptC, ptdep, ptB);
     562             : }
     563             : 
     564             : /* HNF reduce x, apply same transforms to C */
     565             : GEN
     566           0 : mathnfspec(GEN x, GEN *pperm, GEN *pdep, GEN *pB, GEN *pC)
     567             : {
     568           0 :   long i, j, k, l, n, ly, lx = lg(x);
     569             :   GEN z, v1, perm;
     570           0 :   if (lx == 1) return cgetg(1, t_MAT);
     571           0 :   ly = lgcols(x);
     572           0 :   *pperm = perm = identity_perm(ly-1);
     573           0 :   z = cgetg(lx,t_MAT);
     574           0 :   for (i=1; i<lx; i++)
     575             :   {
     576           0 :     GEN C = cgetg(ly,t_COL), D = gel(x,i);
     577           0 :     gel(z,i) = C;
     578           0 :     for (j=1; j<ly; j++)
     579             :     {
     580           0 :       GEN d = gel(D,j);
     581           0 :       if (is_bigint(d)) goto TOOLARGE;
     582           0 :       C[j] = itos(d);
     583             :     }
     584             :   }
     585             :   /*  [ dep |     ]
     586             :    *  [-----|  B  ]
     587             :    *  [  H  |     ]
     588             :    *  [-----|-----]
     589             :    *  [  0  | Id  ] */
     590           0 :   return hnfspec(z,perm, pdep, pB, pC, 0);
     591             : 
     592           0 : TOOLARGE:
     593           0 :   if (lg(*pC) > 1 && lgcols(*pC) > 1)
     594           0 :     pari_err_IMPL("mathnfspec with large entries");
     595           0 :   x = ZM_hnf(x); lx = lg(x);
     596           0 :   v1 = cgetg(ly, t_VECSMALL);
     597           0 :   n = lx - ly;
     598           0 :   for (i = k = l = 1; i < ly; i++)
     599           0 :     if (equali1(gcoeff(x,i,i + n))) v1[l++] = i; else perm[k++] = i;
     600           0 :   setlg(perm, k);
     601           0 :   setlg(v1, l);
     602           0 :   x = rowpermute(x, perm); /* upper part */
     603           0 :   *pperm = vecsmall_concat(perm, v1);
     604           0 :   *pB = vecslice(x, k+n, lx-1);
     605           0 :   setlg(x, k);
     606           0 :   *pdep = rowslice(x, 1, n);
     607           0 :   return n? rowslice(x, n+1, k-1): x; /* H */
     608             : }
     609             : 
     610             : /* add new relations to a matrix treated by hnfspec (extramat / extraC) */
     611             : GEN
     612       91229 : hnfadd_i(GEN H, GEN perm, GEN* ptdep, GEN* ptB, GEN* ptC, /* cf hnfspec */
     613             :        GEN extramat,GEN extraC)
     614             : {
     615       91229 :   GEN matb, extratop, Cnew, permpro, B = *ptB, C = *ptC, dep = *ptdep;
     616             :   long i, lH, lB, li, lig, co, col, nlze;
     617             : 
     618       91229 :   if (lg(extramat) == 1) return H;
     619       91229 :   co = lg(C)-1;
     620       91229 :   lH = lg(H)-1;
     621       91229 :   lB = lg(B)-1;
     622       91229 :   li = lg(perm)-1;
     623       91229 :   lig = li - lB;
     624       91229 :   col = co - lB;
     625       91229 :   nlze = lig - lH;
     626             : 
     627             :  /*               col    co
     628             :   *       [ 0 |dep |     ]
     629             :   *  nlze [--------|  B  ]
     630             :   *       [ 0 | H  |     ]
     631             :   *       [--------|-----] lig
     632             :   *       [   0    | Id  ]
     633             :   *       [        |     ] li */
     634       91229 :   extratop = zm_to_ZM( rowslicepermute(extramat, perm, 1, lig) );
     635       91229 :   if (li != lig)
     636             :   { /* zero out bottom part, using the Id block */
     637       91061 :     GEN A = vecslice(C, col+1, co);
     638       91061 :     GEN c = rowslicepermute(extramat, perm, lig+1, li);
     639       91061 :     extraC   = gsub(extraC, typ(A)==t_MAT? RgM_zm_mul(A, c): RgV_zm_mul(A,c));
     640       91061 :     extratop = ZM_sub(extratop, ZM_zm_mul(B, c));
     641             :   }
     642             : 
     643       91227 :   extramat = shallowconcat(extratop, vconcat(dep, H));
     644       91229 :   Cnew     = shallowconcat(extraC, vecslice(C, col-lH+1, co));
     645       91229 :   if (DEBUGLEVEL>5) err_printf("    1st phase done\n");
     646       91229 :   permpro = ZM_rowrankprofile(extramat, &nlze);
     647       91229 :   extramat = rowpermute(extramat, permpro);
     648       91229 :   *ptB     = rowpermute(B,        permpro);
     649       91229 :   permpro = vecsmallpermute(perm, permpro);
     650      313587 :   for (i=1; i<=lig; i++) perm[i] = permpro[i]; /* perm o= permpro */
     651             : 
     652       91229 :   *ptdep  = rowslice(extramat, 1, nlze);
     653       91229 :   matb    = rowslice(extramat, nlze+1, lig);
     654       91229 :   if (DEBUGLEVEL>5) err_printf("    2nd phase done\n");
     655       91229 :   H = hnffinal(matb,perm,ptdep,ptB,&Cnew);
     656       91229 :   *ptC = shallowconcat(vecslice(C, 1, col-lH), Cnew);
     657       91229 :   return H;
     658             : }
     659             : 
     660             : GEN
     661           0 : hnfadd(GEN H, GEN perm, GEN* ptdep, GEN* ptB, GEN* ptC, /* cf hnfspec */
     662             :        GEN extramat,GEN extraC)
     663             : {
     664           0 :   pari_sp av = avma;
     665           0 :   H = hnfadd_i(H, perm, ptdep, ptB, ptC, ZM_to_zm(extramat), extraC);
     666           0 :   return gc_all(av, 4, &H, ptC, ptdep, ptB);
     667             : }
     668             : 
     669             : /* zero aj = Aij (!= 0)  using  ak = Aik (maybe 0), via linear combination of
     670             :  * A[j] and A[k] of determinant 1. If U != NULL, likewise update its columns */
     671             : static void
     672    39780275 : ZC_elem(GEN aj, GEN ak, GEN A, GEN U, long j, long k)
     673             : {
     674             :   GEN p1,u,v,d;
     675             : 
     676    39780275 :   if (!signe(ak)) {
     677      107182 :     swap(gel(A,j), gel(A,k));
     678      107182 :     if (U) swap(gel(U,j), gel(U,k));
     679    36079861 :     return;
     680             :   }
     681    39673093 :   d = bezout(aj,ak,&u,&v);
     682             :   /* frequent special case (u,v) = (1,0) or (0,1) */
     683    39681439 :   if (!signe(u))
     684             :   { /* ak | aj */
     685    19539170 :     p1 = diviiexact(aj,ak); togglesign(p1);
     686    19538909 :     ZC_lincomb1_inplace(gel(A,j), gel(A,k), p1);
     687    19547656 :     if (U)
     688     2266051 :       ZC_lincomb1_inplace(gel(U,j), gel(U,k), p1);
     689    19547102 :     return;
     690             :   }
     691    20142269 :   if (!signe(v))
     692             :   { /* aj | ak */
     693    16427565 :     p1 = diviiexact(ak,aj); togglesign(p1);
     694    16426049 :     ZC_lincomb1_inplace(gel(A,k), gel(A,j), p1);
     695    16425579 :     swap(gel(A,j), gel(A,k));
     696    16425579 :     if (U) {
     697      372038 :       ZC_lincomb1_inplace(gel(U,k), gel(U,j), p1);
     698      372036 :       swap(gel(U,j), gel(U,k));
     699             :     }
     700    16425577 :     return;
     701             :   }
     702             : 
     703     3714704 :   if (!is_pm1(d)) { aj = diviiexact(aj, d); ak = diviiexact(ak, d); }
     704     3714687 :   p1 = gel(A,k); aj = negi(aj); /* NOT togglesign */
     705     3715000 :   gel(A,k) = ZC_lincomb(u,v, gel(A,j),p1);
     706     3714517 :   gel(A,j) = ZC_lincomb(aj,ak, p1,gel(A,j));
     707     3714448 :   if (U)
     708             :   {
     709      750132 :     p1 = gel(U,k);
     710      750132 :     gel(U,k) = ZC_lincomb(u,v, gel(U,j),p1);
     711      750188 :     gel(U,j) = ZC_lincomb(aj,ak, p1,gel(U,j));
     712             :   }
     713             : }
     714             : 
     715             : INLINE int
     716        2590 : is_RgX(GEN a, long v) { return typ(a) == t_POL && varn(a)==v; }
     717             : /* set u,v such that au + bv = gcd(a,b), divide a,b by the gcd */
     718             : static GEN
     719         882 : gbezout_step(GEN *pa, GEN *pb, GEN *pu, GEN *pv, long vx)
     720             : {
     721         882 :   GEN a = *pa, b = *pb, d, l;
     722         882 :   if (gequal0(a))
     723             :   {
     724           0 :     *pa = gen_0; *pu = gen_0;
     725           0 :     *pb = gen_1; *pv = gen_1; return b;
     726             :   }
     727         882 :   a = is_RgX(a,vx)? RgX_renormalize(a): scalarpol(a, vx);
     728         882 :   b = is_RgX(b,vx)? RgX_renormalize(b): scalarpol(b, vx);
     729         882 :   d = RgX_extgcd(a,b, pu,pv);
     730         882 :   l = pollead(d,vx);
     731         882 :   d = RgX_Rg_div(d,l);
     732         882 :   *pu = RgX_Rg_div(*pu,l);
     733         882 :   *pv = RgX_Rg_div(*pv,l);
     734         882 :   if (degpol(d)) { a = RgX_div(a, d); b = RgX_div(b, d); }
     735         238 :   else if (typ(gel(d,2)) == t_REAL && lg(gel(d,2)) <= 3)
     736             : #if 1
     737             :   { /* possible accuracy problem */
     738           0 :     GEN D = RgX_gcd_simple(a,b);
     739           0 :     if (degpol(D)) {
     740           0 :       D = RgX_normalize(D);
     741           0 :       a = RgX_div(a, D);
     742           0 :       b = RgX_div(b, D);
     743           0 :       d = RgX_extgcd(a,b, pu,pv); /* retry now */
     744           0 :       d = RgX_mul(d, D);
     745             :     }
     746             :   }
     747             : #else
     748             :   { /* less stable */
     749             :     d = RgX_extgcd_simple(a,b, pu,pv);
     750             :     if (degpol(d)) { a = RgX_div(a, d); b = RgX_div(b, d); }
     751             :   }
     752             : #endif
     753         882 :   *pa = a;
     754         882 :   *pb = b; return d;
     755             : }
     756             : static GEN
     757     2698405 : col_mul(GEN x, GEN c)
     758             : {
     759     2698405 :   if (typ(x) == t_INT)
     760             :   {
     761     2696392 :     long s = signe(x);
     762     2696392 :     if (!s) return NULL;
     763     2044298 :     if (is_pm1(x)) return (s > 0)? c: RgC_neg(c);
     764             :   }
     765      308270 :   return RgC_Rg_mul(c, x);
     766             : }
     767             : static void
     768           0 : do_zero(GEN x)
     769             : {
     770           0 :   long i, lx = lg(x);
     771           0 :   for (i=1; i<lx; i++) gel(x,i) = gen_0;
     772           0 : }
     773             : 
     774             : /* (c1, c2) *= [u,-b; v,a] */
     775             : static void
     776      674612 : update(GEN u, GEN v, GEN a, GEN b, GEN *c1, GEN *c2)
     777             : {
     778             :   GEN p1,p2;
     779             : 
     780      674612 :   u = col_mul(u,*c1);
     781      674648 :   v = col_mul(v,*c2);
     782      674683 :   if (u) p1 = v? gadd(u,v): u;
     783        6505 :   else   p1 = v? v: NULL;
     784             : 
     785      674684 :   a = col_mul(a,*c2);
     786      674712 :   b = col_mul(gneg_i(b),*c1);
     787      674711 :   if (a) p2 = b? RgC_add(a,b): a;
     788           0 :   else   p2 = b? b: NULL;
     789             : 
     790      674655 :   if (!p1) do_zero(*c1); else *c1 = p1;
     791      674655 :   if (!p2) do_zero(*c2); else *c2 = p2;
     792      674655 : }
     793             : 
     794             : /* zero aj = Aij (!= 0)  using  ak = Aik (maybe 0), via linear combination of
     795             :  * A[j] and A[k] of determinant 1. If U != NULL, likewise update its columns */
     796             : static void
     797         511 : RgC_elem(GEN aj, GEN ak, GEN A, GEN V, long j, long k, long li, long vx)
     798             : {
     799         511 :   GEN u,v, d = gbezout_step(&aj, &ak, &u, &v, vx);
     800             :   long l;
     801             :   /* (A[,k], A[,j]) *= [v, -aj; u, ak ] */
     802        1778 :   for (l = 1; l < li; l++)
     803             :   {
     804        1267 :     GEN t = gadd(gmul(u,gcoeff(A,l,j)), gmul(v,gcoeff(A,l,k)));
     805        1267 :     gcoeff(A,l,j) = gsub(gmul(ak,gcoeff(A,l,j)), gmul(aj,gcoeff(A,l,k)));
     806        1267 :     gcoeff(A,l,k) = t;
     807             :   }
     808         511 :   gcoeff(A,li,j) = gen_0;
     809         511 :   gcoeff(A,li,k) = d;
     810         511 :   if (V) update(v,u,ak,aj,(GEN*)(V+k),(GEN*)(V+j));
     811         511 : }
     812             : 
     813             : /* reduce A[i,j] mod A[i,j0] for j=j0+1... via column operations */
     814             : static void
     815     4133144 : ZM_reduce(GEN A, GEN U, long i, long j0)
     816             : {
     817     4133144 :   long j, lA = lg(A);
     818     4133144 :   GEN d = gcoeff(A,i,j0);
     819     4133144 :   if (signe(d) < 0)
     820             :   {
     821       29585 :     ZV_neg_inplace(gel(A,j0));
     822       29585 :     if (U) ZV_togglesign(gel(U,j0));
     823       29585 :     d = gcoeff(A,i,j0);
     824             :   }
     825     8228207 :   for (j=j0+1; j<lA; j++)
     826             :   {
     827     4095019 :     GEN q = truedivii(gcoeff(A,i,j), d);
     828     4095040 :     if (!signe(q)) continue;
     829             : 
     830      257595 :     togglesign(q);
     831      257679 :     ZC_lincomb1_inplace(gel(A,j), gel(A,j0), q);
     832      257687 :     if (U) ZC_lincomb1_inplace(gel(U,j), gel(U,j0), q);
     833             :   }
     834     4133188 : }
     835             : 
     836             : /* normalize T as if it were a t_POL in variable v */
     837             : static GEN
     838         364 : normalize_as_RgX(GEN T, long v, GEN *pd)
     839             : {
     840             :   GEN d;
     841         364 :   if (!is_RgX(T,v)) { *pd = T; return gen_1; }
     842         336 :   d = leading_coeff(T);
     843         336 :   while (gequal0(d) || (typ(d) == t_REAL && lg(d) == 3
     844           0 :                         && gexpo(T) - expo(d) > (long)BITS_IN_LONG)) {
     845          14 :     T = normalizepol_lg(T, lg(T)-1);
     846          14 :     if (!signe(T)) { *pd = gen_1; return T; }
     847           0 :     d = leading_coeff(T);
     848             :   }
     849         322 :   if (degpol(T)) T = RgX_Rg_div(T,d); else { d = gel(T,2); T = gen_1; }
     850         322 :   *pd = d; return T;
     851             : }
     852             : /* reduce A[i,j] mod A[i,j0] for j=j0+1... via column operations */
     853             : static void
     854          84 : RgM_reduce(GEN A, GEN U, long i, long j0, long vx)
     855             : {
     856          84 :   long j, lA = lg(A);
     857          84 :   GEN d, T = normalize_as_RgX(gcoeff(A,i,j0), vx, &d);
     858          84 :   if (U && !gequal1(d)) gel(U,j0) = RgC_Rg_div(gel(U,j0), d);
     859          84 :   gcoeff(A,i,j0) = T;
     860             : 
     861         196 :   for (j=j0+1; j<lA; j++)
     862             :   {
     863         112 :     GEN t = gcoeff(A,i,j), q;
     864         112 :     if (gequal0(t)) continue;
     865          14 :     if (T == gen_1)
     866           0 :       q = t;
     867          14 :     else if (is_RgX(t,vx))
     868          14 :       q = RgX_div(t, T);
     869           0 :     else continue;
     870             : 
     871          14 :     if (gequal0(q)) continue;
     872           7 :     gel(A,j) = RgC_sub(gel(A,j), RgC_Rg_mul(gel(A,j0), q));
     873           7 :     if (U) gel(U,j) = RgC_sub(gel(U,j), RgC_Rg_mul(gel(U,j0), q));
     874             :   }
     875          84 : }
     876             : 
     877             : /* A,B square integral in upper HNF, of the same dimension > 0. Return Au
     878             :  * in Z^n (v in Z^n not computed), such that Au + Bv = [1, 0, ..., 0] */
     879             : GEN
     880      743844 : hnfmerge_get_1(GEN A, GEN B)
     881             : {
     882      743844 :   pari_sp av = avma;
     883      743844 :   long j, k, l = lg(A), lb;
     884      743844 :   GEN b, U = cgetg(l + 1, t_MAT), C = cgetg(l + 1, t_VEC);
     885             : 
     886      743917 :   b = gcoeff(B,1,1); lb = lgefint(b);
     887     1510814 :   for (j = 1; j < l; j++)
     888             :   {
     889             :     GEN t;
     890     1510806 :     long c = j+1;
     891     1510806 :     gel(U,j) = col_ei(l-1, j);
     892     1510977 :     gel(U,c) = zerocol(l-1); /* dummy */
     893     1510958 :     gel(C,j) = vecslice(gel(A,j), 1,j);
     894     1510983 :     gel(C,c) = vecslice(gel(B,j), 1,j);
     895     4195733 :     for (k = j; k > 0; k--)
     896             :     {
     897     2685282 :       t = gcoeff(C,k,c);
     898     2685282 :       if (gequal0(t)) continue;
     899     2479668 :       setlg(C[c], k+1);
     900     2479653 :       ZC_elem(t, gcoeff(C,k,k), C, U, c, k);
     901     2478712 :       if (lgefint(gcoeff(C,k,k)) > lb) gel(C,k) = FpC_red(gel(C,k), b);
     902     2478712 :       if (j > 4)
     903             :       {
     904       89324 :         GEN u = gel(U,k);
     905             :         long h;
     906     1337233 :         for (h=1; h<l; h++)
     907     1247909 :           if (lgefint(gel(u,h)) > lb) gel(u,h) = remii(gel(u,h), b);
     908             :       }
     909             :     }
     910     1510451 :     if (j == 1)
     911      743696 :       t = gcoeff(C,1,1);
     912             :     else
     913             :     {
     914             :       GEN u;
     915      766755 :       t = bezout(gcoeff(C,1,1), b, &u, NULL); /* >= 0 */
     916      767175 :       if (signe(u) && !equali1(u)) gel(U,1) = ZC_Z_mul(gel(U,1), u);
     917      767184 :       gcoeff(C,1,1) = t;
     918             :     }
     919     1510880 :     if (equali1(t)) break;
     920             :   }
     921      743813 :   if (j >= l) return NULL;
     922      743813 :   b = lcmii(gcoeff(A,1,1),b);
     923      743815 :   A = FpC_red(ZM_ZC_mul(A,gel(U,1)), b);
     924      743646 :   return gerepileupto(av, FpC_center(A, b, shifti(b,-1)));
     925             : }
     926             : 
     927             : /* remove the first r columns */
     928             : static void
     929      321529 : remove_0cols(long r, GEN *pA, GEN *pB, long remove)
     930             : {
     931      321529 :   GEN A = *pA, B = *pB;
     932      321529 :   long l = lg(A);
     933      321529 :   A += r; A[0] = evaltyp(t_MAT) | _evallg(l-r);
     934      321529 :   if (B && remove == 2) { B += r; B[0] = A[0]; }
     935      321529 :   *pA = A; *pB = B;
     936      321529 : }
     937             : 
     938             : /* Inefficient compared to hnfall. 'remove' = throw away lin.dep columns */
     939             : static GEN
     940       61618 : hnf_i(GEN A, int remove)
     941             : {
     942       61618 :   pari_sp av0 = avma, av;
     943             :   long s, n, m, j, k, li, def, ldef;
     944             : 
     945       61618 :   RgM_dimensions(A, &m, &n);
     946       61618 :   if (!n) return cgetg(1,t_MAT);
     947       61366 :   av = avma;
     948       61366 :   A = RgM_shallowcopy(A);
     949       61366 :   def = n; ldef = (m>n)? m-n: 0;
     950      164422 :   for (li=m; li>ldef; li--)
     951             :   {
     952      409287 :     for (j=def-1; j; j--)
     953             :     {
     954      306231 :       GEN a = gcoeff(A,li,j);
     955      306231 :       if (!signe(a)) continue;
     956             : 
     957             :       /* zero a = Aij  using  b = Aik */
     958      157577 :       k = (j==1)? def: j-1;
     959      157577 :       ZC_elem(a,gcoeff(A,li,k), A,NULL, j,k);
     960      157576 :       if (gc_needed(av,1))
     961             :       {
     962           0 :         if (DEBUGMEM>1) pari_warn(warnmem,"ZM_hnf[1]. li=%ld",li);
     963           0 :         A = gerepilecopy(av, A);
     964             :       }
     965             :     }
     966      103056 :     s = signe(gcoeff(A,li,def));
     967      103056 :     if (s)
     968             :     {
     969      101312 :       if (s < 0) ZV_neg_inplace(gel(A,def));
     970      101312 :       ZM_reduce(A, NULL, li,def);
     971      101312 :       def--;
     972             :     }
     973             :     else
     974        1744 :       if (ldef) ldef--;
     975      103056 :     if (gc_needed(av,1))
     976             :     {
     977           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"ZM_hnf[2]. li=%ld",li);
     978           0 :       A = gerepilecopy(av, A);
     979             :     }
     980             :   }
     981             :   /* rank A = n - def */
     982       61365 :   if (remove) { GEN B = NULL; remove_0cols(def, &A, &B, remove); }
     983       61366 :   return gerepileupto(av0, ZM_copy(A));
     984             : }
     985             : 
     986             : GEN
     987       88742 : ZM_hnf(GEN x) { return lg(x) > 8? ZM_hnfall(x, NULL, 1): hnf_i(x, 1); }
     988             : 
     989             : /* u*z[1..k] mod p, in place */
     990             : static void
     991     5308951 : FpV_Fp_mul_part_ip(GEN z, GEN u, GEN p, long k)
     992             : {
     993             :   long i;
     994     5308951 :   if (is_pm1(u)) {
     995      199517 :     if (signe(u) > 0) {
     996      589888 :       for (i = 1; i <= k; i++)
     997      444974 :         if (signe(gel(z,i))) gel(z,i) = modii(gel(z,i), p);
     998             :     } else {
     999      153081 :       for (i = 1; i <= k; i++)
    1000       98483 :         if (signe(gel(z,i))) gel(z,i) = modii(negi(gel(z,i)), p);
    1001             :     }
    1002             :   }
    1003             :   else {
    1004    19279968 :     for (i = 1; i <= k; i++)
    1005    14171379 :       if (signe(gel(z,i))) gel(z,i) = Fp_mul(u,gel(z,i), p);
    1006             :   }
    1007     5308101 : }
    1008             : static void
    1009    36041211 : FpV_red_part_ipvec(GEN z, GEN p, long k)
    1010             : {
    1011             :   long i;
    1012   109284543 :   for (i = 1; i <= k; i++) gel(z,i) = modii(gel(z,i), gel(p,i));
    1013    36037053 : }
    1014             : 
    1015             : /* return x * U, in echelon form (mod p^m), where (det(U),p) = 1.
    1016             :  * If early_abort is set, return NULL as soon as one pivot is 0 (mod p^m) */
    1017             : GEN
    1018      524903 : ZpM_echelon(GEN x, long early_abort, GEN p, GEN pm)
    1019             : {
    1020      524903 :   pari_sp av0 = avma, av;
    1021             :   long m, li, co, i, j, k, def, ldef;
    1022             : 
    1023      524903 :   co = lg(x); if (co == 1) return cgetg(1,t_MAT);
    1024      524903 :   li = lgcols(x);
    1025      524903 :   av = avma;
    1026      524903 :   x = RgM_shallowcopy(x);
    1027      524903 :   m = Z_pval(pm, p);
    1028             : 
    1029      524904 :   ldef = (li > co)? li - co: 0;
    1030     2502446 :   for (def = co-1,i = li-1; i > ldef; i--)
    1031             :   {
    1032     1978178 :     long vmin = LONG_MAX, kmin = 0;
    1033     1978178 :     GEN umin = gen_0, pvmin, q;
    1034     9653042 :     for (k = 1; k <= def; k++)
    1035             :     {
    1036     8138170 :       GEN u = gcoeff(x,i,k);
    1037             :       long v;
    1038     8138170 :       if (!signe(u)) continue;
    1039     3850409 :       v = Z_pvalrem(u, p, &u);
    1040     3850413 :       if (v >= m) gcoeff(x,i,k) = gen_0;
    1041     2717885 :       else if (v < vmin) {
    1042     1397585 :         vmin = v; kmin = k; umin = u;
    1043     1397585 :         if (!vmin) break;
    1044             :       }
    1045             :     }
    1046     1978182 :     if (!kmin)
    1047             :     {
    1048      923218 :       if (early_abort) return NULL;
    1049      922530 :       gcoeff(x,i,def) = gen_0;
    1050      922530 :       ldef--;
    1051      922530 :       if (ldef < 0) ldef = 0;
    1052      922530 :       continue;
    1053             :     }
    1054     1054964 :     if (kmin != def) swap(gel(x,def), gel(x,kmin));
    1055     1054964 :     q = vmin? powiu(p, m-vmin): pm;
    1056             :     /* pivot has valuation vmin */
    1057     1054962 :     umin = modii(umin, q);
    1058     1054950 :     if (!equali1(umin))
    1059      529538 :       FpV_Fp_mul_part_ip(gel(x,def), Fp_inv(umin,q), pm, i-1);
    1060     1054947 :     gcoeff(x, i, def) = pvmin = powiu(p, vmin);
    1061     7105839 :     for (j = def-1; j; j--)
    1062             :     { /* zero x[i, 1..def-1] using x[i,def] = pvmin */
    1063     6050827 :       GEN t, a = gcoeff(x,i,j) = modii(gcoeff(x,i,j), pm);
    1064     6050050 :       if (!signe(a)) continue;
    1065             : 
    1066     2835720 :       t = diviiexact(a, pvmin); togglesign(t);
    1067     2836523 :       ZC_lincomb1_inplace(gel(x,j), gel(x,def), t);
    1068     2836544 :       if (gc_needed(av,1))
    1069             :       {
    1070          14 :         if (DEBUGMEM>1) pari_warn(warnmem,"ZpM_echelon. i=%ld",i);
    1071          14 :         x = gerepilecopy(av, x); pvmin = gcoeff(x,i,def);
    1072             :       }
    1073             :     }
    1074     1055012 :     def--;
    1075             :   }
    1076      524268 :   if (co > li)
    1077             :   {
    1078           0 :     x += co - li;
    1079           0 :     x[0] = evaltyp(t_MAT) | _evallg(li);
    1080             :   }
    1081      524268 :   return gerepilecopy(av0, x);
    1082             : }
    1083             : GEN
    1084     4658057 : zlm_echelon(GEN x, long early_abort, ulong p, ulong pm)
    1085             : {
    1086     4658057 :   pari_sp av0 = avma;
    1087             :   long li, co, i, j, k, def, ldef;
    1088             :   ulong m;
    1089             : 
    1090     4658057 :   co = lg(x); if (co == 1) return cgetg(1,t_MAT);
    1091     4658057 :   li = lgcols(x);
    1092     4658054 :   x = Flm_copy(x);
    1093     4658154 :   m = u_lval(pm, p);
    1094             : 
    1095     4658177 :   ldef = (li > co)? li - co: 0;
    1096    20765123 :   for (def = co-1,i = li-1; i > ldef; i--)
    1097             :   {
    1098    16292491 :     long vmin = LONG_MAX, kmin = 0;
    1099    16292491 :     ulong umin = 0, pvmin, q;
    1100    49089269 :     for (k = 1; k <= def; k++)
    1101             :     {
    1102    41488600 :       ulong u = ucoeff(x,i,k);
    1103             :       long v;
    1104    41488600 :       if (!u) continue;
    1105    22667003 :       v = u_lvalrem(u, p, &u);
    1106    22666941 :       if (v >= (long) m) ucoeff(x,i,k) = 0;
    1107    22666941 :       else if (v < vmin) {
    1108    16337862 :         vmin = v; kmin = k; umin = u;
    1109    16337862 :         if (!vmin) break;
    1110             :       }
    1111             :     }
    1112    16292429 :     if (!kmin)
    1113             :     {
    1114     1014410 :       if (early_abort) return NULL;
    1115      829695 :       ucoeff(x,i,def) = 0;
    1116      829695 :       ldef--;
    1117      829695 :       if (ldef < 0) ldef = 0;
    1118      829695 :       continue;
    1119             :     }
    1120    15278019 :     if (kmin != def) swap(gel(x,def), gel(x,kmin));
    1121    15278019 :     q = vmin? upowuu(p, m-vmin): pm;
    1122             :     /* pivot has valuation vmin */
    1123    15278041 :     umin %= q;
    1124    15278041 :     if (umin != 1)
    1125     7733531 :       Flv_Fl_mul_part_inplace(gel(x,def), Fl_inv(umin,q), pm, i-1);
    1126    15278080 :     ucoeff(x, i, def) = pvmin = upowuu(p, vmin);
    1127    67515999 :     for (j = def-1; j; j--)
    1128             :     { /* zero x[i, 1..def-1] using x[i,def] = pvmin */
    1129    52238748 :       ulong t, a = ucoeff(x,i,j);
    1130    52238748 :       if (!a) continue;
    1131             : 
    1132    29492546 :       t = Fl_neg(a / pvmin, q);
    1133    29492526 :       Flc_lincomb1_inplace(gel(x,j), gel(x,def), t, pm);
    1134             :     }
    1135    15277251 :     def--;
    1136             :   }
    1137     4472632 :   if (co > li)
    1138             :   {
    1139           0 :     x += co - li;
    1140           0 :     x[0] = evaltyp(t_MAT) | _evallg(li);
    1141             :   }
    1142     4472632 :   return gerepilecopy(av0, x);
    1143             : }
    1144             : 
    1145             : static int
    1146      219982 : ZV_allequal(GEN v)
    1147             : {
    1148      219982 :   long i, l = lg(v);
    1149      219982 :   if (l > 1)
    1150             :   {
    1151      219982 :     GEN x = gel(v,1);
    1152      389417 :     for (i = 2; i < l; i++) if (!equalii(x,gel(v,i))) return 0;
    1153             :   }
    1154      199290 :   return 1;
    1155             : }
    1156             : /* compute optimal D for hnfmod: x upper triangular */
    1157             : static GEN
    1158     4963307 : optimal_D(GEN x, GEN D)
    1159             : {
    1160     4963307 :   long i, n = nbrows(x);
    1161     4963192 :   GEN C = shallowcopy(D);
    1162     4963609 :   gel(C,1) = gcoeff(x,1,1);
    1163     5384738 :   for (i = 2; i < n; i++)
    1164             :   {
    1165     2563106 :     GEN c = mulii(gel(C,i-1), gcoeff(x,i,i));
    1166     2562994 :     if (signe(c) < 0) togglesign(c);
    1167     2562994 :     if (cmpii(c, gel(D,i)) >= 0) break;
    1168      421129 :     gel(C,i) = c;
    1169             :   }
    1170     4963474 :   return C;
    1171             : }
    1172             : 
    1173             : /* D = multiple of det x (usually detint(x)) or vector of positive moduli
    1174             :  * (compute hnf(x | D))
    1175             :  * flag & hnf_MODID: reduce mod D * matid [ otherwise as above ].
    1176             :  * flag & hnf_PART: don't reduce once diagonal is known
    1177             :  * flag & hnf_CENTER: centermod HNF (2|x[i,j]| <] x[i,i]) */
    1178             : GEN
    1179     4965452 : ZM_hnfmodall_i(GEN x, GEN D, long flag)
    1180             : {
    1181             :   pari_sp av;
    1182     4965452 :   const long center = (flag & hnf_CENTER);
    1183             :   long moddiag, modsame, nli, li, co, i, j, k, def, ldef;
    1184             :   GEN u, LDM;
    1185             : 
    1186     4965452 :   co = lg(x);
    1187     4965452 :   if (co == 1)
    1188             :   {
    1189        1568 :     if (typ(D) == t_INT || lg(D) == 1) return cgetg(1,t_MAT);
    1190         105 :     x = diagonal_shallow(D); /* handle flags properly */
    1191         105 :     co = lg(x);
    1192             :   }
    1193     4963989 :   li = lgcols(x);
    1194     4963998 :   if (li == 1)
    1195             :   {
    1196         162 :     if (typ(D) != t_INT && lg(D) != li) pari_err_DIM("ZM_hnfmod");
    1197         162 :     return cgetg(1,t_MAT);
    1198             :   }
    1199     4963875 :   nli = li - 1;
    1200     4963875 :   modsame = typ(D)==t_INT;
    1201     4963875 :   if (!modsame)
    1202             :   {
    1203      219995 :     if (lg(D) != li) pari_err_DIM("ZM_hnfmod");
    1204      219981 :     if (ZV_allequal(D)) { modsame = 1; D = gel(D,1); }
    1205             :   }
    1206     4963862 :   moddiag = (flag & hnf_MODID) || !modsame;
    1207             :   /* modsame: triangularize mod fixed d*Id;
    1208             :    * moddiag: modulo diagonal matrix, else modulo multiple of determinant */
    1209             : 
    1210     4963862 :   if (modsame)
    1211             :   {
    1212     4943187 :     LDM = const_vecsmall(nli, 2*lgefint(D)-2);
    1213     4943958 :     D = const_vec(nli,D);
    1214             :   }
    1215             :   else
    1216             :   {
    1217       20675 :     LDM = cgetg(li, t_VECSMALL);
    1218       68067 :     for (i=1; i<li; i++) LDM[i] = lgefint(gel(D,i));
    1219             :   }
    1220     4964621 :   av = avma;
    1221     4964621 :   x = RgM_shallowcopy(x);
    1222             : 
    1223     4964866 :   ldef = 0;
    1224     4964866 :   if (li > co)
    1225             :   {
    1226      139345 :     ldef = li - co;
    1227      139345 :     if (!moddiag)
    1228           7 :       pari_err_DOMAIN("ZM_hnfmod","nb lines",">", strtoGENstr("nb columns"), x);
    1229             :   }
    1230    19569947 :   for (def = co-1,i = nli; i > ldef; i--,def--)
    1231             :   {
    1232    14598315 :     GEN d = gel(D,i);
    1233    14598315 :     long add_N = modsame;
    1234    63396786 :     for (j = 1; j < def; j++)
    1235             :     {
    1236    48791709 :       GEN p1, p2, b, a = gcoeff(x,i,j) = remii(gcoeff(x,i,j), d);
    1237    61220253 :       if (!signe(a)) continue;
    1238             : 
    1239    32816365 :       k = j+1;
    1240    32816365 :       b = gcoeff(x,i,k) = remii(gcoeff(x,i,k), d);
    1241    32819266 :       if (!signe(b)) { swap(gel(x,j), gel(x,k)); continue; }
    1242    20389462 :       if (add_N)
    1243             :       { /* ensure the moving pivot on row i divides d from now on */
    1244     7277657 :         add_N = 0;
    1245     7277657 :         if (!equali1(a))
    1246             :         { /* x[j] *= u; after this, a = x[i,j] | d */
    1247     5921092 :           GEN u = Fp_invgen(a, d, &a);
    1248             :           long t;
    1249     5922577 :           p1 = gel(x,j);
    1250    20911864 :           for (t = 1; t < i; t++) gel(p1,t) = mulii(gel(p1,t), u);
    1251     5921157 :           FpV_red_part_ipvec(p1, D, i-1);
    1252     5920525 :           gel(p1,i) = a;
    1253     5920525 :           if (2*lg(a) < lg(b))
    1254             :           { /* reduce x[i,k] mod x[i,j]: helps ZC_elem */
    1255       26591 :             GEN r, q = dvmdii(b, a, &r);
    1256       26591 :             togglesign(q);
    1257       26591 :             ZC_lincomb1_inplace_i(gel(x,k), gel(x,j), q, i-1);
    1258       26591 :             FpV_red_part_ipvec(gel(x,k), D, i-1);
    1259       26591 :             gcoeff(x,i,k) = b = r;
    1260             :           }
    1261             :         }
    1262             :       }
    1263    20388844 :       ZC_elem(a,b, x, NULL, j,k);
    1264    20394582 :       p1 = gel(x,j);
    1265    20394582 :       p2 = gel(x,k);
    1266             :       /* prevent coeffs explosion */
    1267   119427044 :       for (k = 1; k < i; k++)
    1268             :       {
    1269    99032461 :         if (lgefint(gel(p1,k)) > LDM[k]) gel(p1,k) = remii(gel(p1,k),gel(D,k));
    1270    99032462 :         if (lgefint(gel(p2,k)) > LDM[k]) gel(p2,k) = remii(gel(p2,k),gel(D,k));
    1271             :       }
    1272             :     }
    1273    14605077 :     if (gc_needed(av,2))
    1274             :     {
    1275          24 :       if (DEBUGMEM>1) pari_warn(warnmem,"ZM_hnfmod[1]. i=%ld",i);
    1276          24 :       x = gerepilecopy(av, x);
    1277             :     }
    1278    14605077 :     if (moddiag && !signe(gcoeff(x,i,def)))
    1279             :     { /* missing pivot on line i, insert column */
    1280      968180 :       GEN a = cgetg(co + 1, t_MAT);
    1281     5176069 :       for (k = 1; k <= def; k++) gel(a,k) = gel(x,k);
    1282      968182 :       gel(a,k++) = Rg_col_ei(gel(D,i), nli, i);
    1283     4153807 :       for (     ; k <= co;  k++) gel(a,k) = gel(x,k-1);
    1284      968187 :       ldef--; if (ldef < 0) ldef = 0;
    1285      968187 :       co++; def++; x = a;
    1286             :     }
    1287             :   }
    1288     4971632 :   if (co < li)
    1289             :   { /* implies moddiag, add missing diag(D) components */
    1290      107287 :     GEN a = cgetg(li+1, t_MAT);
    1291      258305 :     for (k = 1; k <= li-co; k++) gel(a,k) = Rg_col_ei(gel(D,k), nli, k);
    1292      308025 :     for (i = 1; i < co; i++) gel(a,k-1+i) = gel(x,i);
    1293      107286 :     gel(a,li) = zerocol(nli); x = a;
    1294             :   }
    1295             :   else
    1296             :   {
    1297     4864345 :     x += co - li;
    1298     4864345 :     x[0] = evaltyp(t_MAT) | _evallg(li); /* kill 0 columns */
    1299     4864345 :     if (moddiag) x = shallowconcat(x, zerocol(nli));
    1300             :   }
    1301     4964582 :   if (moddiag)
    1302             :   { /* x[li]: extra column, an accumulator discarded at the end */
    1303             :     GEN D2;
    1304     4799014 :     gcoeff(x,1,1) = gcdii(gcoeff(x,1,1), gel(D,1));
    1305     4797880 :     D2 = optimal_D(x,D);
    1306             :     /* add up missing diag(D) components */
    1307    19092590 :     for (i = nli; i > 0; i--)
    1308             :     {
    1309    14294131 :       gcoeff(x, i, li) = gel(D,i);
    1310    55115572 :       for (j = i; j > 0; j--)
    1311             :       {
    1312    40823048 :         GEN a = gcoeff(x, j, li);
    1313    40823048 :         if (!signe(a)) continue;
    1314    15060180 :         ZC_elem(a, gcoeff(x,j,j), x, NULL, li,j);
    1315    15059835 :         FpV_red_part_ipvec(gel(x,li), D, j-1);
    1316    15057325 :         FpV_red_part_ipvec(gel(x,j),  D, j-1);
    1317             :       }
    1318    14292524 :       if (gc_needed(av,1))
    1319             :       {
    1320          32 :         if (DEBUGMEM>1) pari_warn(warnmem,"ZM_hnfmod[2]. i=%ld", i);
    1321          32 :         gerepileall(av, 2, &x, &D2);
    1322             :       }
    1323             :     }
    1324     4798459 :     D = D2;
    1325             :   }
    1326             :   else
    1327             :   {
    1328      165568 :     GEN b = gel(D,1);
    1329      616781 :     for (i = nli; i > 0; i--)
    1330             :     {
    1331      451221 :       GEN d = bezout(gcoeff(x,i,i),b, &u,NULL);
    1332      451230 :       gcoeff(x,i,i) = d;
    1333      451230 :       FpV_Fp_mul_part_ip(gel(x,i), u, b, i-1);
    1334      451220 :       if (i > 1) b = diviiexact(b,d);
    1335             :     }
    1336      165560 :     D = optimal_D(x,D);
    1337             :   }
    1338     4964024 :   x[0] = evaltyp(t_MAT) | _evallg(li); /* kill 0 columns, discard accumulator */
    1339     4964024 :   if (flag & hnf_PART) return x;
    1340             : 
    1341    19711882 :   for (i = nli; i > 0; i--)
    1342             :   {
    1343    14751287 :     GEN diag = gcoeff(x,i,i);
    1344    14751287 :     if (signe(diag) < 0) { gel(x,i) = ZC_neg(gel(x,i)); diag = gcoeff(x,i,i); }
    1345    14751292 :     if (i != nli)
    1346    27177786 :       for (j = 1; j < i; j++) gcoeff(x,j,i) = remii(gcoeff(x,j,i), gel(D,j));
    1347    41926105 :     for (j = i+1; j < li; j++)
    1348             :     {
    1349    27176199 :       GEN b = gcoeff(x,i,j) = remii(gcoeff(x,i,j), gel(D,i));
    1350    27175004 :       GEN r, q = truedvmdii(b, diag, &r);
    1351             :       /* ensure -diag/2 <= r < diag/2 */
    1352    27174809 :       if (center && signe(r) && abscmpii(shifti(r,1),diag) >= 0)
    1353      620139 :       { r = subii(r,diag); q = addiu(q,1); }
    1354    27176496 :       if (!signe(q)) continue;
    1355    10883201 :       togglesign(q);
    1356    10883185 :       ZC_lincomb1_inplace_i(gel(x,j), gel(x,i), q, i-1);
    1357    10881576 :       gcoeff(x,i,j) = r;
    1358             :     }
    1359    14749906 :     if (gc_needed(av,1))
    1360             :     {
    1361          40 :       if (DEBUGMEM>1) pari_warn(warnmem,"ZM_hnfmod[3]. i=%ld", i);
    1362          40 :       gerepileall(av, 2, &x, &D);
    1363             :     }
    1364             :   }
    1365     4960595 :   return x;
    1366             : }
    1367             : GEN
    1368     4915881 : ZM_hnfmodall(GEN x, GEN dm, long flag)
    1369             : {
    1370     4915881 :   pari_sp av = avma;
    1371     4915881 :   return gerepilecopy(av, ZM_hnfmodall_i(x, dm, flag));
    1372             : }
    1373             : GEN
    1374      165561 : ZM_hnfmod(GEN x, GEN d) { return ZM_hnfmodall(x,d,0); }
    1375             : GEN
    1376     4340901 : ZM_hnfmodid(GEN x, GEN d)
    1377     4340901 : { return ZM_hnfmodall(x,d,hnf_MODID); }
    1378             : /* return the column echelon form of x with 1's as pivots,
    1379             :  * P contains the row indices containing the pivots in increasing order */
    1380             : static GEN
    1381     3144716 : FpM_echelon(GEN x, GEN *pP, GEN p)
    1382             : {
    1383             :   pari_sp av;
    1384             :   long iP, li, co, i, j, k, def, ldef;
    1385             :   GEN P;
    1386             : 
    1387     3144716 :   co = lg(x); if (co == 1) { *pP = cgetg(1,t_VECSMALL); return cgetg(1,t_MAT); }
    1388     3144716 :   li = lgcols(x);
    1389     3144702 :   iP = 1;
    1390     3144702 :   *pP = P = cgetg(li, t_VECSMALL);
    1391     3144702 :   av = avma;
    1392     3144702 :   x = FpM_red(x, p);
    1393             : 
    1394     3143974 :   ldef = (li > co)? li - co: 0;
    1395    11987126 :   for (def = co-1,i = li-1; i > ldef; i--)
    1396             :   {
    1397     8840555 :     GEN u = NULL;
    1398    13640182 :     for (k = def; k; k--)
    1399             :     {
    1400    10231828 :       u = gcoeff(x,i,k);
    1401    10231828 :       if (signe(u)) break;
    1402             :     }
    1403     8840555 :     if (!k)
    1404             :     {
    1405     3410452 :       if (--ldef < 0) ldef = 0;
    1406     3410452 :       continue;
    1407             :     }
    1408     5430103 :     P[iP++] = i;
    1409     5430103 :     if (k != def) swap(gel(x,def), gel(x,k));
    1410     5430103 :     if (!equali1(u))
    1411     4327971 :       FpV_Fp_mul_part_ip(gel(x,def), Fp_inv(u,p), p, i-1);
    1412     5432216 :     gcoeff(x, i, def) = gen_1;
    1413    17217626 :     for (j = def-1; j; j--)
    1414             :     { /* zero x[i, 1..def-1] using x[i,def] = 1*/
    1415    11785629 :       GEN xj = gel(x,j), u = gel(xj,i);
    1416    11785629 :       if (!signe(u)) continue;
    1417             : 
    1418     7455236 :       ZC_lincomb1_inplace(xj, gel(x,def), negi(u));
    1419    38509823 :       for (k = 1; k < i; k++) gel(xj,k) = modii(gel(xj,k), p);
    1420             :     }
    1421     5431997 :     if (gc_needed(av,2))
    1422             :     {
    1423           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"FpM_echelon. i=%ld",i);
    1424           0 :       x = gerepilecopy(av, x);
    1425             :     }
    1426     5432700 :     def--;
    1427             :   }
    1428             :   /* rank = iP - 1 */
    1429     3146571 :   setlg(P, iP); vecsmall_sort(P);
    1430     3144611 :   if (co > iP) x += co - iP;
    1431     3144611 :   x[0] = evaltyp(t_MAT) | _evallg(iP);
    1432     3144611 :   return x;
    1433             : }
    1434             : /* given x square of maximal rank with 1 or p on diagonal from hnfmodid
    1435             :  * (=> a column containing p has its other entries at 0 ), return the HNF */
    1436             : static GEN
    1437     3144843 : FpM_hnfend(pari_sp av, GEN x, GEN p)
    1438             : {
    1439     3144843 :   long i, l = lgcols(x);
    1440    11987549 :   for (i = l-1; i > 0; i--)
    1441             :   {
    1442     8844030 :     GEN diag = gcoeff(x,i,i);
    1443             :     long j;
    1444     8844030 :     if (is_pm1(diag))
    1445    11161870 :       for (j = i+1; j < l; j++)
    1446             :       {
    1447     5728377 :         GEN xj = gel(x,j), b = gel(xj,i);
    1448             :         long k;
    1449     5728377 :         if (!signe(b)) continue;
    1450     4027344 :         ZC_lincomb1_inplace(xj, gel(x,i), negi(b));
    1451    19564161 :         for (k=1; k<i; k++)
    1452    15536994 :           if (lgefint(gel(xj,k)) > 3) gel(xj,k) = remii(gel(xj,k), p);
    1453             :       }
    1454             :     else
    1455     9810523 :       for (j = i+1; j < l; j++) gcoeff(x,i,j) = modii(gcoeff(x,i,j), p);
    1456     8843485 :     if (gc_needed(av,2))
    1457             :     {
    1458           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"FpM_hnfend. i=%ld",i);
    1459           0 :       x = gerepilecopy(av, x);
    1460             :     }
    1461             :   }
    1462     3143519 :   return x;
    1463             : }
    1464             : GEN
    1465     3144712 : ZM_hnfmodprime(GEN x, GEN p)
    1466             : {
    1467     3144712 :   pari_sp av = avma;
    1468             :   GEN P, y;
    1469             :   long l, lP, i;
    1470     3144712 :   if (lg(x) == 1) return cgetg(1, t_MAT);
    1471     3144712 :   l = lgcols(x);
    1472     3144706 :   x = FpM_echelon(x, &P, p);
    1473     3144612 :   lP = lg(P); /* rank = lP-1 */
    1474     3144612 :   if (lP == l) { set_avma(av); return matid(l-1); }
    1475     3144612 :   y = scalarmat_shallow(p, l-1);
    1476     8578240 :   for (i = 1; i < lP; i++) gel(y,P[i]) = gel(x,i);
    1477     3144829 :   return gerepilecopy(av, FpM_hnfend(av,y,p));
    1478             : }
    1479             : 
    1480             : static GEN
    1481      408045 : allhnfmod(GEN x, GEN dm, int flag)
    1482             : {
    1483      408045 :   if (typ(x)!=t_MAT) pari_err_TYPE("allhnfmod",x);
    1484      408045 :   RgM_check_ZM(x, "allhnfmod");
    1485      408040 :   if (isintzero(dm)) return ZM_hnf(x);
    1486      408042 :   return ZM_hnfmodall(x, dm, flag);
    1487             : }
    1488             : GEN
    1489          14 : hnfmod(GEN x, GEN d)
    1490             : {
    1491          14 :   if (typ(d) != t_INT) pari_err_TYPE("mathnfmod",d);
    1492          14 :   return allhnfmod(x, d, 0);
    1493             : }
    1494             : GEN
    1495      408032 : hnfmodid(GEN x, GEN d)
    1496             : {
    1497      408032 :   switch(typ(d))
    1498             :   {
    1499      390846 :     case t_INT: break;
    1500       17186 :     case t_VEC: case t_COL:
    1501       17186 :       if (RgV_is_ZV(d)) break;
    1502           0 :     default: pari_err_TYPE("mathnfmodid",d);
    1503             :   }
    1504      408032 :   return allhnfmod(x, d, hnf_MODID);
    1505             : }
    1506             : 
    1507             : /* M a ZM in HNF. Normalize with *centered* residues */
    1508             : GEN
    1509       25160 : ZM_hnfcenter(GEN M)
    1510             : {
    1511       25160 :   long i, j, k, N = lg(M)-1;
    1512       25160 :   pari_sp av = avma;
    1513             : 
    1514       88799 :   for (j=N-1; j>0; j--) /* skip last line */
    1515             :   {
    1516       63640 :     GEN Mj = gel(M,j), a = gel(Mj,j);
    1517      213030 :     for (k = j+1; k <= N; k++)
    1518             :     {
    1519      149391 :       GEN Mk = gel(M,k), q = diviiround(gel(Mk,j), a);
    1520      149392 :       long s = signe(q);
    1521      149392 :       if (!s) continue;
    1522       68433 :       if (is_pm1(q))
    1523             :       {
    1524       32109 :         if (s < 0)
    1525        8604 :           for (i = 1; i <= j; i++) gel(Mk,i) = addii(gel(Mk,i), gel(Mj,i));
    1526             :         else
    1527       93827 :           for (i = 1; i <= j; i++) gel(Mk,i) = subii(gel(Mk,i), gel(Mj,i));
    1528             :       }
    1529             :       else
    1530      204311 :         for (i = 1; i <= j; i++) gel(Mk,i) = subii(gel(Mk,i), mulii(q,gel(Mj,i)));
    1531       68434 :       if (gc_needed(av,1))
    1532             :       {
    1533           0 :         if (DEBUGMEM>1) pari_warn(warnmem,"ZM_hnfcenter, j = %ld",j);
    1534           0 :         M = gerepilecopy(av, M);
    1535             :       }
    1536             :     }
    1537             :   }
    1538       25159 :   return M;
    1539             : }
    1540             : 
    1541             : /***********************************************************************/
    1542             : /*                                                                     */
    1543             : /*                 HNFLLL (Havas, Majewski, Mathews)                   */
    1544             : /*                                                                     */
    1545             : /***********************************************************************/
    1546             : 
    1547             : static void
    1548     2051644 : Minus(long j, GEN lambda)
    1549             : {
    1550     2051644 :   long k, n = lg(lambda);
    1551             : 
    1552    11758668 :   for (k=1  ; k<j; k++) togglesign_safe(&gcoeff(lambda,k,j));
    1553    17436265 :   for (k=j+1; k<n; k++) togglesign_safe(&gcoeff(lambda,j,k));
    1554     2051620 : }
    1555             : 
    1556             : /* index of first nonzero entry */
    1557             : static long
    1558   105856800 : findi(GEN M)
    1559             : {
    1560   105856800 :   long i, n = lg(M);
    1561  2612490978 :   for (i=1; i<n; i++)
    1562  2581062678 :     if (signe(gel(M,i))) return i;
    1563    31428300 :   return 0;
    1564             : }
    1565             : 
    1566             : static long
    1567   105846526 : findi_normalize(GEN Aj, GEN B, long j, GEN lambda)
    1568             : {
    1569   105846526 :   long r = findi(Aj);
    1570   105874481 :   if (r && signe(gel(Aj,r)) < 0)
    1571             :   {
    1572     2051475 :     ZV_togglesign(Aj); if (B) ZV_togglesign(gel(B,j));
    1573     2051632 :     Minus(j,lambda);
    1574             :   }
    1575   105878107 :   return r;
    1576             : }
    1577             : 
    1578             : static void
    1579     6447878 : subzi(GEN *a, GEN b)
    1580             : {
    1581     6447878 :   pari_sp av = avma;
    1582     6447878 :   b = subii(*a, b);
    1583     6447584 :   if (lgefint(b)<=lg(*a) && isonstack(*a)) { affii(b,*a); set_avma(av); }
    1584     1301170 :   else *a = b;
    1585     6448373 : }
    1586             : 
    1587             : static void
    1588   239246434 : addzi(GEN *a, GEN b)
    1589             : {
    1590   239246434 :   pari_sp av = avma;
    1591   239246434 :   b = addii(*a, b);
    1592   239328923 :   if (lgefint(b)<=lg(*a) && isonstack(*a)) { affii(b,*a); set_avma(av); }
    1593    17125025 :   else *a = b;
    1594   239673573 : }
    1595             : 
    1596             : /* x + y*z */
    1597             : static void
    1598   517651228 : addmulzii(GEN *x, GEN y, GEN z)
    1599             : {
    1600   517651228 :   long lx = lgefint(*x);
    1601             :   pari_sp av;
    1602             :   GEN t;
    1603   517651228 :   long ly = lgefint(y), lz;
    1604   517651228 :   if (ly == 2) return;
    1605   236058465 :   lz = lgefint(z);
    1606   236058465 :   av = avma; (void)new_chunk(lx+ly+lz); /*HACK*/
    1607   237620622 :   t = mulii(z, y);
    1608   236779849 :   set_avma(av); return addzi(x,t);
    1609             : }
    1610             : 
    1611             : static void
    1612    21833592 : ZC_lincomb1z_i(GEN X, GEN Y, GEN v, long n)
    1613             : {
    1614             :   long i;
    1615   500540030 :   for (i = n; i; i--) addmulzii(&gel(X,i), gel(Y,i), v);
    1616    21542974 : }
    1617             : /* X <- X + v Y (elementary col operation) */
    1618             : static void
    1619    21824475 : ZC_lincomb1z(GEN X, GEN Y, GEN v)
    1620             : {
    1621    21824475 :   if (lgefint(v) != 2) return ZC_lincomb1z_i(X, Y, v, lg(X)-1);
    1622             : }
    1623             : 
    1624             : static void
    1625    52908011 : reduce2(GEN A, GEN B, long k, long j, long *row0, long *row1, GEN lambda, GEN D)
    1626             : {
    1627             :   GEN q;
    1628             :   long i;
    1629             : 
    1630    52908011 :   *row0 = findi_normalize(gel(A,j), B,j,lambda);
    1631    52943382 :   *row1 = findi_normalize(gel(A,k), B,k,lambda);
    1632    52978433 :   if (*row0)
    1633    30721260 :     q = truedivii(gcoeff(A,*row0,k), gcoeff(A,*row0,j));
    1634             :   else
    1635             :   {
    1636    22257173 :     pari_sp btop = avma;
    1637    22257173 :     int c = abscmpii(shifti(gcoeff(lambda,j,k), 1), gel(D,j));
    1638    22216740 :     set_avma(btop);
    1639    22238947 :     if (c > 0)
    1640     7497879 :       q = diviiround(gcoeff(lambda,j,k), gel(D,j));
    1641             :     else
    1642    14741068 :       return;
    1643             :   }
    1644             : 
    1645    38211022 :   if (signe(q))
    1646             :   {
    1647    14665805 :     GEN Lk = gel(lambda,k), Lj = gel(lambda,j);
    1648    14665805 :     togglesign_safe(&q);
    1649    14668318 :     if (*row0) ZC_lincomb1z(gel(A,k),gel(A,j),q);
    1650    14668821 :     if (B) ZC_lincomb1z(gel(B,k),gel(B,j),q);
    1651    14636226 :     addmulzii(&gel(Lk,j), q, gel(D,j));
    1652    14631484 :     if (is_pm1(q))
    1653             :     {
    1654     5404243 :       if (signe(q) > 0)
    1655             :       {
    1656     5531236 :         for (i=1; i<j; i++)
    1657     3464372 :           if (signe(gel(Lj,i))) addzi(&gel(Lk,i), gel(Lj,i));
    1658             :       }
    1659             :       else
    1660             :       {
    1661    11202336 :         for (i=1; i<j; i++)
    1662     7865365 :           if (signe(gel(Lj,i))) subzi(&gel(Lk,i), gel(Lj,i));
    1663             :       }
    1664             :     }
    1665             :     else
    1666             :     {
    1667    41145452 :       for (i=1; i<j; i++)
    1668    31886045 :         if (signe(gel(Lj,i))) addmulzii(&gel(Lk,i), q, gel(Lj,i));
    1669             :     }
    1670             :   }
    1671             : }
    1672             : static void
    1673    47783412 : affii2_or_copy_gc(pari_sp av, GEN x, GEN *y, GEN x1, GEN *y1)
    1674             : {
    1675    47783412 :   long l = lg(*y), l1 = lg(*y1);
    1676    47783412 :   if(x==*y && x1==*y1) return;
    1677    47783412 :   if (lgefint(x) <= l && isonstack(*y) && lgefint(x1) <= l1 && isonstack(*y1))
    1678             :   {
    1679    46443514 :     affii(x,*y);
    1680    46476902 :     affii(x1,*y1);
    1681    46509589 :     set_avma(av);
    1682             :   }
    1683     1343764 :   else if (lgefint(x) <= l && isonstack(*y))
    1684             :   {
    1685      950339 :     affii(x,*y);
    1686      950369 :     *y1 = gerepilecopy(av, x1);
    1687             :   }
    1688      393435 :   else if (lgefint(x1) <= l1 && isonstack(*y1))
    1689             :   {
    1690      346709 :     affii(x1,*y1);
    1691      346709 :     *y = gerepilecopy(av, x);
    1692             :   }
    1693             :   else
    1694             :   {
    1695       46726 :     *y  = cgeti(2*lg(x));
    1696       52752 :     *y1 = cgeti(2*lg(x1));
    1697       52752 :     affii(x,*y);
    1698       52752 :     affii(x1,*y1);
    1699       52752 :     gerepileall(av,2,y,y1);
    1700             :   }
    1701             : }
    1702             : static void
    1703     8460602 : affii_or_copy_gc(pari_sp av, GEN x, GEN *y)
    1704             : {
    1705     8460602 :   long l = lg(*y);
    1706     8460602 :   if (lgefint(x) <= l && isonstack(*y))
    1707             :   {
    1708     6816760 :     affii(x,*y);
    1709     6819417 :     set_avma(av);
    1710             :   }
    1711             :   else
    1712     1643508 :     *y = gerepileuptoint(av, x);
    1713     8472836 : }
    1714             : static void
    1715     8472264 : hnfswap(GEN A, GEN B, long k, GEN lambda, GEN D)
    1716             : {
    1717             :   pari_sp av;
    1718     8472264 :   GEN t, p1, p2, Lk = gel(lambda,k);
    1719     8472264 :   long i,j,n = lg(A);
    1720             : 
    1721     8472264 :   swap(gel(A,k), gel(A,k-1));
    1722     8472264 :   if (B) swap(gel(B,k), gel(B,k-1));
    1723    46218251 :   for (j=k-2; j; j--) swap(gcoeff(lambda,j,k-1), gel(Lk,j));
    1724   106060599 :   for (i=k+1; i<n; i++)
    1725             :   {
    1726    97575273 :     pari_sp btop = avma;
    1727    97575273 :     GEN Li = gel(lambda,i);
    1728    97575273 :     if (signe(gel(Li,k-1))==0 && signe(gel(Li,k))==0) continue;
    1729    47825379 :     p1 = mulii(gel(Li,k-1), gel(D,k));
    1730    47804392 :     p2 = mulii(gel(Li,k), gel(Lk,k-1));
    1731    47796033 :     t = subii(p1,p2);
    1732    47796419 :     p1 = mulii(gel(Li,k), gel(D,k-2));
    1733    47792459 :     p2 = mulii(gel(Li,k-1), gel(Lk,k-1));
    1734    47793578 :     affii2_or_copy_gc(btop, diviiexact(addii(p1,p2), gel(D,k-1)),
    1735    47792118 :                       &gel(Li,k-1), diviiexact(t, gel(D,k-1)), &gel(Li,k));
    1736             :   }
    1737     8485326 :   av = avma;
    1738     8485326 :   p1 = mulii(gel(D,k-2), gel(D,k));
    1739     8467075 :   p2 = sqri(gel(Lk,k-1));
    1740     8462812 :   affii_or_copy_gc(av, diviiexact(addii(p1,p2), gel(D,k-1)), &gel(D,k-1));
    1741     8469825 : }
    1742             : 
    1743             : /* reverse row order in matrix A, IN PLACE */
    1744             : static GEN
    1745      462014 : reverse_rows(GEN A)
    1746             : {
    1747      462014 :   long i, j, h, n = lg(A);
    1748      462014 :   if (n == 1) return A;
    1749      462014 :   h = lgcols(A);
    1750     3856898 :   for (j=1; j<n; j++)
    1751             :   {
    1752     3394878 :     GEN c = gel(A,j);
    1753             :     /* start at (h-1) >>1 : if h = 2i even, no need to swap c[i] and itself */
    1754     9554410 :     for (i=(h-1)>>1; i; i--) swap(gel(c,i), gel(c,h-i));
    1755             :   }
    1756      462020 :   return A;
    1757             : }
    1758             : /* decide whether to swap */
    1759             : static int
    1760     4614223 : must_swap(long k, GEN lambda, GEN D)
    1761             : {
    1762     4614223 :   pari_sp av = avma;
    1763     4614223 :   GEN z = addii(mulii(gel(D,k-2),gel(D,k)), sqri(gcoeff(lambda,k-1,k)));
    1764     4612395 :   return gc_bool(av, cmpii(z, sqri(gel(D,k-1))) < 0);
    1765             : }
    1766             : 
    1767             : GEN
    1768      231018 : ZM_hnflll(GEN A, GEN *ptB, int remove)
    1769             : {
    1770      231018 :   pari_sp av = avma;
    1771             :   long n, k, kmax;
    1772             :   GEN B, lambda, D;
    1773             : 
    1774      231018 :   n = lg(A);
    1775      231018 :   A = reverse_rows(ZM_copy(A)); /* ZM_copy for in place findi_normalize() */
    1776      231017 :   B = ptB? matid(n-1): NULL;
    1777      231018 :   D = const_vec(n, gen_1) + 1;
    1778      231017 :   lambda = zeromatcopy(n-1,n-1);
    1779      215910 :   k = kmax = 2;
    1780    15738378 :   while (k < n)
    1781             :   {
    1782             :     long row0, row1;
    1783             :     int do_swap;
    1784    15507367 :     reduce2(A,B,k,k-1,&row0,&row1,lambda,D);
    1785    15500566 :     if      (row0) do_swap = (!row1 || row0 <= row1);
    1786     5557112 :     else if (row1) do_swap = 0;
    1787     4358392 :     else do_swap = must_swap(k,lambda,D);
    1788    15521305 :     if (do_swap)
    1789             :     {
    1790     8226494 :       hnfswap(A,B,k,lambda,D);
    1791     8230266 :       if (k > 2) k--;
    1792             :     }
    1793             :     else
    1794             :     {
    1795             :       long i;
    1796    44707056 :       for (i=k-2; i; i--)
    1797             :       {
    1798             :         long row0, row1;
    1799    37414854 :         reduce2(A,B,k,i,&row0,&row1,lambda,D);
    1800             :       }
    1801     7292202 :       if (++k > kmax) kmax = k;
    1802             :     }
    1803    15522468 :     if (gc_needed(av,3))
    1804             :     {
    1805           0 :       GEN b = D-1;
    1806           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"hnflll, kmax = %ld / %ld",kmax,n-1);
    1807           0 :       gerepileall(av, B? 4: 3, &A, &lambda, &b, &B);
    1808           0 :       if (gc_needed(av,1)) paristack_resize(0); /* avoid desperation GC */
    1809           0 :       D = b+1;
    1810             :     }
    1811             :   }
    1812             :   /* handle trivial case: return negative diag coefficient otherwise */
    1813      231011 :   if (n == 2) (void)findi_normalize(gel(A,1), B,1,lambda);
    1814      231011 :   A = reverse_rows(A);
    1815      231016 :   if (remove)
    1816             :   {
    1817             :     long i;
    1818        7546 :     for (i = 1; i < n; i++)
    1819        7322 :       if (!ZV_equal0(gel(A,i))) break;
    1820        2282 :     remove_0cols(i-1, &A, &B, remove);
    1821             :   }
    1822      231016 :   gerepileall(av, B? 2: 1, &A, &B);
    1823      231018 :   if (B) *ptB = B;
    1824      231018 :   return A;
    1825             : }
    1826             : 
    1827             : GEN
    1828           7 : hnflll(GEN x)
    1829             : {
    1830           7 :   GEN z = cgetg(3, t_VEC);
    1831           7 :   gel(z,1) = ZM_hnflll(x, &gel(z,2), 1);
    1832           7 :   return z;
    1833             : }
    1834             : 
    1835             : /* Variation on HNFLLL: Extended GCD */
    1836             : 
    1837             : static void
    1838      950510 : reduce1(GEN A, GEN B, long k, long j, GEN lambda, GEN D)
    1839             : {
    1840             :   GEN q;
    1841             :   long i;
    1842             : 
    1843      950510 :   if (signe(gel(A,j)))
    1844      172934 :     q = diviiround(gel(A,k),gel(A,j));
    1845      777576 :   else if (abscmpii(shifti(gcoeff(lambda,j,k), 1), gel(D,j)) > 0)
    1846       96207 :     q = diviiround(gcoeff(lambda,j,k), gel(D,j));
    1847             :   else
    1848      681369 :     return;
    1849             : 
    1850      269141 :   if (signe(q))
    1851             :   {
    1852      232549 :     GEN Lk = gel(lambda,k), Lj = gel(lambda,j);
    1853      232549 :     togglesign_safe(&q);
    1854      232549 :     gel(A,k) = addmulii(gel(A,k), q, gel(A,j));
    1855      232549 :     ZC_lincomb1_inplace(gel(B,k),gel(B,j),q);
    1856      232549 :     gel(Lk,j) = addmulii(gel(Lk,j), q, gel(D,j));
    1857      421761 :     for (i=1; i<j; i++)
    1858      189212 :       if (signe(gel(Lj,i))) gel(Lk,i) = addmulii(gel(Lk,i), q, gel(Lj,i));
    1859             :   }
    1860             : }
    1861             : 
    1862             : static GEN
    1863      109221 : ZV_gcdext_i(GEN A)
    1864             : {
    1865      109221 :   long k, n = lg(A);
    1866             :   GEN B, lambda, D;
    1867             : 
    1868      109221 :   if (n == 1) retmkvec2(gen_1, cgetg(1,t_MAT));
    1869      109221 :   A = leafcopy(A);
    1870      109221 :   B = matid(n-1);
    1871      109221 :   lambda = zeromatcopy(n-1,n-1);
    1872      109221 :   D = const_vec(n, gen_1) + 1;
    1873      109221 :   k = 2;
    1874      698177 :   while (k < n)
    1875             :   {
    1876             :     int do_swap;
    1877             : 
    1878      588956 :     reduce1(A,B,k,k-1,lambda,D);
    1879      588956 :     if    (signe(gel(A,k-1))) do_swap = 1;
    1880      416022 :     else if (signe(gel(A,k))) do_swap = 0;
    1881      223530 :     else do_swap = must_swap(k,lambda,D);
    1882      588956 :     if (do_swap)
    1883             :     {
    1884      241222 :       hnfswap(A,B,k,lambda,D);
    1885      241222 :       if (k > 2) k--;
    1886             :     }
    1887             :     else
    1888             :     {
    1889             :       long i;
    1890      709288 :       for (i=k-2; i; i--) reduce1(A,B,k,i,lambda,D);
    1891      347734 :       k++;
    1892             :     }
    1893             :   }
    1894      109221 :   if (signe(gel(A,n-1)) < 0)
    1895             :   {
    1896       13402 :     gel(A,n-1) = negi(gel(A,n-1));
    1897       13402 :     ZV_togglesign(gel(B,n-1));
    1898             :   }
    1899      109221 :   return mkvec2(gel(A,n-1), B);
    1900             : }
    1901             : GEN
    1902      109207 : ZV_extgcd(GEN A)
    1903             : {
    1904      109207 :   pari_sp av = avma;
    1905      109207 :   return gerepilecopy(av, ZV_gcdext_i(A));
    1906             : }
    1907             : /* as ZV_extgcd, transforming the gcd into a t_MAT, for mathnf0 */
    1908             : static GEN
    1909          21 : ZV_hnfgcdext(GEN A)
    1910             : {
    1911          21 :   pari_sp av = avma;
    1912             :   GEN z;
    1913          21 :   if (lg(A) == 1) retmkvec2(cgetg(1,t_MAT),cgetg(1,t_MAT));
    1914          14 :   z = ZV_gcdext_i(A);
    1915          14 :   gel(z,1) = mkmat(mkcol(gel(z,1)));
    1916          14 :   return gerepilecopy(av, z);
    1917             : }
    1918             : 
    1919             : /* HNF with permutation. */
    1920             : GEN
    1921         385 : ZM_hnfperm(GEN A, GEN *ptU, GEN *ptperm)
    1922             : {
    1923             :   GEN U, c, l, perm, d, p, q, b;
    1924         385 :   pari_sp av = avma, av1;
    1925             :   long r, t, i, j, j1, k, m, n;
    1926             : 
    1927         385 :   n = lg(A)-1;
    1928         385 :   if (!n)
    1929             :   {
    1930           7 :     if (ptU) *ptU = cgetg(1,t_MAT);
    1931           7 :     if (ptperm) *ptperm = cgetg(1,t_VEC);
    1932           7 :     return cgetg(1, t_MAT);
    1933             :   }
    1934         378 :   m = nbrows(A);
    1935         378 :   c = zero_zv(m);
    1936         378 :   l = zero_zv(n);
    1937         378 :   perm = cgetg(m+1, t_VECSMALL);
    1938         378 :   av1 = avma;
    1939         378 :   A = RgM_shallowcopy(A);
    1940         378 :   U = ptU? matid(n): NULL;
    1941             :   /* U base change matrix : A0*U = A all along */
    1942        1722 :   for (r=0, k=1; k <= n; k++)
    1943             :   {
    1944        3976 :     for (j=1; j<k; j++)
    1945             :     {
    1946        2632 :       if (!l[j]) continue;
    1947        2478 :       t=l[j]; b=gcoeff(A,t,k);
    1948        2478 :       if (!signe(b)) continue;
    1949             : 
    1950         532 :       ZC_elem(b,gcoeff(A,t,j), A,U,k,j);
    1951         532 :       d = gcoeff(A,t,j);
    1952         532 :       if (signe(d) < 0)
    1953             :       {
    1954           0 :         ZV_neg_inplace(gel(A,j));
    1955           0 :         if (U) ZV_togglesign(gel(U,j));
    1956           0 :         d = gcoeff(A,t,j);
    1957             :       }
    1958        1372 :       for (j1=1; j1<j; j1++)
    1959             :       {
    1960         840 :         if (!l[j1]) continue;
    1961         819 :         q = truedivii(gcoeff(A,t,j1),d);
    1962         819 :         if (!signe(q)) continue;
    1963             : 
    1964         329 :         togglesign(q);
    1965         329 :         ZC_lincomb1_inplace(gel(A,j1), gel(A,j), q);
    1966         329 :         if (U) ZC_lincomb1_inplace(gel(U,j1), gel(U,j), q);
    1967             :       }
    1968             :     }
    1969        4585 :     t = m; while (t && (c[t] || !signe(gcoeff(A,t,k)))) t--;
    1970        1344 :     if (t)
    1971             :     {
    1972        1225 :       p = gcoeff(A,t,k);
    1973       19971 :       for (i=t-1; i; i--)
    1974             :       {
    1975       18746 :         q = gcoeff(A,i,k);
    1976       18746 :         if (signe(q) && abscmpii(p,q) > 0) { p = q; t = i; }
    1977             :       }
    1978        1225 :       perm[++r] = l[k] = t; c[t] = k;
    1979        1225 :       if (signe(p) < 0)
    1980             :       {
    1981         142 :         ZV_neg_inplace(gel(A,k));
    1982         142 :         if (U) ZV_togglesign(gel(U,k));
    1983         142 :         p = gcoeff(A,t,k);
    1984             :       }
    1985             :       /* p > 0 */
    1986        3493 :       for (j=1; j<k; j++)
    1987             :       {
    1988        2268 :         if (!l[j]) continue;
    1989        2233 :         q = truedivii(gcoeff(A,t,j),p);
    1990        2233 :         if (!signe(q)) continue;
    1991             : 
    1992         238 :         togglesign(q);
    1993         238 :         ZC_lincomb1_inplace(gel(A,j), gel(A,k), q);
    1994         238 :         if (U) ZC_lincomb1_inplace(gel(U,j), gel(U,k), q);
    1995             :       }
    1996             :     }
    1997        1344 :     if (gc_needed(av1,1))
    1998             :     {
    1999           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"hnfperm, k=%ld",k);
    2000           0 :       gerepileall(av1, U? 2: 1, &A, &U);
    2001             :     }
    2002             :   }
    2003         378 :   if (r < m)
    2004             :   {
    2005        5131 :     for (i=1,k=r; i<=m; i++)
    2006        4816 :       if (!c[i]) perm[++k] = i;
    2007             :   }
    2008             : 
    2009             : /* We have A0*U=A, U in Gl(n,Z)
    2010             :  * basis for Im(A):  columns of A s.t l[j]>0 (r   cols)
    2011             :  * basis for Ker(A): columns of U s.t l[j]=0 (n-r cols) */
    2012         378 :   p = cgetg(r+1,t_MAT);
    2013        2814 :   for (i=1; i<=m/2; i++) lswap(perm[i], perm[m+1-i]);
    2014         378 :   if (U)
    2015             :   {
    2016          84 :     GEN u = cgetg(n+1,t_MAT);
    2017         378 :     for (t=1,k=r,j=1; j<=n; j++)
    2018         294 :       if (l[j])
    2019             :       {
    2020         182 :         u[k + n-r] = U[j];
    2021         182 :         gel(p,k--) = vecpermute(gel(A,j), perm);
    2022             :       }
    2023             :       else
    2024         112 :         u[t++] = U[j];
    2025          84 :     *ptU = u;
    2026          84 :     if (ptperm) *ptperm = perm;
    2027          84 :     gerepileall(av, ptperm? 3: 2, &p, ptU, ptperm);
    2028             :   }
    2029             :   else
    2030             :   {
    2031        1344 :     for (k=r,j=1; j<=n; j++)
    2032        1050 :       if (l[j]) gel(p,k--) = vecpermute(gel(A,j), perm);
    2033         294 :     if (ptperm) *ptperm = perm;
    2034         294 :     gerepileall(av, ptperm? 2: 1, &p, ptperm);
    2035             :   }
    2036         378 :   return p;
    2037             : }
    2038             : 
    2039             : GEN
    2040         238 : ZM_hnf_knapsack(GEN x)
    2041             : {
    2042         238 :   GEN t, perm, H = ZM_hnfperm(x,NULL,&perm);
    2043         238 :   long i,j, l = lg(H), h = lgcols(H);
    2044        3213 :   for (i=1; i<h; i++)
    2045             :   {
    2046        3066 :     int fl = 0;
    2047       16697 :     for (j=1; j<l; j++)
    2048             :     {
    2049       13722 :       t = gcoeff(H,i,j);
    2050       13722 :       if (signe(t))
    2051             :       {
    2052        3101 :         if (!is_pm1(t) || fl) return NULL;
    2053        3010 :         fl = 1;
    2054             :       }
    2055             :     }
    2056             :   }
    2057         147 :   return rowpermute(H, perm_inv(perm));
    2058             : }
    2059             : 
    2060             : GEN
    2061          14 : hnfperm(GEN A)
    2062             : {
    2063          14 :   GEN y = cgetg(4, t_VEC);
    2064          14 :   gel(y,1) = ZM_hnfperm(A, &gel(y,2), &gel(y,3));
    2065          14 :   return y;
    2066             : }
    2067             : 
    2068             : /* Hermite Normal Form, with base change matrix if ptB != NULL.
    2069             :  * If 'remove' = 1, remove 0 columns (do NOT update *ptB accordingly)
    2070             :  * If 'remove' = 2, remove 0 columns and update *ptB accordingly */
    2071             : GEN
    2072      725293 : ZM_hnfall_i(GEN A, GEN *ptB, long remove)
    2073             : {
    2074             :   pari_sp av;
    2075             :   long m, n, r, i, j, k, li;
    2076             :   GEN B, c, h, a;
    2077             : 
    2078      725293 :   RgM_dimensions(A, &m,&n);
    2079      725292 :   if (!n)
    2080             :   {
    2081           7 :     if (ptB) *ptB = cgetg(1,t_MAT);
    2082           7 :     return cgetg(1,t_MAT);
    2083             :   }
    2084      725285 :   c = zero_zv(m);
    2085      725277 :   h = const_vecsmall(n, m);
    2086      725275 :   av = avma;
    2087      725275 :   A = RgM_shallowcopy(A);
    2088      725287 :   B = ptB? matid(n): NULL;
    2089      725287 :   r = n+1;
    2090     2031511 :   for (li=m; li; li--)
    2091             :   {
    2092     2924693 :     for (j=1; j<r; j++)
    2093             :     {
    2094     3559540 :       for (i=h[j]; i>li; i--)
    2095             :       {
    2096      635327 :         a = gcoeff(A,i,j);
    2097      635327 :         k = c[i];
    2098             :         /* zero a = Aij  using  Aik */
    2099      635327 :         if (signe(a)) ZC_elem(a,gcoeff(A,i,k), A,B,j,k);
    2100      635267 :         ZM_reduce(A,B, i,k); /* ensure reduced entries even if a = 0 */
    2101             :       }
    2102     2924213 :       if (gc_needed(av,1) && (j & 0x7f) == 0)
    2103             :       {
    2104           0 :         if (DEBUGMEM>1)
    2105           0 :           pari_warn(warnmem,"ZM_hnfall[1], li = %ld, j = %ld", li, j);
    2106           0 :         gerepileall(av, B? 2: 1, &A, &B);
    2107             :       }
    2108     2924271 :       if (signe( gcoeff(A,li,j) )) break;
    2109     1618446 :       h[j] = li-1;
    2110             :     }
    2111     1306233 :     if (j == r) continue;
    2112     1305877 :     r--;
    2113     1305877 :     if (j < r) /* A[j] != 0 */
    2114             :     {
    2115      846680 :       swap(gel(A,j), gel(A,r));
    2116      846680 :       if (B) swap(gel(B,j), gel(B,r));
    2117      846680 :       h[j] = h[r]; h[r] = li; c[li] = r;
    2118             :     }
    2119     1305877 :     if (signe(gcoeff(A,li,r)) < 0)
    2120             :     {
    2121      302886 :       ZV_neg_inplace(gel(A,r));
    2122      302909 :       if (B) ZV_togglesign(gel(B,r));
    2123             :     }
    2124     1305898 :     ZM_reduce(A,B, li,r);
    2125     1305868 :     if (gc_needed(av,1))
    2126             :     {
    2127           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"ZM_hnfall[2], li = %ld", li);
    2128           0 :       gerepileall(av, B? 2: 1, &A, &B);
    2129             :     }
    2130             :   }
    2131             : 
    2132      725264 :   if (DEBUGLEVEL>5) err_printf("\nhnfall, final phase: ");
    2133      725277 :   r--; /* first r cols are in the image the n-r (independent) last ones */
    2134     2505785 :   for (j=1; j<=r; j++)
    2135             :   {
    2136     3871754 :     for (i=h[j]; i; i--)
    2137             :     {
    2138     2091245 :       a = gcoeff(A,i,j);
    2139     2091245 :       k = c[i];
    2140     2091245 :       if (signe(a)) ZC_elem(a,gcoeff(A,i,k), A,B, j,k);
    2141     2091208 :       ZM_reduce(A,B, i,k); /* ensure reduced entries, even if a = 0 */
    2142             :     }
    2143     1780509 :     if (gc_needed(av,1) && (j & 0x7f) == 0)
    2144             :     {
    2145           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"ZM_hnfall[3], j = %ld", j);
    2146           0 :       gerepileall(av, B? 2: 1, &A, &B);
    2147             :     }
    2148             :   }
    2149      725246 :   if (DEBUGLEVEL>5) err_printf("\n");
    2150      725246 :   if (remove) remove_0cols(r, &A, &B, remove);
    2151      725271 :   if (ptB) *ptB = B;
    2152      725271 :   return A;
    2153             : }
    2154             : GEN
    2155       28601 : ZM_hnfall(GEN A, GEN *ptB, long remove)
    2156             : {
    2157       28601 :   pari_sp av = avma;
    2158       28601 :   A = ZM_hnfall_i(A, ptB, remove);
    2159       28601 :   return gc_all(av, ptB? 2: 1, &A, ptB);
    2160             : }
    2161             : 
    2162             : GEN
    2163          14 : hnfall(GEN x)
    2164             : {
    2165          14 :   GEN z = cgetg(3, t_VEC);
    2166          14 :   gel(z,1) = ZM_hnfall(x, (GEN*)(z+2), 1);
    2167          14 :   return z;
    2168             : }
    2169             : GEN
    2170          14 : hnf(GEN x) { return mathnf0(x,0); }
    2171             : 
    2172             : /* C = A^(-1)t where A and C are integral, A is upper triangular, t t_INT */
    2173             : GEN
    2174      237211 : hnf_invscale(GEN A, GEN t)
    2175             : {
    2176      237211 :   long n = lg(A)-1, i,j,k;
    2177      237211 :   GEN m, c = cgetg(n+1,t_MAT);
    2178             : 
    2179      237210 :   if (!n) return c;
    2180      879877 :   for (k=1; k<=n; k++)
    2181             :   { /* cf hnf_divscale with B = id, thus b = e_k */
    2182      642702 :     GEN u = cgetg(n+1, t_COL);
    2183      642694 :     pari_sp av = avma;
    2184      642694 :     gel(c,k) = u;
    2185      642694 :     gel(u,n) = k == n? gerepileuptoint(av, diviiexact(t, gcoeff(A,n,n))): gen_0;
    2186     2175216 :     for (i=n-1; i>0; i--)
    2187             :     {
    2188     1532549 :       av = avma; m = i == k? t: gen_0;
    2189     5800018 :       for (j=i+1; j<=n; j++) m = subii(m, mulii(gcoeff(A,i,j),gel(u,j)));
    2190     1532357 :       gel(u,i) = gerepileuptoint(av, diviiexact(m, gcoeff(A,i,i)));
    2191             :     }
    2192             :   }
    2193      237175 :   return c;
    2194             : }
    2195             : 
    2196             : /* C = A^(-1)(tB) where A, B, C are integral, A is upper triangular, t t_INT */
    2197             : GEN
    2198      201726 : hnf_divscale(GEN A, GEN B, GEN t)
    2199             : {
    2200      201726 :   long n = lg(A)-1, i,j,k;
    2201      201726 :   GEN m, c = cgetg(n+1,t_MAT);
    2202             : 
    2203      201725 :   if (!n) return c;
    2204      889858 :   for (k=1; k<=n; k++)
    2205             :   {
    2206      688132 :     GEN u = cgetg(n+1, t_COL), b = gel(B,k);
    2207      688133 :     pari_sp av = avma;
    2208      688133 :     gel(c,k) = u; m = mulii(gel(b,n),t);
    2209      688126 :     gel(u,n) = gerepileuptoint(av, diviiexact(m, gcoeff(A,n,n)));
    2210     3929717 :     for (i=n-1; i>0; i--)
    2211             :     {
    2212     3241584 :       av = avma; m = mulii(gel(b,i),t);
    2213    30151501 :       for (j=i+1; j<=n; j++) m = subii(m, mulii(gcoeff(A,i,j),gel(u,j)));
    2214     3241554 :       gel(u,i) = gerepileuptoint(av, diviiexact(m, gcoeff(A,i,i)));
    2215             :     }
    2216             :   }
    2217      201726 :   return c;
    2218             : }
    2219             : 
    2220             : /* A, B integral upper HNF. A^(-1) B integral ? */
    2221             : int
    2222         133 : hnfdivide(GEN A, GEN B)
    2223             : {
    2224         133 :   pari_sp av = avma;
    2225         133 :   long n = lg(A)-1, i,j,k;
    2226             :   GEN u, b, m, r;
    2227             : 
    2228         133 :   if (!n) return 1;
    2229         133 :   if (lg(B)-1 != n) pari_err_DIM("hnfdivide");
    2230         133 :   u = cgetg(n+1, t_COL);
    2231         483 :   for (k=1; k<=n; k++)
    2232             :   {
    2233         364 :     b = gel(B,k);
    2234         364 :     m = gel(b,k);
    2235         364 :     gel(u,k) = dvmdii(m, gcoeff(A,k,k), &r);
    2236         364 :     if (r != gen_0) return gc_long(av, 0);
    2237         826 :     for (i=k-1; i>0; i--)
    2238             :     {
    2239         476 :       m = gel(b,i);
    2240        1337 :       for (j=i+1; j<=k; j++) m = subii(m, mulii(gcoeff(A,i,j),gel(u,j)));
    2241         476 :       m = dvmdii(m, gcoeff(A,i,i), &r);
    2242         476 :       if (r != gen_0) return gc_long(av, 0);
    2243         476 :       gel(u,i) = m;
    2244             :     }
    2245             :   }
    2246         119 :   return gc_long(av, 1);
    2247             : }
    2248             : 
    2249             : /* A upper HNF, b integral vector. Return A^(-1) b if integral,
    2250             :  * NULL otherwise. Assume #A[,1] = #b. */
    2251             : GEN
    2252      366153 : hnf_invimage(GEN A, GEN b)
    2253             : {
    2254      366153 :   pari_sp av = avma;
    2255      366153 :   long n = lg(A)-1, m, i, k;
    2256             :   GEN u, r;
    2257             : 
    2258      366153 :   if (!n) return lg(b)==1? cgetg(1,t_COL):NULL;
    2259      366111 :   m = nbrows(A); /* m >= n */
    2260      366110 :   u = cgetg(n+1, t_COL);
    2261      830300 :   for (i = n, k = m; k > 0; k--)
    2262             :   {
    2263      830299 :     pari_sp av2 = avma;
    2264             :     long j;
    2265      830299 :     GEN t = gel(b,k), Aki = gcoeff(A,k,i);
    2266      830299 :     if (typ(t) != t_INT) pari_err_TYPE("hnf_invimage",t);
    2267     2000200 :     for (j=i+1; j<=n; j++) t = subii(t, mulii(gcoeff(A,k,j),gel(u,j)));
    2268      830253 :     if (!signe(Aki))
    2269             :     {
    2270           0 :       if (signe(t)) return gc_NULL(av);
    2271           0 :       set_avma(av2); gel(u,i) = gen_0; continue;
    2272             :     }
    2273      830253 :     t = dvmdii(t, Aki, &r);
    2274      830144 :     if (r != gen_0) return gc_NULL(av);
    2275      607798 :     gel(u,i) = gerepileuptoint(av2, t);
    2276      607930 :     if (--i == 0) break;
    2277             :   }
    2278             :   /* If there is a solution, it must be u. Check remaining equations */
    2279      287507 :   for (; k > 0; k--)
    2280             :   {
    2281      143755 :     pari_sp av2 = avma;
    2282             :     long j;
    2283      143755 :     GEN t = gel(b,k);
    2284      143755 :     if (typ(t) != t_INT) pari_err_TYPE("hnf_invimage",t);
    2285      550086 :     for (j=1; j<=n; j++) t = subii(t, mulii(gcoeff(A,k,j),gel(u,j)));
    2286      143752 :     if (signe(t)) return gc_NULL(av);
    2287      143752 :     set_avma(av2);
    2288             :   }
    2289      143752 :   return u;
    2290             : }
    2291             : 
    2292             : /* A upper HNF, B integral matrix or column. Return A^(-1) B if integral,
    2293             :  * NULL otherwise */
    2294             : GEN
    2295      314903 : hnf_solve(GEN A, GEN B)
    2296             : {
    2297             :   pari_sp av;
    2298             :   long i, l;
    2299             :   GEN C;
    2300             : 
    2301      314903 :   if (typ(B) == t_COL) return hnf_invimage(A, B);
    2302      252187 :   av = avma; C = cgetg_copy(B, &l);
    2303      351471 :   for (i = 1; i < l; i++)
    2304             :   {
    2305      266467 :     GEN c = hnf_invimage(A, gel(B,i));
    2306      266435 :     if (!c) return gc_NULL(av);
    2307       99262 :     gel(C,i) = c;
    2308             :   }
    2309       85004 :   return C;
    2310             : }
    2311             : 
    2312             : /***************************************************************/
    2313             : /**                                                           **/
    2314             : /**               SMITH NORMAL FORM REDUCTION                 **/
    2315             : /**                                                           **/
    2316             : /***************************************************************/
    2317             : 
    2318             : static GEN
    2319           0 : trivsmith(long all)
    2320             : {
    2321             :   GEN z;
    2322           0 :   if (!all) return cgetg(1,t_VEC);
    2323           0 :   z=cgetg(4,t_VEC);
    2324           0 :   gel(z,1) = cgetg(1,t_MAT);
    2325           0 :   gel(z,2) = cgetg(1,t_MAT);
    2326           0 :   gel(z,3) = cgetg(1,t_MAT); return z;
    2327             : }
    2328             : 
    2329             : static void
    2330          48 : snf_pile1(pari_sp av, GEN *x, GEN *U)
    2331             : {
    2332             :   GEN *gptr[2];
    2333          48 :   int c = 1; gptr[0]=x;
    2334          48 :   if (*U) gptr[c++] = U;
    2335          48 :   gerepilemany(av,gptr,c);
    2336          48 : }
    2337             : static void
    2338      984751 : snf_pile(pari_sp av, GEN *x, GEN *U, GEN *V)
    2339             : {
    2340             :   GEN *gptr[3];
    2341      984751 :   int c = 1; gptr[0]=x;
    2342      984751 :   if (*U) gptr[c++] = U;
    2343      984751 :   if (*V) gptr[c++] = V;
    2344      984751 :   gerepilemany(av,gptr,c);
    2345      984831 : }
    2346             : 
    2347             : static GEN
    2348      702717 : bezout_step(GEN *pa, GEN *pb, GEN *pu, GEN *pv)
    2349             : {
    2350      702717 :   GEN a = *pa, b = *pb, d;
    2351      702717 :   if (absequalii(a,b))
    2352             :   {
    2353      418362 :     long sa = signe(a), sb = signe(b);
    2354      418362 :     *pv = gen_0;
    2355      418362 :     if (sb == sa) {
    2356      410694 :       *pa = *pb = gen_1;
    2357      410694 :       if (sa > 0) {*pu=gen_1; return a;} else {*pu=gen_m1; return absi(a);}
    2358             :     }
    2359        7668 :     if (sa > 0) { *pa = *pu = gen_1; *pb = gen_m1; return a; }
    2360        1162 :     *pa = *pu = gen_m1; *pb = gen_1; return b;
    2361             :   }
    2362      284364 :   d = bezout(a,b, pu,pv);
    2363      284379 :   *pa = diviiexact(a, d);
    2364      284374 :   *pb = diviiexact(b, d); return d;
    2365             : }
    2366             : 
    2367             : static int
    2368      520227 : negcmpii(void *E, GEN x, GEN y) { (void)E; return -cmpii(x,y); }
    2369             : 
    2370             : /* x square of maximal rank; does b = x[i,i] divide all entries in
    2371             :  * x[1..i-1, 1..i-1] ? If so, return 0; else the index of a problematic row */
    2372             : static long
    2373     1076989 : ZM_snf_no_divide(GEN x, long i)
    2374             : {
    2375     1076989 :   GEN b = gcoeff(x,i,i);
    2376             :   long j, k;
    2377             : 
    2378     1076989 :   if (is_pm1(b)) return 0;
    2379      912836 :   for (k = 1; k < i; k++)
    2380     2267851 :     for (j = 1; j < i; j++)
    2381     1791614 :       if (!dvdii(gcoeff(x,k,j),b)) return k;
    2382      283925 :   return 0;
    2383             : }
    2384             : 
    2385             : static void
    2386     1388353 : ZM_redpart(GEN x, GEN p, long I)
    2387             : {
    2388     1388353 :   long l = lgefint(p), i, j;
    2389     5914458 :   for (i = 1; i <= I; i++)
    2390    27223639 :     for (j = 1; j <= I; j++)
    2391             :     {
    2392    22697534 :       GEN c = gcoeff(x,i,j);
    2393    22697534 :       if (lgefint(c) > l) gcoeff(x,i,j) = remii(c, p);
    2394             :     }
    2395     1388353 : }
    2396             : static void
    2397      862913 : ZMrow_divexact_inplace(GEN M, long i, GEN c)
    2398             : {
    2399      862913 :   long j, l = lg(M);
    2400     3419943 :   for (j = 1; j < l; j++) gcoeff(M,i,j) = diviiexact(gcoeff(M,i,j), c);
    2401      863025 : }
    2402             : 
    2403             : /* Return the SNF D of matrix X. If ptU/ptV non-NULL set them to U/V
    2404             :  * to that D = UXV */
    2405             : GEN
    2406      722034 : ZM_snfall_i(GEN x, GEN *ptU, GEN *ptV, long flag)
    2407             : {
    2408      722034 :   pari_sp av0 = avma, av;
    2409      722034 :   const long return_vec = flag & 1;
    2410             :   long i, j, k, m0, m, n0, n;
    2411      722034 :   GEN u, v, U, V, V0, mdet, A = NULL, perm = NULL;
    2412             : 
    2413      722034 :   n0 = n = lg(x)-1;
    2414      722034 :   if (!n) {
    2415       41182 :     if (ptU) *ptU = cgetg(1,t_MAT);
    2416       41182 :     if (ptV) *ptV = cgetg(1,t_MAT);
    2417       41182 :     return cgetg(1, return_vec? t_VEC: t_MAT);
    2418             :   }
    2419      680852 :   m0 = m = nbrows(x);
    2420             : 
    2421      680861 :   U = V = V0 = NULL; /* U = TRANSPOSE of row transform matrix [act on columns]*/
    2422      680861 :   if (m == n && ZM_ishnf(x))
    2423             :   {
    2424      611225 :     mdet = ZM_det_triangular(x); /* != 0 */
    2425      611153 :     if (ptV) *ptV = matid(n);
    2426             :   }
    2427             :   else
    2428             :   {
    2429       69673 :     mdet = ZM_detmult(x);
    2430       69691 :     if (!signe(mdet))
    2431          77 :       x = ZM_hnfperm(x, ptV, ptU? &perm: NULL);
    2432             :     else
    2433             :     { /* m <= n */
    2434       69614 :       if (!ptV)
    2435        9520 :         x = ZM_hnfmod(x,mdet);
    2436       60094 :       else if (m == n)
    2437             :       {
    2438       60066 :         GEN H = ZM_hnfmod(x,mdet);
    2439       60066 :         *ptV = ZM_gauss(x,H);
    2440       60066 :         x = H;
    2441             :       }
    2442             :       else
    2443          28 :         x = ZM_hnfperm(x, ptV, ptU? &perm: NULL);
    2444       69614 :       mdet = ZM_det_triangular(x); /* > 0 */
    2445             :     }
    2446       69690 :     n = lg(x)-1; /* n independent columns */
    2447       69690 :     if (ptV)
    2448             :     {
    2449       60114 :       V = *ptV;
    2450       60114 :       if (n != n0)
    2451             :       {
    2452          28 :         V0 = vecslice(V, 1, n0 - n); /* kernel */
    2453          28 :         V  = vecslice(V, n0-n+1, n0);
    2454             :       }
    2455             :     }
    2456       69690 :     if (!signe(mdet))
    2457             :     {
    2458          77 :       if (n)
    2459             :       {
    2460          70 :         x = ZM_snfall_i(shallowtrans(x), ptV, ptU, return_vec); /* swap V,U */
    2461          70 :         if (!return_vec && n != m) x = shallowtrans(x);
    2462          70 :         if (ptV) V = ZM_mul(V, shallowtrans(*ptV));
    2463          70 :         if (ptU) U = *ptU; /* TRANSPOSE */
    2464             :       }
    2465             :       else /* 0 matrix */
    2466             :       {
    2467           7 :         x = cgetg(1,t_MAT);
    2468           7 :         if (ptV) V = cgetg(1, t_MAT);
    2469           7 :         if (ptU) U = matid(m);
    2470             :       }
    2471          77 :       goto THEEND;
    2472             :     }
    2473             :   }
    2474      680824 :   if (ptV || ptU) U = matid(n); /* we will compute V in terms of U */
    2475      680831 :   if (DEBUGLEVEL>7) err_printf("starting SNF loop");
    2476             : 
    2477             :   /* square, maximal rank n */
    2478      680831 :   A = x; x = shallowcopy(x); av = avma;
    2479     1605175 :   for (i = n; i > 1; i--)
    2480             :   {
    2481      924360 :     if (DEBUGLEVEL>7) err_printf("\ni = %ld: ",i);
    2482             :     for(;;)
    2483      703484 :     {
    2484     1627821 :       int c = 0;
    2485             :       GEN a, b;
    2486     4669389 :       for (j = i-1; j >= 1; j--)
    2487             :       {
    2488     3041715 :         b = gcoeff(x,i,j); if (!signe(b)) continue;
    2489      267051 :         a = gcoeff(x,i,i);
    2490      267051 :         ZC_elem(b, a, x,NULL, j,i);
    2491      267040 :         if (gc_needed(av,1))
    2492             :         {
    2493           0 :           if (DEBUGMEM>1) pari_warn(warnmem,"[1]: ZM_snfall i = %ld", i);
    2494           0 :           snf_pile1(av, &x,&U);
    2495             :         }
    2496             :       }
    2497     1627674 :       if (DEBUGLEVEL>7) err_printf("; ");
    2498     4669362 :       for (j=i-1; j>=1; j--)
    2499             :       {
    2500             :         GEN d;
    2501     3041671 :         b = gcoeff(x,j,i); if (!signe(b)) continue;
    2502      702716 :         a = gcoeff(x,i,i);
    2503      702716 :         d = bezout_step(&a, &b, &u, &v);
    2504     4190228 :         for (k = 1; k < i; k++)
    2505             :         {
    2506     3487596 :           GEN t = addii(mulii(u,gcoeff(x,i,k)),mulii(v,gcoeff(x,j,k)));
    2507     3487658 :           gcoeff(x,j,k) = subii(mulii(a,gcoeff(x,j,k)),
    2508     3487657 :                                 mulii(b,gcoeff(x,i,k)));
    2509     3487520 :           gcoeff(x,i,k) = t;
    2510             :         }
    2511      702632 :         gcoeff(x,j,i) = gen_0;
    2512      702632 :         gcoeff(x,i,i) = d;
    2513      702632 :         if (U) update(u,v,a,b,(GEN*)(U+i),(GEN*)(U+j));
    2514      702674 :         if (gc_needed(av,1))
    2515             :         {
    2516          48 :           if (DEBUGMEM>1) pari_warn(warnmem,"[2]: ZM_snfall, i = %ld", i);
    2517          48 :           snf_pile1(av, &x,&U);
    2518             :         }
    2519      702674 :         c = 1;
    2520             :       }
    2521     1627691 :       if (!c)
    2522             :       {
    2523     1076990 :         k = ZM_snf_no_divide(x, i);
    2524     1076950 :         if (!k) break;
    2525             : 
    2526             :         /* x[k,j] != 0 mod b */
    2527      582458 :         for (j = 1; j <= i; j++)
    2528      429859 :           gcoeff(x,i,j) = addii(gcoeff(x,i,j),gcoeff(x,k,j));
    2529      152599 :         if (U) gel(U,i) = gadd(gel(U,i),gel(U,k));
    2530             :       }
    2531      703297 :       ZM_redpart(x, mdet, i);
    2532      703446 :       if (U && (flag & 2)) ZM_redpart(U, mdet, n);
    2533      703480 :       if (gc_needed(av,1))
    2534             :       {
    2535           0 :         if (DEBUGMEM>1) pari_warn(warnmem,"[3]: ZM_snfall");
    2536           0 :         snf_pile1(av, &x,&U);
    2537             :       }
    2538             :     }
    2539             :   }
    2540      680815 :   if (DEBUGLEVEL>7) err_printf("\n");
    2541     2284978 :   for (k = n; k; k--)
    2542             :   {
    2543     1604605 :     GEN d = gcdii(gcoeff(x,k,k), mdet);
    2544     1604386 :     gcoeff(x,k,k) = d;
    2545     1604386 :     if (!is_pm1(d)) mdet = diviiexact(mdet,d);
    2546             :   }
    2547      680373 : THEEND:
    2548      680450 :   if (U) U = shallowtrans(U);
    2549      680690 :   if (ptV && A)
    2550             :   { /* U A V = D => D^(-1) U A = V^(-1) */
    2551      666524 :     long l = lg(x);
    2552      666524 :     GEN W = ZM_mul(U, A);
    2553     1529231 :     for (i = 1; i < l; i++)
    2554             :     {
    2555     1336856 :       GEN c = gcoeff(x,i,i);
    2556     1336856 :       if (is_pm1(c)) break; /* only 1 from now on */
    2557      862625 :       ZMrow_divexact_inplace(W, i, c);
    2558             :     }
    2559      666628 :     if (flag & 2)
    2560             :     {
    2561      642289 :       W = FpM_red(W, gcoeff(x,1,1));
    2562      642380 :       W = matinvmod(W, gcoeff(x,1,1));
    2563             :     }
    2564             :     else
    2565       24339 :       W = ZM_inv(W, NULL);
    2566      666557 :     V = V? ZM_mul(V, W): W;
    2567             :   }
    2568      680711 :   if (return_vec)
    2569             :   {
    2570      656662 :     long l = lg(x)-1;
    2571      656662 :     if (typ(x) == t_MAT) x = RgM_diagonal_shallow(x);
    2572      656650 :     if (m0 > l) x = shallowconcat(zerovec(m0-l), x);
    2573             :   }
    2574             : 
    2575      680699 :   if (V0)
    2576             :   { /* add kernel */
    2577          28 :     if (!return_vec) x = shallowconcat(zeromat(m,n0-n), x);
    2578          28 :     if (ptV) V = shallowconcat(V0, V);
    2579             :   }
    2580      680699 :   if (perm && U) U = vecpermute(U, perm_inv(perm));
    2581      680699 :   snf_pile(av0, &x,&U,&V);
    2582      680943 :   if (ptU) *ptU = U;
    2583      680943 :   if (ptV) *ptV = V;
    2584      680943 :   return x;
    2585             : }
    2586             : GEN
    2587       63972 : ZM_snfall(GEN x, GEN *U, GEN *V) { return ZM_snfall_i(x, U, V, 0); }
    2588             : GEN
    2589       10843 : ZM_snf(GEN x) { return ZM_snfall_i(x, NULL,NULL, 1); }
    2590             : GEN
    2591         385 : smith(GEN x) { return ZM_snfall_i(x, NULL,NULL, 1); }
    2592             : GEN
    2593          35 : smithall(GEN x)
    2594             : {
    2595          35 :   GEN z = cgetg(4, t_VEC);
    2596          35 :   gel(z,3) = ZM_snfall_i(x, (GEN*)(z+1),(GEN*)(z+2), 0);
    2597          35 :   return z;
    2598             : }
    2599             : 
    2600             : GEN
    2601      214281 : ZV_snfclean(GEN d)
    2602             : {
    2603      214281 :   long c, l = lg(d);
    2604      226650 :   for (c = 1; c < l; c++) { GEN t = gel(d,c); if (is_pm1(t)) break; }
    2605      214281 :   return c == l? d: vec_shorten(d, c-1);
    2606             : }
    2607             : void
    2608      947637 : ZM_snfclean(GEN d, GEN u, GEN v)
    2609             : {
    2610      947637 :   long i, c, l = lg(d);
    2611             : 
    2612      947637 :   if (typ(d) == t_VEC)
    2613     2432167 :     for (c=1; c<l; c++) { GEN t = gel(d,c); if (is_pm1(t)) break; }
    2614             :   else
    2615             :   {
    2616           0 :     for (c=1; c<l; c++) { GEN t = gcoeff(d,c,c); if (is_pm1(t)) break; }
    2617           0 :     if (c < l) for (i = 1; i < c; i++) setlg(gel(d,i), c);
    2618             :   }
    2619      947629 :   setlg(d, c);
    2620     2927746 :   if (u) for (i=1; i<l; i++) setlg(gel(u,i), c);
    2621      947613 :   if (v) setlg(v, c);
    2622      947613 : }
    2623             : 
    2624             : /* Assume z was computed by [g]smithall(). Remove the 1s on the diagonal */
    2625             : GEN
    2626         693 : smithclean(GEN z)
    2627             : {
    2628             :   long i, j, h, l, c, d;
    2629             :   GEN U, V, y, D, t;
    2630             : 
    2631         693 :   if (typ(z) != t_VEC) pari_err_TYPE("smithclean",z);
    2632         693 :   l = lg(z); if (l == 1) return cgetg(1,t_VEC);
    2633         686 :   U = gel(z,1);
    2634         686 :   if (l != 4 || typ(U) != t_MAT)
    2635             :   { /* assume z = vector of elementary divisors */
    2636        1827 :     for (c=1; c<l; c++)
    2637        1519 :       if (gequal1(gel(z,c))) break;
    2638         665 :     return gcopy_lg(z, c);
    2639             :   }
    2640          21 :   V = gel(z,2);
    2641          21 :   D = gel(z,3);
    2642          21 :   l = lg(D);
    2643          21 :   if (l == 1) return gcopy(z);
    2644          21 :   h = lgcols(D);
    2645          21 :   if (h > l)
    2646             :   { /* D = vconcat(zero matrix, diagonal matrix) */
    2647          21 :     for (c=1+h-l, d=1; c<h; c++,d++)
    2648          21 :       if (gequal1(gcoeff(D,c,d))) break;
    2649             :   }
    2650           7 :   else if (h < l)
    2651             :   { /* D = concat(zero matrix, diagonal matrix) */
    2652           7 :     for (c=1, d=1+l-h; d<l; c++,d++)
    2653           7 :       if (gequal1(gcoeff(D,c,d))) break;
    2654             :   }
    2655             :   else
    2656             :   { /* D diagonal */
    2657           0 :     for (c=1; c<l; c++)
    2658           0 :       if (gequal1(gcoeff(D,c,c))) break;
    2659           0 :     d = c;
    2660             :   }
    2661             :   /* U was (h-1)x(h-1), V was (l-1)x(l-1), D was (h-1)x(l-1) */
    2662          21 :   y = cgetg(4,t_VEC);
    2663             :   /* truncate U to (c-1) x (h-1) */
    2664          21 :   gel(y,1) = t = cgetg(h,t_MAT);
    2665          77 :   for (j=1; j<h; j++) gel(t,j) = gcopy_lg(gel(U,j), c);
    2666             :   /* truncate V to (l-1) x (d-1) */
    2667          21 :   gel(y,2) = gcopy_lg(V, d);
    2668          21 :   gel(y,3) = t = zeromatcopy(c-1, d-1);
    2669             :   /* truncate D to a (c-1) x (d-1) matrix */
    2670          21 :   if (d > 1)
    2671             :   {
    2672          14 :     if (h > l)
    2673             :     {
    2674          14 :       for (i=1+h-l, j=1; i<c; i++,j++)
    2675           7 :         gcoeff(t,i,j) = gcopy(gcoeff(D,i,j));
    2676             :     }
    2677           7 :     else if (h < l)
    2678             :     {
    2679           7 :       for (i=1, j=1+l-h; j<d; i++,j++)
    2680           0 :         gcoeff(t,i,j) = gcopy(gcoeff(D,i,j));
    2681             :     }
    2682             :     else
    2683             :     {
    2684           0 :       for (j=1; j<d; j++)
    2685           0 :         gcoeff(t,j,j) = gcopy(gcoeff(D,j,j));
    2686             :     }
    2687             :   }
    2688          21 :   return y;
    2689             : }
    2690             : 
    2691             : /* does b = x[i,i] divide all entries in x[1..i-1,1..i-1] ? If so, return 0;
    2692             :  * else return the index of a problematic row */
    2693             : static long
    2694         196 : gsnf_no_divide(GEN x, long i, long vx)
    2695             : {
    2696         196 :   GEN b = gcoeff(x,i,i);
    2697             :   long j, k;
    2698             : 
    2699         196 :   if (gequal0(b))
    2700             :   {
    2701          14 :     for (k = 1; k < i; k++)
    2702          14 :       for (j = 1; j < i; j++)
    2703          14 :         if (!gequal0(gcoeff(x,k,j))) return k;
    2704           0 :     return 0;
    2705             :   }
    2706             : 
    2707         182 :   if (!is_RgX(b,vx) || degpol(b)<=0) return 0;
    2708         217 :   for (k = 1; k < i; k++)
    2709         378 :     for (j = 1; j < i; j++)
    2710             :     {
    2711         266 :       GEN z = gcoeff(x,k,j), r;
    2712         266 :       if (!is_RgX(z,vx)) z = scalarpol(z, vx);
    2713         266 :       r = RgX_rem(z, b);
    2714         266 :       if (signe(r) && (! isinexactreal(r) ||
    2715           0 :              gexpo(r) > 16 + gexpo(b) - gprecision(r))
    2716          35 :          ) return k;
    2717             :     }
    2718          70 :   return 0;
    2719             : }
    2720             : 
    2721             : /* Hermite Normal Form, with base change matrix if ptB != NULL.
    2722             :  * If 'remove' = 1, remove 0 columns (do NOT update *ptB accordingly)
    2723             :  * If 'remove' = 2, remove 0 columns and update *ptB accordingly */
    2724             : GEN
    2725          42 : RgM_hnfall(GEN A, GEN *pB, long remove)
    2726             : {
    2727             :   pari_sp av;
    2728             :   long li, j, k, m, n, def, ldef;
    2729             :   GEN B;
    2730          42 :   long vx = gvar(A);
    2731             : 
    2732          42 :   n = lg(A)-1;
    2733          42 :   if (vx==NO_VARIABLE || !n)
    2734             :   {
    2735           0 :     RgM_check_ZM(A, "mathnf0");
    2736           0 :     return ZM_hnfall(A, pB, remove);
    2737             :   }
    2738          42 :   m = nbrows(A);
    2739          42 :   av = avma;
    2740          42 :   A = RgM_shallowcopy(A);
    2741          42 :   B = pB? matid(n): NULL;
    2742          42 :   def = n; ldef = (m>n)? m-n: 0;
    2743         126 :   for (li=m; li>ldef; li--)
    2744             :   {
    2745             :     GEN d, T;
    2746         714 :     for (j=def-1; j; j--)
    2747             :     {
    2748         630 :       GEN a = gcoeff(A,li,j);
    2749         630 :       if (gequal0(a)) continue;
    2750             : 
    2751         511 :       k = (j==1)? def: j-1;
    2752         511 :       RgC_elem(a,gcoeff(A,li,k), A,B, j,k, li, vx);
    2753             :     }
    2754          84 :     T = normalize_as_RgX(gcoeff(A,li,def), vx, &d);
    2755          84 :     if (gequal0(T))
    2756           0 :     { if (ldef) ldef--; }
    2757             :     else
    2758             :     {
    2759          84 :       if (!gequal1(d))
    2760             :       {
    2761           7 :         gel(A,def) = RgC_Rg_div(gel(A,def), d);
    2762           7 :         if (B) gel(B, def) = RgC_Rg_div(gel(B, def), d);
    2763             :       }
    2764          84 :       RgM_reduce(A, B, li, def, vx);
    2765          84 :       def--;
    2766             :     }
    2767          84 :     if (gc_needed(av,1))
    2768             :     {
    2769           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"ghnfall");
    2770           0 :       gerepileall(av, B? 2: 1, &A, &B);
    2771             :     }
    2772             :   }
    2773             :   /* rank A = n - def */
    2774          42 :   if (remove) remove_0cols(def, &A, &B, remove);
    2775          42 :   gerepileall(av, B? 2: 1, &A, &B);
    2776          42 :   if (B) *pB = B;
    2777          42 :   return A;
    2778             : }
    2779             : 
    2780             : static GEN
    2781          49 : RgXM_snf(GEN x,long all)
    2782             : {
    2783             :   pari_sp av;
    2784             :   long i, j, k, n;
    2785             :   GEN z, u, v, U, V;
    2786          49 :   long vx = gvar(x);
    2787          49 :   n = lg(x)-1; if (!n) return trivsmith(all);
    2788          49 :   if (vx==NO_VARIABLE) pari_err_TYPE("RgXM_snf",x);
    2789          49 :   if (lgcols(x) != n+1) pari_err_DIM("gsmithall");
    2790          49 :   av = avma;
    2791          49 :   x = RgM_shallowcopy(x);
    2792          49 :   if (all) { U = matid(n); V = matid(n); }
    2793         196 :   for (i=n; i>=2; i--)
    2794             :   {
    2795             :     for(;;)
    2796         168 :     {
    2797             :       GEN a, b, d;
    2798         315 :       int c = 0;
    2799        1092 :       for (j=i-1; j>=1; j--)
    2800             :       {
    2801         777 :         b = gcoeff(x,i,j); if (gequal0(b)) continue;
    2802         196 :         a = gcoeff(x,i,i);
    2803         196 :         d = gbezout_step(&b, &a, &v, &u, vx);
    2804         700 :         for (k = 1; k < i; k++)
    2805             :         {
    2806         504 :           GEN t = gadd(gmul(u,gcoeff(x,k,i)),gmul(v,gcoeff(x,k,j)));
    2807         504 :           gcoeff(x,k,j) = gsub(gmul(a,gcoeff(x,k,j)),gmul(b,gcoeff(x,k,i)));
    2808         504 :           gcoeff(x,k,i) = t;
    2809             :         }
    2810         196 :         gcoeff(x,i,j) = gen_0;
    2811         196 :         gcoeff(x,i,i) = d;
    2812         196 :         if (all) update(u,v,a,b,(GEN*)(V+i),(GEN*)(V+j));
    2813             :       }
    2814        1092 :       for (j=i-1; j>=1; j--)
    2815             :       {
    2816         777 :         b = gcoeff(x,j,i); if (gequal0(b)) continue;
    2817         175 :         a = gcoeff(x,i,i);
    2818         175 :         d = gbezout_step(&b, &a, &v, &u, vx);
    2819         651 :         for (k = 1; k < i; k++)
    2820             :         {
    2821         476 :           GEN t = gadd(gmul(u,gcoeff(x,i,k)),gmul(v,gcoeff(x,j,k)));
    2822         476 :           gcoeff(x,j,k) = gsub(gmul(a,gcoeff(x,j,k)),gmul(b,gcoeff(x,i,k)));
    2823         476 :           gcoeff(x,i,k) = t;
    2824             :         }
    2825         175 :         gcoeff(x,j,i) = gen_0;
    2826         175 :         gcoeff(x,i,i) = d;
    2827         175 :         if (all) update(u,v,a,b,(GEN*)(U+i),(GEN*)(U+j));
    2828         175 :         c = 1;
    2829             :       }
    2830         315 :       if (!c)
    2831             :       {
    2832         196 :         k = gsnf_no_divide(x, i, vx);
    2833         196 :         if (!k) break;
    2834             : 
    2835         245 :         for (j=1; j<=i; j++)
    2836         196 :           gcoeff(x,i,j) = gadd(gcoeff(x,i,j),gcoeff(x,k,j));
    2837          49 :         if (all) gel(U,i) = gadd(gel(U,i),gel(U,k));
    2838             :       }
    2839         168 :       if (gc_needed(av,1))
    2840             :       {
    2841           0 :         if (DEBUGMEM>1) pari_warn(warnmem,"gsmithall");
    2842           0 :         gerepileall(av, all? 3: 1, &x, &U, &V);
    2843             :       }
    2844             :     }
    2845             :   }
    2846         245 :   for (k=1; k<=n; k++)
    2847             :   {
    2848         196 :     GEN d, T = normalize_as_RgX(gcoeff(x,k,k), vx, &d);
    2849         196 :     if (gequal0(T)) continue;
    2850         182 :     if (all && !gequal1(d)) gel(V,k) = RgC_Rg_div(gel(V,k), d);
    2851         182 :     gcoeff(x,k,k) = T;
    2852             :   }
    2853          49 :   z = all? mkvec3(shallowtrans(U), V, x): RgM_diagonal_shallow(x);
    2854          49 :   return gerepilecopy(av, z);
    2855             : }
    2856             : 
    2857             : GEN
    2858         469 : matsnf0(GEN x,long flag)
    2859             : {
    2860         469 :   pari_sp av = avma;
    2861         469 :   if (flag > 7) pari_err_FLAG("matsnf");
    2862         469 :   if (typ(x) == t_VEC && flag & 4) return smithclean(x);
    2863         469 :   if (typ(x)!=t_MAT) pari_err_TYPE("matsnf",x);
    2864         469 :   if (RgM_is_ZM(x)) x = flag&1 ? smithall(x): smith(x);
    2865          49 :   else              x = RgXM_snf(x, flag&1);
    2866         469 :   if (flag & 4) x = gerepileupto(av, smithclean(x));
    2867         469 :   return x;
    2868             : }
    2869             : GEN
    2870           0 : gsmith(GEN x) { return RgXM_snf(x,0); }
    2871             : GEN
    2872           0 : gsmithall(GEN x) { return RgXM_snf(x,1); }
    2873             : 
    2874             : /* H is a relation matrix, either in HNF or a t_VEC (diagonal HNF) */
    2875             : static GEN
    2876      947644 : snf_group(GEN H, GEN D, GEN *newU, GEN *newUi)
    2877             : {
    2878             :   long i, j, l;
    2879             : 
    2880      947644 :   ZM_snfclean(D, newU? *newU: NULL, newUi? *newUi: NULL);
    2881      947613 :   l = lg(D);
    2882      947613 :   if (newU) {
    2883      825727 :     GEN U = *newU;
    2884     2129457 :     for (i = 1; i < l; i++)
    2885             :     {
    2886     1304268 :       GEN d = gel(D,i), d2 = shifti(d, 1);
    2887     5096484 :       for (j = 1; j < lg(U); j++)
    2888     3792754 :         gcoeff(U,i,j) = centermodii(gcoeff(U,i,j), d, d2);
    2889             :     }
    2890      825189 :     *newU = U;
    2891             :   }
    2892      947075 :   if (newUi && l > 1)
    2893             :   { /* UHV=D -> U^-1 = (HV)D^-1 -> U^-1 = H(VD^-1 mod 1) mod H */
    2894             :     /* Ui = ZM_inv(U, NULL); setlg(Ui, l); */
    2895      859971 :     GEN V = *newUi, Ui;
    2896      859971 :     int Hvec = (typ(H) == t_VEC);
    2897     2343367 :     for (i = 1; i < l; i++) gel(V,i) = FpC_red(gel(V,i), gel(D,i));
    2898      859943 :     if (!Hvec)
    2899             :     {
    2900      556506 :       if (ZM_isdiagonal(H)) { H = RgM_diagonal_shallow(H); Hvec = 1; }
    2901             :     }
    2902      860134 :     Ui = Hvec? ZM_diag_mul(H, V): ZM_mul(H, V);
    2903     2343212 :     for (i = 1; i < l; i++) gel(Ui,i) = ZC_Z_divexact(gel(Ui,i), gel(D,i));
    2904      859887 :     *newUi = Hvec? ZM_ZV_mod(Ui, H): ZM_hnfrem(Ui, H);
    2905             :   }
    2906      947163 :   return D;
    2907             : }
    2908             : /* H relation matrix among row of generators g in HNF.  Let URV = D its SNF,
    2909             :  * newU R newV = newD its clean SNF (no 1 in Dnew). Return the diagonal of
    2910             :  * newD, newU and newUi such that  1/U = (newUi, ?).
    2911             :  * Rationale: let (G,0) = g Ui be the new generators then
    2912             :  * 0 = G U R --> G D = 0,  g = G newU  and  G = g newUi */
    2913             : GEN
    2914      643689 : ZM_snf_group(GEN H, GEN *newU, GEN *newUi)
    2915             : {
    2916      643689 :   GEN D = ZM_snfall_i(H, newU, newUi, 1 + 2);
    2917      643773 :   return snf_group(H, D, newU, newUi);
    2918             : }
    2919             : 
    2920             : /* D a ZV: SNF for matdiagonal(D). Faster because we only ensure elementary
    2921             :  * divisors condition: d[n] | ... | d[1] and need not convert D to matrix form*/
    2922             : GEN
    2923      303890 : ZV_snfall(GEN D, GEN *pU, GEN *pV)
    2924             : {
    2925      303890 :   pari_sp av = avma;
    2926      303890 :   long j, n = lg(D)-1;
    2927      303890 :   GEN U = pU? matid(n): NULL;
    2928      303888 :   GEN V = pV? matid(n): NULL;
    2929             :   GEN p;
    2930             : 
    2931      303891 :   D = leafcopy(D);
    2932      966666 :   for (j = n; j > 0; j--)
    2933             :   {
    2934      662774 :     GEN b = gel(D,j);
    2935      662774 :     if (signe(b) < 0)
    2936             :     {
    2937           0 :       gel(D,j) = negi(b);
    2938           0 :       if (V) ZV_togglesign(gel(V,j));
    2939             :     }
    2940             :   }
    2941             :   /* entries are nonnegative integers */
    2942      303892 :   p = gen_indexsort(D, NULL, &negcmpii);
    2943      303889 :   D = vecpermute(D, p);
    2944      303888 :   if (U) U = vecpermute(U, p);
    2945      303888 :   if (V) V = vecpermute(V, p);
    2946             :   /* entries are sorted by decreasing value */
    2947      966622 :   for (j = n; j > 0; j--)
    2948             :   {
    2949      662729 :     GEN b = gel(D,j);
    2950             :     long i;
    2951     1184273 :     for (i = j-1; i > 0; i--)
    2952             :     { /* invariant: a >= b. If au+bv = d is a Bezout relation, A=a/d and B=b/d
    2953             :        * we have [B,-A;u,v]*diag(a,b)*[1-u*A,1; -u*A,1]] = diag(Ab, d) */
    2954      533733 :       GEN a = gel(D,i), u,v, d = bezout(a,b, &u,&v), A, Wi, Wj;
    2955      533725 :       if (equalii(d,b)) continue;
    2956       70388 :       A = diviiexact(a,d);
    2957       70388 :       if (V)
    2958             :       {
    2959       70332 :         GEN t = mulii(u,A);
    2960       70333 :         Wi = ZC_lincomb(subui(1,t), negi(t), gel(V,i), gel(V,j));
    2961       70330 :         Wj = ZC_add(gel(V,i), gel(V,j));
    2962       70330 :         gel(V,i) = Wi;
    2963       70330 :         gel(V,j) = Wj;
    2964             :       }
    2965       70386 :       if (U)
    2966             :       {
    2967       70386 :         GEN B = diviiexact(b,d);
    2968       70385 :         Wi = ZC_lincomb(B, negi(A), gel(U,i), gel(U,j));
    2969       70389 :         Wj = ZC_lincomb(u, v, gel(U,i), gel(U,j));
    2970       70387 :         gel(U,i) = Wi;
    2971       70387 :         gel(U,j) = Wj;
    2972             :       }
    2973       70387 :       gel(D,i) = mulii(A,b); /* lcm(a,b) */
    2974       70387 :       gel(D,j) = d; /* gcd(a,b) */
    2975       70387 :       b = gel(D,j); if (equali1(b)) break;
    2976             :     }
    2977             :   }
    2978      303893 :   snf_pile(av, &D,&U,&V);
    2979      303897 :   if (U) *pU = shallowtrans(U);
    2980      303892 :   if (V) *pV = V;
    2981      303892 :   return D;
    2982             : }
    2983             : GEN
    2984      303890 : ZV_snf_group(GEN d, GEN *newU, GEN *newUi)
    2985             : {
    2986      303890 :   GEN D = ZV_snfall(d, newU, newUi);
    2987      303892 :   return snf_group(d, D, newU, newUi);
    2988             : }
    2989             : 
    2990             : /* D a vector of elementary divisors. Truncate (setlg) to leave out trivial
    2991             :  * entries (= 1) */
    2992             : void
    2993           0 : ZV_snf_trunc(GEN D)
    2994             : {
    2995           0 :   long i, l = lg(D);
    2996           0 :   for (i = 1; i < l; i++)
    2997           0 :     if (is_pm1(gel(D,i))) { setlg(D,i); break; }
    2998           0 : }
    2999             : 
    3000             : long
    3001           0 : zv_snf_rank(GEN D, ulong p)
    3002             : {
    3003           0 :   long i, l = lg(D);
    3004           0 :   if (!p) return l - 1;
    3005           0 :   for (i = 1; i < l; i++)
    3006           0 :     if (D[i] % p) break;
    3007           0 :   return i - 1;
    3008             : }
    3009             : long
    3010          49 : ZV_snf_rank_u(GEN D, ulong p)
    3011             : {
    3012          49 :   long i, l = lg(D);
    3013          49 :   while (l > 1 && D[l-1] == 1) l--;
    3014          49 :   if (!p) return l - 1;
    3015          49 :   if (p == 2)
    3016             :   {
    3017          49 :     for (i = 1; i < l; i++)
    3018          42 :       if (mpodd(gel(D,i))) break;
    3019             :   }
    3020          35 :   else if (!(p & (p-1)))
    3021             :   { /* power of 2 */
    3022          14 :     long n = vals(p);
    3023          28 :     for (i = 1; i < l; i++)
    3024          28 :       if (umodi2n(gel(D,i), n)) break;
    3025             :   }
    3026             :   else
    3027             :   {
    3028          49 :     for (i = 1; i < l; i++)
    3029          42 :       if (umodiu(gel(D,i), p)) break;
    3030             :   }
    3031          49 :   return i - 1;
    3032             : }
    3033             : long
    3034          91 : ZV_snf_rank(GEN D, GEN p)
    3035             : {
    3036          91 :   long i, l = lg(D);
    3037          91 :   if (lgefint(p) == 3) return ZV_snf_rank_u(D, p[2]);
    3038          77 :   while (l > 1 && equali1(gel(D, l-1))) l--;
    3039          42 :   if (!signe(p)) return l - 1;
    3040          77 :   for (i = 1; i < l; i++)
    3041          70 :     if (!dvdii(gel(D,i), p)) break;
    3042          21 :   return i - 1;
    3043             : }
    3044             : long
    3045         154 : snfrank(GEN D, GEN p)
    3046             : {
    3047             :   long i, l;
    3048         154 :   if (typ(D) != t_VEC) pari_err_TYPE("snfrank", D);
    3049         154 :   if (!p) p = gen_0;
    3050         154 :   l = lg(D);
    3051         154 :   if (l == 4 && typ(gel(D,3)) == t_MAT)
    3052             :   { /* from matsnf(,1) */
    3053          14 :     pari_sp av = avma;
    3054             :     long z;
    3055             :     GEN v;
    3056          14 :     D = gel(D,3); l = lg(D);
    3057          14 :     if (l == 1) return 0;
    3058          14 :     z = lgcols(D) - l; /* missing columns of 0s */
    3059          14 :     if (z < 0) pari_err_TYPE("snfrank", D);
    3060          14 :     v = cgetg(l, t_VEC);
    3061          35 :     for (i = 1; i < l; i++) gel(v, i) = gcoeff(D, i + z, i);
    3062          14 :     return gc_long(av, z + snfrank(v, p));
    3063             :   }
    3064         140 :   switch(typ(p))
    3065             :   {
    3066          98 :     case t_INT:
    3067          98 :       if (RgV_is_ZV(D)) return ZV_snf_rank(D, p);
    3068           7 :       if (!signe(p)) break; /* allow p = 0 */
    3069           0 :       pari_err_TYPE("snfrank", D);
    3070          42 :     case t_POL: break;
    3071           0 :     default: pari_err_TYPE("snfrank", p);
    3072             :   }
    3073         175 :   while (l > 1 && gequal1(gel(D, l-1))) l--;
    3074          49 :   if (gequal0(p)) return l - 1;
    3075         112 :   for (i = 1; i < l; i++)
    3076          91 :     if (!gdvd(gel(D,i), p)) break;
    3077          42 :   return i - 1;
    3078             : }

Generated by: LCOV version 1.16