Code coverage tests

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

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

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

LCOV - code coverage report
Current view: top level - basemath - hnf_snf.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.10.0 lcov report (development 21348-d75f58f) Lines: 1493 1639 91.1 %
Date: 2017-11-20 06:21:05 Functions: 84 89 94.4 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2000  The PARI group.
       2             : 
       3             : This file is part of the PARI/GP package.
       4             : 
       5             : PARI/GP is free software; you can redistribute it and/or modify it under the
       6             : terms of the GNU General Public License as published by the Free Software
       7             : Foundation. It is distributed in the hope that it will be useful, but WITHOUT
       8             : ANY WARRANTY WHATSOEVER.
       9             : 
      10             : Check the License for details. You should have received a copy of it, along
      11             : with the package; see the file 'COPYING'. If not, write to the Free Software
      12             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      13             : 
      14             : #include "pari.h"
      15             : #include "paripriv.h"
      16             : /**************************************************************/
      17             : /**                                                          **/
      18             : /**                HERMITE NORMAL FORM REDUCTION             **/
      19             : /**                                                          **/
      20             : /**************************************************************/
      21             : static GEN ZV_hnfgcdext(GEN A);
      22             : static GEN
      23          14 : hnfallgen(GEN x)
      24             : {
      25          14 :   GEN z = cgetg(3, t_VEC);
      26          14 :   gel(z,1) = RgM_hnfall(x, (GEN*)(z+2), 1);
      27          14 :   return z;
      28             : }
      29             : GEN
      30         168 : mathnf0(GEN x, long flag)
      31             : {
      32         168 :   switch(typ(x))
      33             :   {
      34             :     case t_VEC:
      35          70 :       if (RgV_is_ZV(x))
      36          42 :         switch (flag)
      37             :         {
      38             :           case 0:
      39          14 :             if (lg(x) == 1) return cgetg(1, t_MAT);
      40           7 :             retmkmat(mkcol(ZV_content(x)));
      41             :           case 1:
      42             :           case 4:
      43          21 :             return ZV_hnfgcdext(x);
      44             :         }
      45          35 :       x = gtomat(x); break;
      46          98 :     case t_MAT: break;
      47           0 :     default: pari_err_TYPE("mathnf0",x);
      48             :   }
      49             : 
      50         133 :   switch(flag)
      51             :   {
      52          84 :     case 0: case 2: return RgM_is_ZM(x)? ZM_hnf(x): RgM_hnfall(x,NULL,1);
      53          28 :     case 1: case 3: return RgM_is_ZM(x)? hnfall(x): hnfallgen(x);
      54           7 :     case 4: RgM_check_ZM(x, "mathnf0"); return hnflll(x);
      55          14 :     case 5: RgM_check_ZM(x, "mathnf0"); return hnfperm(x);
      56           0 :     default: pari_err_FLAG("mathnf");
      57             :   }
      58             :   return NULL; /* LCOV_EXCL_LINE */
      59             : }
      60             : 
      61             : /*******************************************************************/
      62             : /*                                                                 */
      63             : /*                SPECIAL HNF (FOR INTERNAL USE !!!)               */
      64             : /*                                                                 */
      65             : /*******************************************************************/
      66             : static int
      67     3104417 : count(GEN mat, long row, long len, long *firstnonzero)
      68             : {
      69     3104417 :   long j, n = 0;
      70             : 
      71   213794778 :   for (j=1; j<=len; j++)
      72             :   {
      73   211372590 :     long p = mael(mat,j,row);
      74   211372590 :     if (p)
      75             :     {
      76     8350390 :       if (labs(p)!=1) return -1;
      77     7668161 :       n++; *firstnonzero=j;
      78             :     }
      79             :   }
      80     2422188 :   return n;
      81             : }
      82             : 
      83             : static int
      84      103224 : count2(GEN mat, long row, long len)
      85             : {
      86             :   long j;
      87      944652 :   for (j=len; j; j--)
      88      911275 :     if (labs(mael(mat,j,row)) == 1) return j;
      89       33377 :   return 0;
      90             : }
      91             : 
      92             : static GEN
      93       41858 : hnffinal(GEN matgen,GEN perm,GEN* ptdep,GEN* ptB,GEN* ptC)
      94             : {
      95             :   GEN p1,p2,U,H,Hnew,Bnew,Cnew,diagH1;
      96       41858 :   GEN B = *ptB, C = *ptC, dep = *ptdep, depnew;
      97             :   pari_sp av;
      98             :   long i,j,k,s,i1,j1,zc;
      99       41858 :   long co = lg(C);
     100       41858 :   long col = lg(matgen)-1;
     101             :   long lnz, nlze, lig;
     102             : 
     103       41858 :   if (col == 0) return matgen;
     104       41501 :   lnz = nbrows(matgen); /* was called lnz-1 - nr in hnfspec */
     105       41501 :   nlze = nbrows(dep);
     106       41501 :   lig = nlze + lnz;
     107             :   /* H: lnz x lnz [disregarding initial 0 cols], U: col x col */
     108       41501 :   H = ZM_hnflll(matgen, &U, 0);
     109       41501 :   H += (lg(H)-1 - lnz); H[0] = evaltyp(t_MAT) | evallg(lnz+1);
     110             :   /* Only keep the part above the H (above the 0s is 0 since the dep rows
     111             :    * are dependent from the ones in matgen) */
     112       41501 :   zc = col - lnz; /* # of 0 columns, correspond to units */
     113       41501 :   if (nlze) { dep = ZM_mul(dep,U); dep += zc; }
     114             : 
     115       41501 :   diagH1 = new_chunk(lnz+1); /* diagH1[i] = 0 iff H[i,i] != 1 (set later) */
     116             : 
     117       41501 :   av = avma;
     118       41501 :   Cnew = cgetg(co, typ(C));
     119       41501 :   setlg(C, col+1); p1 = gmul(C,U);
     120       41501 :   for (j=1; j<=col; j++) gel(Cnew,j) = gel(p1,j);
     121       41501 :   for (   ; j<co ; j++)  gel(Cnew,j) = gel(C,j);
     122             : 
     123             :   /* Clean up B using new H */
     124      197551 :   for (s=0,i=lnz; i; i--)
     125             :   {
     126      156050 :     GEN Di = gel(dep,i), Hi = gel(H,i);
     127      156050 :     GEN h = gel(Hi,i); /* H[i,i] */
     128      156050 :     if ( (diagH1[i] = is_pm1(h)) ) { h = NULL; s++; }
     129     7536529 :     for (j=col+1; j<co; j++)
     130             :     {
     131     7380479 :       GEN z = gel(B,j-col);
     132     7380479 :       p1 = gel(z,i+nlze);
     133     7380479 :       if (h) p1 = truedivii(p1,h);
     134     7380479 :       if (!signe(p1)) continue;
     135     4390136 :       for (k=1; k<=nlze; k++) gel(z,k) = subii(gel(z,k), mulii(p1, gel(Di,k)));
     136     4390136 :       for (   ; k<=lig;  k++) gel(z,k) = subii(gel(z,k), mulii(p1, gel(Hi,k-nlze)));
     137     4390136 :       gel(Cnew,j) = gsub(gel(Cnew,j), gmul(p1, gel(Cnew,i+zc)));
     138             :     }
     139      156050 :     if (gc_needed(av,3))
     140             :     {
     141          47 :       if(DEBUGMEM>1) pari_warn(warnmem,"hnffinal, i = %ld",i);
     142          47 :       gerepileall(av, 2, &Cnew, &B);
     143             :     }
     144             :   }
     145       41501 :   p1 = cgetg(lnz+1,t_VEC); p2 = perm + nlze;
     146      197551 :   for (i1=0, j1=lnz-s, i=1; i<=lnz; i++) /* push the 1 rows down */
     147      156050 :     if (diagH1[i])
     148       74180 :       gel(p1,++j1) = gel(p2,i);
     149             :     else
     150       81870 :       gel(p2,++i1) = gel(p2,i);
     151       41501 :   for (i=i1+1; i<=lnz; i++) gel(p2,i) = gel(p1,i);
     152             : 
     153             :   /* s = # extra redundant generators taken from H
     154             :    *          zc  col-s  co   zc = col - lnz
     155             :    *       [ 0 |dep |     ]    i = nlze + lnz - s = lig - s
     156             :    *  nlze [--------|  B' ]
     157             :    *       [ 0 | H' |     ]    H' = H minus the s rows with a 1 on diagonal
     158             :    *     i [--------|-----] lig-s           (= "1-rows")
     159             :    *       [   0    | Id  ]
     160             :    *       [        |     ] li */
     161       41501 :   lig -= s; col -= s; lnz -= s;
     162       41501 :   Hnew = cgetg(lnz+1,t_MAT);
     163       41501 :   depnew = cgetg(lnz+1,t_MAT); /* only used if nlze > 0 */
     164       41501 :   Bnew = cgetg(co-col,t_MAT);
     165       41501 :   C = shallowcopy(Cnew);
     166      197551 :   for (j=1,i1=j1=0; j<=lnz+s; j++)
     167             :   {
     168      156050 :     GEN z = gel(H,j);
     169      156050 :     if (diagH1[j])
     170             :     { /* hit exactly s times */
     171       74180 :       i1++; C[i1+col] = Cnew[j+zc];
     172       74180 :       p1 = cgetg(lig+1,t_COL); gel(Bnew,i1) = p1;
     173       74180 :       for (i=1; i<=nlze; i++) gel(p1,i) = gcoeff(dep,i,j);
     174       74180 :       p1 += nlze;
     175             :     }
     176             :     else
     177             :     {
     178       81870 :       j1++; C[j1+zc] = Cnew[j+zc];
     179       81870 :       p1 = cgetg(lnz+1,t_COL); gel(Hnew,j1) = p1;
     180       81870 :       depnew[j1] = dep[j];
     181             :     }
     182      782989 :     for (i=k=1; k<=lnz; i++)
     183      626939 :       if (!diagH1[i]) p1[k++] = z[i];
     184             :   }
     185     1264975 :   for (j=s+1; j<co-col; j++)
     186             :   {
     187     1223474 :     GEN z = gel(B,j-s);
     188     1223474 :     p1 = cgetg(lig+1,t_COL); gel(Bnew,j) = p1;
     189     1223474 :     for (i=1; i<=nlze; i++) gel(p1,i) = gel(z,i);
     190     1223474 :     z += nlze; p1 += nlze;
     191     5720181 :     for (i=k=1; k<=lnz; i++)
     192     4496707 :       if (!diagH1[i]) gel(p1,k++) = gel(z,i);
     193             :   }
     194       41501 :   *ptdep = depnew;
     195       41501 :   *ptC = C;
     196       41501 :   *ptB = Bnew; return Hnew;
     197             : }
     198             : 
     199             : /* for debugging */
     200             : static void
     201           0 : p_mat(GEN mat, GEN perm, long k)
     202             : {
     203           0 :   pari_sp av = avma;
     204           0 :   perm = vecslice(perm, k+1, lg(perm)-1);
     205           0 :   err_printf("Permutation: %Ps\n",perm);
     206           0 :   if (DEBUGLEVEL > 6)
     207           0 :     err_printf("matgen = %Ps\n", zm_to_ZM( rowpermute(mat, perm) ));
     208           0 :   avma = av;
     209           0 : }
     210             : 
     211             : static GEN
     212      636585 : col_dup(long l, GEN col)
     213             : {
     214      636585 :   GEN c = new_chunk(l);
     215      636585 :   memcpy(c,col,l * sizeof(long)); return c;
     216             : }
     217             : 
     218             : /* permutation giving imagecompl(x') | image(x'), x' = transpose of x */
     219             : static GEN
     220       41501 : ZM_rowrankprofile(GEN x, long *nlze)
     221             : {
     222       41501 :   pari_sp av = avma;
     223             :   GEN d, y;
     224             :   long i, j, k, l, r;
     225             : 
     226       41501 :   x = shallowtrans(x); l = lg(x);
     227       41501 :   (void)new_chunk(l); /* HACK */
     228       41501 :   d = ZM_pivots(x,&r); avma = av;
     229       41501 :   *nlze = r;
     230       41501 :   if (!d) return identity_perm(l-1);
     231       38605 :   y = cgetg(l,t_VECSMALL);
     232      251847 :   for (i = j = 1, k = r+1; i<l; i++)
     233      213242 :     if (d[i]) y[k++] = i; else y[j++] = i;
     234       38605 :   return y;
     235             : }
     236             : 
     237             : /* HNF reduce a relation matrix (column operations + row permutation)
     238             : ** Input:
     239             : **   mat = (li-1) x (co-1) matrix of long
     240             : **   C   = r x (co-1) matrix of GEN
     241             : **   perm= permutation vector (length li-1), indexing the rows of mat: easier
     242             : **     to maintain perm than to copy rows. For columns we can do it directly
     243             : **     using e.g. swap(mat[i], mat[j])
     244             : **   k0 = integer. The k0 first lines of mat are dense, the others are sparse.
     245             : ** Output: cf ASCII art in the function body
     246             : **
     247             : ** row permutations applied to perm
     248             : ** column operations applied to C. IN PLACE
     249             : **/
     250             : GEN
     251       24638 : hnfspec_i(GEN mat0, GEN perm, GEN* ptdep, GEN* ptB, GEN* ptC, long k0)
     252             : {
     253             :   pari_sp av;
     254             :   long co, n, s, nlze, lnz, nr, i, j, k, lk0, col, lig, *p;
     255             :   GEN mat;
     256             :   GEN p1, p2, matb, matbnew, vmax, matt, T, extramat, B, C, H, dep, permpro;
     257       24638 :   const long li = lg(perm); /* = lgcols(mat0) */
     258       24638 :   const long CO = lg(mat0);
     259             : 
     260       24638 :   n = 0; /* -Wall */
     261             : 
     262       24638 :   C = *ptC; co = CO;
     263       24638 :   if (co > 300 && co > 1.5 * li)
     264             :   { /* treat the rest at the end */
     265           0 :     co = (long)(1.2 * li);
     266           0 :     setlg(C, co);
     267             :   }
     268             : 
     269       24638 :   if (DEBUGLEVEL>5)
     270             :   {
     271           0 :     err_printf("Entering hnfspec\n");
     272           0 :     p_mat(mat0,perm,0);
     273             :   }
     274       24638 :   matt = cgetg(co, t_MAT); /* dense part of mat (top) */
     275       24638 :   mat = cgetg(co, t_MAT);
     276      661223 :   for (j = 1; j < co; j++)
     277             :   {
     278      636585 :     GEN matj = col_dup(li, gel(mat0,j));
     279      636585 :     p1 = cgetg(k0+1,t_COL); gel(matt,j) = p1; gel(mat,j) = matj;
     280      636585 :     for (i=1; i<=k0; i++) gel(p1,i) = stoi(matj[perm[i]]);
     281             :   }
     282       24638 :   av = avma;
     283             : 
     284       24638 :   i = lig = li-1; col = co-1; lk0 = k0;
     285       24638 :   T = (k0 || (lg(C) > 1 && lgcols(C) > 1))? matid(col): NULL;
     286             :   /* Look for lines with a single non-0 entry, equal to 1 in absolute value */
     287     2423153 :   while (i > lk0 && col)
     288     2373877 :     switch( count(mat,perm[i],col,&n) )
     289             :     {
     290             :       case 0: /* move zero lines between k0+1 and lk0 */
     291       17349 :         lk0++; lswap(perm[i], perm[lk0]);
     292       17349 :         i = lig; continue;
     293             : 
     294             :       case 1: /* move trivial generator between lig+1 and li */
     295      131364 :         lswap(perm[i], perm[lig]);
     296      131364 :         if (T) swap(gel(T,n), gel(T,col));
     297      131364 :         swap(gel(mat,n), gel(mat,col)); p = gel(mat,col);
     298      131364 :         if (p[perm[lig]] < 0) /* = -1 */
     299             :         { /* convert relation -g = 0 to g = 0 */
     300      113202 :           for (i=lk0+1; i<lig; i++) p[perm[i]] = -p[perm[i]];
     301      113202 :           if (T)
     302             :           {
     303      113202 :             p1 = gel(T,col);
     304     5851376 :             for (i=1; ; i++) /* T is a permuted identity: single non-0 entry */
     305     5851376 :               if (signe(gel(p1,i))) { togglesign_safe(&gel(p1,i)); break; }
     306     5738174 :           }
     307             :         }
     308      131364 :         lig--; col--; i = lig; continue;
     309             : 
     310     2225164 :       default: i--;
     311             :     }
     312       24638 :   if (DEBUGLEVEL>5) { err_printf("    after phase1:\n"); p_mat(mat,perm,0); }
     313             : 
     314             : #define absmax(s,z) {long _z; _z = labs(z); if (_z > s) s = _z;}
     315             :   /* Get rid of all lines containing only 0 and +/- 1, keeping track of column
     316             :    * operations in T. Leave the rows 1..lk0 alone [up to k0, coefficient
     317             :    * explosion, between k0+1 and lk0, row is 0] */
     318       24638 :   s = 0;
     319      268098 :   while (lig > lk0 && col && s < (long)(HIGHBIT>>1))
     320             :   {
     321      752023 :     for (i=lig; i>lk0; i--)
     322      730540 :       if (count(mat,perm[i],col,&n) > 0) break;
     323      240305 :     if (i == lk0) break;
     324             : 
     325             :     /* only 0, +/- 1 entries, at least 2 of them non-zero */
     326      218822 :     lswap(perm[i], perm[lig]);
     327      218822 :     swap(gel(mat,n), gel(mat,col)); p = gel(mat,col);
     328      218822 :     if (T) swap(gel(T,n), gel(T,col));
     329      218822 :     if (p[perm[lig]] < 0)
     330             :     {
     331      155889 :       for (i=lk0+1; i<=lig; i++) p[perm[i]] = -p[perm[i]];
     332      155889 :       if (T) ZV_togglesign(gel(T,col));
     333             :     }
     334     9565590 :     for (j=1; j<col; j++)
     335             :     {
     336     9346768 :       GEN matj = gel(mat,j);
     337             :       long t;
     338     9346768 :       if (! (t = matj[perm[lig]]) ) continue;
     339      591714 :       if (t == 1) {
     340      281295 :         for (i=lk0+1; i<=lig; i++) absmax(s, matj[perm[i]] -= p[perm[i]]);
     341             :       }
     342             :       else { /* t = -1 */
     343      310419 :         for (i=lk0+1; i<=lig; i++) absmax(s, matj[perm[i]] += p[perm[i]]);
     344             :       }
     345      591714 :       if (T) ZC_lincomb1_inplace(gel(T,j), gel(T,col), stoi(-t));
     346             :     }
     347      218822 :     lig--; col--;
     348      218822 :     if (gc_needed(av,3))
     349             :     {
     350           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"hnfspec[1]");
     351           0 :       if (T) T = gerepilecopy(av, T); else avma = av;
     352             :     }
     353             :   }
     354             :   /* As above with lines containing a +/- 1 (no other assumption).
     355             :    * Stop when single precision becomes dangerous */
     356       24638 :   vmax = cgetg(co,t_VECSMALL);
     357      311037 :   for (j=1; j<=col; j++)
     358             :   {
     359      286399 :     GEN matj = gel(mat,j);
     360      286399 :     for (s=0, i=lk0+1; i<=lig; i++) absmax(s, matj[i]);
     361      286399 :     vmax[j] = s;
     362             :   }
     363      119123 :   while (lig > lk0 && col)
     364             :   {
     365      109765 :     for (i=lig; i>lk0; i--)
     366      103224 :       if ( (n = count2(mat,perm[i],col)) ) break;
     367       76388 :     if (i == lk0) break;
     368             : 
     369       69847 :     lswap(vmax[n], vmax[col]);
     370       69847 :     lswap(perm[i], perm[lig]);
     371       69847 :     swap(gel(mat,n), gel(mat,col)); p = gel(mat,col);
     372       69847 :     if (T) swap(gel(T,n), gel(T,col));
     373       69847 :     if (p[perm[lig]] < 0)
     374             :     {
     375       31811 :       for (i=lk0+1; i<=lig; i++) p[perm[i]] = -p[perm[i]];
     376       31811 :       if (T) ZV_togglesign(gel(T,col));
     377             :     }
     378     1170892 :     for (j=1; j<col; j++)
     379             :     {
     380     1101045 :       GEN matj = gel(mat,j);
     381             :       long t;
     382     1101045 :       if (! (t = matj[perm[lig]]) ) continue;
     383      591130 :       if (vmax[col] && (ulong)labs(t) >= (HIGHBIT-vmax[j]) / vmax[col])
     384           0 :         goto END2;
     385             : 
     386      591130 :       for (s=0, i=lk0+1; i<=lig; i++) absmax(s, matj[perm[i]] -= t*p[perm[i]]);
     387      591130 :       vmax[j] = s;
     388      591130 :       if (T) ZC_lincomb1_inplace(gel(T,j), gel(T,col), stoi(-t));
     389             :     }
     390       69847 :     lig--; col--;
     391       69847 :     if (gc_needed(av,3))
     392             :     {
     393           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"hnfspec[2]");
     394           0 :       gerepileall(av, T? 2: 1, &vmax, &T);
     395             :     }
     396             :   }
     397             : 
     398             : END2: /* clean up mat: remove everything to the right of the 1s on diagonal */
     399             :   /* go multiprecision first */
     400       24638 :   matb = cgetg(co,t_MAT); /* bottom part (complement of matt) */
     401      661223 :   for (j=1; j<co; j++)
     402             :   {
     403      636585 :     GEN matj = gel(mat,j);
     404      636585 :     p1 = cgetg(li-k0,t_COL); gel(matb,j) = p1;
     405      636585 :     p1 -= k0;
     406      636585 :     for (i=k0+1; i<li; i++) gel(p1,i) = stoi(matj[perm[i]]);
     407             :   }
     408       24638 :   if (DEBUGLEVEL>5)
     409             :   {
     410           0 :     err_printf("    after phase2:\n");
     411           0 :     p_mat(mat,perm,lk0);
     412             :   }
     413      420326 :   for (i=li-2; i>lig; i--)
     414             :   {
     415      395688 :     long h, i0 = i - k0, k = i + co-li;
     416      395688 :     GEN Bk = gel(matb,k);
     417    13321423 :     for (j=k+1; j<co; j++)
     418             :     {
     419    12925735 :       GEN Bj = gel(matb,j), v = gel(Bj,i0);
     420    12925735 :       s = signe(v); if (!s) continue;
     421             : 
     422     2618816 :       gel(Bj,i0) = gen_0;
     423     2618816 :       if (is_pm1(v))
     424             :       {
     425     1492395 :         if (s > 0) /* v = 1 */
     426      740043 :         { for (h=1; h<i0; h++) gel(Bj,h) = subii(gel(Bj,h), gel(Bk,h)); }
     427             :         else /* v = -1 */
     428      752352 :         { for (h=1; h<i0; h++) gel(Bj,h) = addii(gel(Bj,h), gel(Bk,h)); }
     429             :       }
     430             :       else {
     431     1126421 :         for (h=1; h<i0; h++) gel(Bj,h) = subii(gel(Bj,h), mulii(v,gel(Bk,h)));
     432             :       }
     433     2618816 :       if (T) ZC_lincomb1_inplace(gel(T,j), gel(T,k), negi(v));
     434     2618816 :       if (gc_needed(av,3))
     435             :       {
     436          19 :         if(DEBUGMEM>1) pari_warn(warnmem,"hnfspec[3], (i,j) = %ld,%ld", i,j);
     437          19 :         for (h=1; h<co; h++) setlg(matb[h], i0+1); /* bottom can be forgotten */
     438          19 :         gerepileall(av, T? 2: 1, &matb, &T);
     439          19 :         Bk = gel(matb,k);
     440             :       }
     441             :     }
     442             :   }
     443       24638 :   for (j=1; j<co; j++) setlg(matb[j], lig-k0+1); /* bottom can be forgotten */
     444       24638 :   gerepileall(av, T? 2: 1, &matb, &T);
     445       24638 :   if (DEBUGLEVEL>5) err_printf("    matb cleaned up (using Id block)\n");
     446             : 
     447       24638 :   nlze = lk0 - k0;  /* # of 0 rows */
     448       24638 :   lnz = lig-nlze+1; /* 1 + # of non-0 rows (!= 0...0 1 0 ... 0) */
     449       24638 :   if (T) matt = ZM_mul(matt,T); /* update top rows */
     450       24638 :   extramat = cgetg(col+1,t_MAT); /* = new C minus the 0 rows */
     451      241190 :   for (j=1; j<=col; j++)
     452             :   {
     453      216552 :     GEN z = gel(matt,j);
     454      216552 :     GEN t = (gel(matb,j)) + nlze - k0;
     455      216552 :     p2=cgetg(lnz,t_COL); gel(extramat,j) = p2;
     456      216552 :     for (i=1; i<=k0; i++) gel(p2,i) = gel(z,i); /* top k0 rows */
     457      216552 :     for (   ; i<lnz; i++) gel(p2,i) = gel(t,i); /* other non-0 rows */
     458             :   }
     459       24638 :   if (!col) {
     460         357 :     permpro = identity_perm(lnz);
     461         357 :     nr = lnz;
     462             :   }
     463             :   else
     464       24281 :     permpro = ZM_rowrankprofile(extramat, &nr);
     465             :   /* lnz = lg(permpro) */
     466       24638 :   if (nlze)
     467             :   { /* put the nlze 0 rows (trivial generators) at the top */
     468        3881 :     p1 = new_chunk(lk0+1);
     469        3881 :     for (i=1; i<=nlze; i++) p1[i] = perm[i + k0];
     470        3881 :     for (   ; i<=lk0; i++)  p1[i] = perm[i - nlze];
     471        3881 :     for (i=1; i<=lk0; i++)  perm[i] = p1[i];
     472             :   }
     473             :   /* sort other rows according to permpro (nr redundant generators first) */
     474       24638 :   p1 = new_chunk(lnz); p2 = perm + nlze;
     475       24638 :   for (i=1; i<lnz; i++) p1[i] = p2[permpro[i]];
     476       24638 :   for (i=1; i<lnz; i++) p2[i] = p1[i];
     477             :   /* perm indexes the rows of mat
     478             :    *   |_0__|__redund__|__dense__|__too big__|_____done______|
     479             :    *   0  nlze                              lig             li
     480             :    *         \___nr___/ \___k0__/
     481             :    *         \____________lnz ______________/
     482             :    *
     483             :    *               col   co
     484             :    *       [dep     |     ]
     485             :    *    i0 [--------|  B  ] (i0 = nlze + nr)
     486             :    *       [matbnew |     ] matbnew has maximal rank = lnz-1 - nr
     487             :    * mat = [--------|-----] lig
     488             :    *       [   0    | Id  ]
     489             :    *       [        |     ] li */
     490             : 
     491       24638 :   matbnew = cgetg(col+1,t_MAT); /* dense+toobig, maximal rank. For hnffinal */
     492       24638 :   dep    = cgetg(col+1,t_MAT); /* rows dependent from the ones in matbnew */
     493      241190 :   for (j=1; j<=col; j++)
     494             :   {
     495      216552 :     GEN z = gel(extramat,j);
     496      216552 :     p1 = cgetg(nlze+nr+1,t_COL); gel(dep,j) = p1;
     497      216552 :     p2 = cgetg(lnz-nr,t_COL); gel(matbnew,j) = p2;
     498      216552 :     for (i=1; i<=nlze; i++) gel(p1,i) = gen_0;
     499      216552 :     p1 += nlze; for (i=1; i<=nr; i++) p1[i] = z[permpro[i]];
     500      216552 :     p2 -= nr;   for (   ; i<lnz; i++) p2[i] = z[permpro[i]];
     501             :   }
     502             : 
     503             :   /* redundant generators in terms of the genuine generators
     504             :    * (x_i) = - (g_i) B */
     505       24638 :   B = cgetg(co-col,t_MAT);
     506      444671 :   for (j=col+1; j<co; j++)
     507             :   {
     508      420033 :     GEN y = gel(matt,j);
     509      420033 :     GEN z = gel(matb,j);
     510      420033 :     p1=cgetg(lig+1,t_COL); gel(B,j-col) = p1;
     511      420033 :     for (i=1; i<=nlze; i++) gel(p1,i) = gel(z,i);
     512      420033 :     p1 += nlze; z += nlze-k0;
     513     3626573 :     for (k=1; k<lnz; k++)
     514             :     {
     515     3206540 :       i = permpro[k];
     516     3206540 :       gel(p1,k) = (i <= k0)? gel(y,i): gel(z,i);
     517             :     }
     518             :   }
     519       24638 :   if (T) C = typ(C)==t_MAT? RgM_mul(C,T): RgV_RgM_mul(C,T);
     520       24638 :   gerepileall(av, 4, &matbnew, &B, &dep, &C);
     521       24638 :   *ptdep = dep;
     522       24638 :   *ptB = B;
     523       24638 :   H = hnffinal(matbnew, perm, ptdep, ptB, &C);
     524       24638 :   if (CO > co)
     525             :   { /* treat the rest, N cols at a time (hnflll slow otherwise) */
     526           0 :     const long N = 300;
     527           0 :     long a, L = CO - co, l = minss(L, N); /* L columns to add */
     528           0 :     GEN CC = *ptC, m0 = mat0;
     529           0 :     setlg(CC, CO); /* restore */
     530           0 :     CC += co-1;
     531           0 :     m0 += co-1;
     532           0 :     for (a = l;;)
     533             :     {
     534             :       GEN MAT, emb;
     535           0 :       gerepileall(av, 4, &H,&C,ptB,ptdep);
     536           0 :       MAT = cgetg(l + 1, t_MAT);
     537           0 :       emb = cgetg(l + 1, typ(C));
     538           0 :       for (j = 1 ; j <= l; j++)
     539             :       {
     540           0 :         gel(MAT,j) = gel(m0,j);
     541           0 :         emb[j] = CC[j];
     542             :       }
     543           0 :       H = hnfadd_i(H, perm, ptdep, ptB, &C, MAT, emb);
     544           0 :       if (a == L) break;
     545           0 :       CC += l;
     546           0 :       m0 += l;
     547           0 :       a += l; if (a > L) { l = L - (a - l); a = L; }
     548           0 :     }
     549             :   }
     550       24638 :   *ptC = C; return H;
     551             : }
     552             : 
     553             : GEN
     554         322 : hnfspec(GEN mat, GEN perm, GEN* ptdep, GEN* ptB, GEN* ptC, long k0)
     555             : {
     556         322 :   pari_sp av = avma;
     557         322 :   GEN H = hnfspec_i(mat, perm, ptdep, ptB, ptC, k0);
     558         322 :   gerepileall(av, 4, ptC, ptdep, ptB, &H); return H;
     559             : }
     560             : 
     561             : /* HNF reduce x, apply same transforms to C */
     562             : GEN
     563         322 : mathnfspec(GEN x, GEN *ptperm, GEN *ptdep, GEN *ptB, GEN *ptC)
     564             : {
     565         322 :   long i,j,k,ly,lx = lg(x);
     566             :   GEN z, perm;
     567         322 :   if (lx == 1) return cgetg(1, t_MAT);
     568         322 :   ly = lgcols(x);
     569         322 :   *ptperm = perm = identity_perm(ly-1);
     570         322 :   z = cgetg(lx,t_MAT);
     571        2604 :   for (i=1; i<lx; i++)
     572             :   {
     573        2282 :     GEN C = cgetg(ly,t_COL), D = gel(x,i);
     574        2282 :     gel(z,i) = C;
     575       88340 :     for (j=1; j<ly; j++)
     576             :     {
     577       86058 :       GEN d = gel(D,j);
     578       86058 :       if (is_bigint(d)) goto TOOLARGE;
     579       86058 :       C[j] = itos(d);
     580             :     }
     581             :   }
     582             :   /*  [ dep |     ]
     583             :    *  [-----|  B  ]
     584             :    *  [  H  |     ]
     585             :    *  [-----|-----]
     586             :    *  [  0  | Id  ] */
     587         322 :   return hnfspec(z,perm, ptdep, ptB, ptC, 0);
     588             : 
     589             : TOOLARGE:
     590           0 :   if (lg(*ptC) > 1 && lgcols(*ptC) > 1)
     591           0 :     pari_err_IMPL("mathnfspec with large entries");
     592           0 :   x = ZM_hnf(x); lx = lg(x); j = ly; k = 0;
     593           0 :   for (i=1; i<ly; i++)
     594             :   {
     595           0 :     if (equali1(gcoeff(x,i,i + lx-ly)))
     596           0 :       perm[--j] = i;
     597             :     else
     598           0 :       perm[++k] = i;
     599             :   }
     600           0 :   setlg(perm,k+1);
     601           0 :   x = rowpermute(x, perm); /* upper part */
     602           0 :   setlg(perm,ly);
     603           0 :   *ptB = vecslice(x, j+lx-ly, lx-1);
     604           0 :   setlg(x, j);
     605           0 :   *ptdep = rowslice(x, 1, lx-ly);
     606           0 :   return rowslice(x, lx-ly+1, k); /* H */
     607             : }
     608             : 
     609             : /* add new relations to a matrix treated by hnfspec (extramat / extraC) */
     610             : GEN
     611       17220 : hnfadd_i(GEN H, GEN perm, GEN* ptdep, GEN* ptB, GEN* ptC, /* cf hnfspec */
     612             :        GEN extramat,GEN extraC)
     613             : {
     614       17220 :   GEN matb, extratop, Cnew, permpro, B = *ptB, C = *ptC, dep = *ptdep;
     615             :   long i, lH, lB, li, lig, co, col, nlze;
     616             : 
     617       17220 :   if (lg(extramat) == 1) return H;
     618       17220 :   co = lg(C)-1;
     619       17220 :   lH = lg(H)-1;
     620       17220 :   lB = lg(B)-1;
     621       17220 :   li = lg(perm)-1;
     622       17220 :   lig = li - lB;
     623       17220 :   col = co - lB;
     624       17220 :   nlze = lig - lH;
     625             : 
     626             :  /*               col    co
     627             :   *       [ 0 |dep |     ]
     628             :   *  nlze [--------|  B  ]
     629             :   *       [ 0 | H  |     ]
     630             :   *       [--------|-----] lig
     631             :   *       [   0    | Id  ]
     632             :   *       [        |     ] li */
     633       17220 :   extratop = zm_to_ZM( rowslicepermute(extramat, perm, 1, lig) );
     634       17220 :   if (li != lig)
     635             :   { /* zero out bottom part, using the Id block */
     636       17220 :     GEN A = vecslice(C, col+1, co);
     637       17220 :     GEN c = rowslicepermute(extramat, perm, lig+1, li);
     638       17220 :     extraC   = gsub(extraC, typ(A)==t_MAT? RgM_zm_mul(A, c): RgV_zm_mul(A,c));
     639       17220 :     extratop = ZM_sub(extratop, ZM_zm_mul(B, c));
     640             :   }
     641             : 
     642       17220 :   extramat = shallowconcat(extratop, vconcat(dep, H));
     643       17220 :   Cnew     = shallowconcat(extraC, vecslice(C, col-lH+1, co));
     644       17220 :   if (DEBUGLEVEL>5) err_printf("    1st phase done\n");
     645       17220 :   permpro = ZM_rowrankprofile(extramat, &nlze);
     646       17220 :   extramat = rowpermute(extramat, permpro);
     647       17220 :   *ptB     = rowpermute(B,        permpro);
     648       17220 :   permpro = vecsmallpermute(perm, permpro);
     649       17220 :   for (i=1; i<=lig; i++) perm[i] = permpro[i]; /* perm o= permpro */
     650             : 
     651       17220 :   *ptdep  = rowslice(extramat, 1, nlze);
     652       17220 :   matb    = rowslice(extramat, nlze+1, lig);
     653       17220 :   if (DEBUGLEVEL>5) err_printf("    2nd phase done\n");
     654       17220 :   H = hnffinal(matb,perm,ptdep,ptB,&Cnew);
     655       17220 :   *ptC = shallowconcat(vecslice(C, 1, col-lH), Cnew);
     656       17220 :   return H;
     657             : }
     658             : 
     659             : GEN
     660           0 : hnfadd(GEN H, GEN perm, GEN* ptdep, GEN* ptB, GEN* ptC, /* cf hnfspec */
     661             :        GEN extramat,GEN extraC)
     662             : {
     663           0 :   pari_sp av = avma;
     664           0 :   H = hnfadd_i(H, perm, ptdep, ptB, ptC, ZM_to_zm(extramat), extraC);
     665           0 :   gerepileall(av, 4, ptC, ptdep, ptB, &H); return H;
     666             : }
     667             : 
     668             : /* zero aj = Aij (!= 0)  using  ak = Aik (maybe 0), via linear combination of
     669             :  * A[j] and A[k] of determinant 1. If U != NULL, likewise update its columns */
     670             : static void
     671    19792543 : ZC_elem(GEN aj, GEN ak, GEN A, GEN U, long j, long k)
     672             : {
     673             :   GEN p1,u,v,d;
     674             : 
     675    19792543 :   if (!signe(ak)) {
     676     3628238 :     swap(gel(A,j), gel(A,k));
     677     3628238 :     if (U) swap(gel(U,j), gel(U,k));
     678    20451733 :     return;
     679             :   }
     680    16164305 :   d = bezout(aj,ak,&u,&v);
     681             :   /* frequent special case (u,v) = (1,0) or (0,1) */
     682    16164305 :   if (!signe(u))
     683             :   { /* ak | aj */
     684     6196980 :     p1 = diviiexact(aj,ak); togglesign(p1);
     685     6196980 :     ZC_lincomb1_inplace(gel(A,j), gel(A,k), p1);
     686     6196980 :     if (U)
     687      583977 :       ZC_lincomb1_inplace(gel(U,j), gel(U,k), p1);
     688     6196980 :     return;
     689             :   }
     690     9967325 :   if (!signe(v))
     691             :   { /* aj | ak */
     692     6998277 :     p1 = diviiexact(ak,aj); togglesign(p1);
     693     6998277 :     ZC_lincomb1_inplace(gel(A,k), gel(A,j), p1);
     694     6998277 :     swap(gel(A,j), gel(A,k));
     695     6998277 :     if (U) {
     696      105899 :       ZC_lincomb1_inplace(gel(U,k), gel(U,j), p1);
     697      105899 :       swap(gel(U,j), gel(U,k));
     698             :     }
     699     6998277 :     return;
     700             :   }
     701             : 
     702     2969048 :   if (!is_pm1(d)) { aj = diviiexact(aj, d); ak = diviiexact(ak, d); }
     703     2969048 :   p1 = gel(A,k); aj = negi(aj); /* NOT togglesign */
     704     2969048 :   gel(A,k) = ZC_lincomb(u,v, gel(A,j),p1);
     705     2969048 :   gel(A,j) = ZC_lincomb(aj,ak, p1,gel(A,j));
     706     2969048 :   if (U)
     707             :   {
     708      346894 :     p1 = gel(U,k);
     709      346894 :     gel(U,k) = ZC_lincomb(u,v, gel(U,j),p1);
     710      346894 :     gel(U,j) = ZC_lincomb(aj,ak, p1,gel(U,j));
     711             :   }
     712             : }
     713             : 
     714             : INLINE int
     715        1001 : is_RgX(GEN a, long v) { return typ(a) == t_POL && varn(a)==v; }
     716             : /* set u,v such that au + bv = gcd(a,b), divide a,b by the gcd */
     717             : static GEN
     718         266 : gbezout_step(GEN *pa, GEN *pb, GEN *pu, GEN *pv, long vx)
     719             : {
     720         266 :   GEN a = *pa, b = *pb, d;
     721         266 :   if (gequal0(a))
     722             :   {
     723           0 :     *pa = gen_0; *pu = gen_0;
     724           0 :     *pb = gen_1; *pv = gen_1; return b;
     725             :   }
     726         266 :   a = is_RgX(a,vx)? RgX_renormalize(a): scalarpol(a, vx);
     727         266 :   b = is_RgX(b,vx)? RgX_renormalize(b): scalarpol(b, vx);
     728         266 :   d = RgX_extgcd(a,b, pu,pv);
     729         266 :   if (degpol(d)) { a = RgX_div(a, d); b = RgX_div(b, d); }
     730         154 :   else if (typ(gel(d,2)) == t_REAL && lg(gel(d,2)) <= 3)
     731             : #if 1
     732             :   { /* possible accuracy problem */
     733           0 :     GEN D = RgX_gcd_simple(a,b);
     734           0 :     if (degpol(D)) {
     735           0 :       D = RgX_Rg_div(D, leading_coeff(D));
     736           0 :       a = RgX_div(a, D);
     737           0 :       b = RgX_div(b, D);
     738           0 :       d = RgX_extgcd(a,b, pu,pv); /* retry now */
     739           0 :       d = RgX_mul(d, D);
     740             :     }
     741             :   }
     742             : #else
     743             :   { /* less stable */
     744             :     d = RgX_extgcd_simple(a,b, pu,pv);
     745             :     if (degpol(d)) { a = RgX_div(a, d); b = RgX_div(b, d); }
     746             :   }
     747             : #endif
     748         266 :   *pa = a;
     749         266 :   *pb = b; return d;
     750             : }
     751             : static GEN
     752      236720 : col_mul(GEN x, GEN c)
     753             : {
     754      236720 :   if (typ(x) == t_INT)
     755             :   {
     756      236720 :     long s = signe(x);
     757      236720 :     if (!s) return NULL;
     758      177848 :     if (is_pm1(x)) return (s > 0)? c: RgC_neg(c);
     759             :   }
     760       32151 :   return RgC_Rg_mul(c, x);
     761             : }
     762             : static void
     763           0 : do_zero(GEN x)
     764             : {
     765           0 :   long i, lx = lg(x);
     766           0 :   for (i=1; i<lx; i++) gel(x,i) = gen_0;
     767           0 : }
     768             : 
     769             : /* (c1, c2) *= [u,-b; v,a] */
     770             : static void
     771       59180 : update(GEN u, GEN v, GEN a, GEN b, GEN *c1, GEN *c2)
     772             : {
     773             :   GEN p1,p2;
     774             : 
     775       59180 :   u = col_mul(u,*c1);
     776       59180 :   v = col_mul(v,*c2);
     777       59180 :   if (u) p1 = v? gadd(u,v): u;
     778         294 :   else   p1 = v? v: NULL;
     779             : 
     780       59180 :   a = col_mul(a,*c2);
     781       59180 :   b = col_mul(gneg_i(b),*c1);
     782       59180 :   if (a) p2 = b? RgC_add(a,b): a;
     783           0 :   else   p2 = b? b: NULL;
     784             : 
     785       59180 :   if (!p1) do_zero(*c1); else *c1 = p1;
     786       59180 :   if (!p2) do_zero(*c2); else *c2 = p2;
     787       59180 : }
     788             : 
     789             : /* zero aj = Aij (!= 0)  using  ak = Aik (maybe 0), via linear combination of
     790             :  * A[j] and A[k] of determinant 1. If U != NULL, likewise update its columns */
     791             : static void
     792           7 : RgC_elem(GEN aj, GEN ak, GEN A, GEN V, long j, long k, long li, long vx)
     793             : {
     794           7 :   GEN u,v, d = gbezout_step(&aj, &ak, &u, &v, vx);
     795             :   long l;
     796             :   /* (A[,k], A[,j]) *= [v, -aj; u, ak ] */
     797          14 :   for (l = 1; l < li; l++)
     798             :   {
     799           7 :     GEN t = gadd(gmul(u,gcoeff(A,l,j)), gmul(v,gcoeff(A,l,k)));
     800           7 :     gcoeff(A,l,j) = gsub(gmul(ak,gcoeff(A,l,j)), gmul(aj,gcoeff(A,l,k)));
     801           7 :     gcoeff(A,l,k) = t;
     802             :   }
     803           7 :   gcoeff(A,li,j) = gen_0;
     804           7 :   gcoeff(A,li,k) = d;
     805           7 :   if (V) update(v,u,ak,aj,(GEN*)(V+k),(GEN*)(V+j));
     806           7 : }
     807             : 
     808             : /* reduce A[i,j] mod A[i,j0] for j=j0+1... via column operations */
     809             : static void
     810     2000785 : ZM_reduce(GEN A, GEN U, long i, long j0)
     811             : {
     812     2000785 :   long j, lA = lg(A);
     813     2000785 :   GEN d = gcoeff(A,i,j0);
     814     2000785 :   if (signe(d) < 0)
     815             :   {
     816        4041 :     ZV_neg_inplace(gel(A,j0));
     817        4041 :     if (U) ZV_togglesign(gel(U,j0));
     818        4041 :     d = gcoeff(A,i,j0);
     819             :   }
     820     5172587 :   for (j=j0+1; j<lA; j++)
     821             :   {
     822     3171802 :     GEN q = truedivii(gcoeff(A,i,j), d);
     823     3171802 :     if (!signe(q)) continue;
     824             : 
     825      553021 :     togglesign(q);
     826      553021 :     ZC_lincomb1_inplace(gel(A,j), gel(A,j0), q);
     827      553021 :     if (U) ZC_lincomb1_inplace(gel(U,j), gel(U,j0), q);
     828             :   }
     829     2000785 : }
     830             : 
     831             : /* normalize T as if it were a t_POL in variable v */
     832             : static GEN
     833         210 : normalize_as_RgX(GEN T, long v, GEN *pd)
     834             : {
     835             :   GEN d;
     836         210 :   if (!is_RgX(T,v)) { *pd = T; return gen_1; }
     837         175 :   d = leading_coeff(T);
     838         350 :   while (gequal0(d) || (typ(d) == t_REAL && lg(d) == 3
     839           0 :                         && gexpo(T) - expo(d) > (long)BITS_IN_LONG)) {
     840          14 :     T = normalizepol_lg(T, lg(T)-1);
     841          14 :     if (!signe(T)) { *pd = gen_1; return T; }
     842           0 :     d = leading_coeff(T);
     843             :   }
     844         161 :   if (degpol(T)) T = RgX_Rg_div(T,d); else { d = gel(T,2); T = gen_1; }
     845         161 :   *pd = d; return T;
     846             : }
     847             : /* reduce A[i,j] mod A[i,j0] for j=j0+1... via column operations */
     848             : static void
     849          42 : RgM_reduce(GEN A, GEN U, long i, long j0, long vx)
     850             : {
     851          42 :   long j, lA = lg(A);
     852          42 :   GEN d, T = normalize_as_RgX(gcoeff(A,i,j0), vx, &d);
     853          42 :   if (U && !gequal1(d)) gel(U,j0) = RgC_Rg_div(gel(U,j0), d);
     854          42 :   gcoeff(A,i,j0) = T;
     855             : 
     856          49 :   for (j=j0+1; j<lA; j++)
     857             :   {
     858           7 :     GEN t = gcoeff(A,i,j), q;
     859           7 :     if (gequal0(t)) continue;
     860           7 :     if (T == gen_1)
     861           0 :       q = t;
     862           7 :     else if (is_RgX(t,vx))
     863           7 :       q = RgX_div(t, T);
     864           0 :     else continue;
     865             : 
     866           7 :     if (gequal0(q)) continue;
     867           0 :     gel(A,j) = RgC_sub(gel(A,j), RgC_Rg_mul(gel(A,j0), q));
     868           0 :     if (U) gel(U,j) = RgC_sub(gel(U,j), RgC_Rg_mul(gel(U,j0), q));
     869             :   }
     870          42 : }
     871             : 
     872             : /* A,B square integral in upper HNF, of the same dimension > 0. Return Au
     873             :  * in Z^n (v in Z^n not computed), such that Au + Bv = [1, 0, ..., 0] */
     874             : GEN
     875      123284 : hnfmerge_get_1(GEN A, GEN B)
     876             : {
     877      123284 :   pari_sp av = avma;
     878      123284 :   long j, k, l = lg(A), lb;
     879      123284 :   GEN b, U = cgetg(l + 1, t_MAT), C = cgetg(l + 1, t_VEC);
     880             : 
     881      123284 :   b = gcoeff(B,1,1); lb = lgefint(b);
     882      253083 :   for (j = 1; j < l; j++)
     883             :   {
     884             :     GEN t;
     885      253083 :     long c = j+1;
     886      253083 :     gel(U,j) = col_ei(l-1, j);
     887      253083 :     gel(U,c) = zerocol(l-1); /* dummy */
     888      253083 :     gel(C,j) = vecslice(gel(A,j), 1,j);
     889      253083 :     gel(C,c) = vecslice(gel(B,j), 1,j);
     890      731004 :     for (k = j; k > 0; k--)
     891             :     {
     892      477921 :       t = gcoeff(C,k,c);
     893      477921 :       if (gequal0(t)) continue;
     894      418036 :       setlg(C[c], k+1);
     895      418036 :       ZC_elem(t, gcoeff(C,k,k), C, U, c, k);
     896      418036 :       if (lgefint(gcoeff(C,k,k)) > lb) gel(C,k) = FpC_red(gel(C,k), b);
     897      418036 :       if (j > 4)
     898             :       {
     899       43820 :         GEN u = gel(U,k);
     900             :         long h;
     901      614173 :         for (h=1; h<l; h++)
     902      570353 :           if (lgefint(gel(u,h)) > lb) gel(u,h) = remii(gel(u,h), b);
     903             :       }
     904             :     }
     905      253083 :     if (j == 1)
     906      123284 :       t = gcoeff(C,1,1);
     907             :     else
     908             :     {
     909             :       GEN u;
     910      129799 :       t = bezout(gcoeff(C,1,1), b, &u, NULL); /* >= 0 */
     911      129799 :       if (signe(u) && !equali1(u)) gel(U,1) = ZC_Z_mul(gel(U,1), u);
     912      129799 :       gcoeff(C,1,1) = t;
     913             :     }
     914      253083 :     if (equali1(t)) break;
     915             :   }
     916      123284 :   if (j >= l) return NULL;
     917      123284 :   b = lcmii(gcoeff(A,1,1),b);
     918      123284 :   A = FpC_red(ZM_ZC_mul(A,gel(U,1)), b);
     919      123284 :   return gerepileupto(av, FpC_center(A, b, shifti(b,-1)));
     920             : }
     921             : 
     922             : /* remove the first r columns */
     923             : static void
     924       96236 : remove_0cols(long r, GEN *pA, GEN *pB, long remove)
     925             : {
     926       96236 :   GEN A = *pA, B = *pB;
     927       96236 :   long l = lg(A);
     928       96236 :   A += r; A[0] = evaltyp(t_MAT) | evallg(l-r);
     929       96236 :   if (B && remove == 2) { B += r; B[0] = A[0]; }
     930       96236 :   *pA = A; *pB = B;
     931       96236 : }
     932             : 
     933             : /* Inefficient compared to hnfall. 'remove' = throw away lin.dep columns */
     934             : static GEN
     935       17413 : hnf_i(GEN A, int remove)
     936             : {
     937       17413 :   pari_sp av0 = avma, av;
     938             :   long s, n, m, j, k, li, def, ldef;
     939             : 
     940       17413 :   RgM_dimensions(A, &m, &n);
     941       17413 :   if (!n) return cgetg(1,t_MAT);
     942       17308 :   av = avma;
     943       17308 :   A = RgM_shallowcopy(A);
     944       17308 :   def = n; ldef = (m>n)? m-n: 0;
     945       47370 :   for (li=m; li>ldef; li--)
     946             :   {
     947       74757 :     for (j=def-1; j; j--)
     948             :     {
     949       44695 :       GEN a = gcoeff(A,li,j);
     950       44695 :       if (!signe(a)) continue;
     951             : 
     952             :       /* zero a = Aij  using  b = Aik */
     953       24511 :       k = (j==1)? def: j-1;
     954       24511 :       ZC_elem(a,gcoeff(A,li,k), A,NULL, j,k);
     955       24511 :       if (gc_needed(av,1))
     956             :       {
     957           0 :         if (DEBUGMEM>1) pari_warn(warnmem,"ZM_hnf[1]. li=%ld",li);
     958           0 :         A = gerepilecopy(av, A);
     959             :       }
     960             :     }
     961       30062 :     s = signe(gcoeff(A,li,def));
     962       30062 :     if (s)
     963             :     {
     964       30048 :       if (s < 0) ZV_neg_inplace(gel(A,def));
     965       30048 :       ZM_reduce(A, NULL, li,def);
     966       30048 :       def--;
     967             :     }
     968             :     else
     969          14 :       if (ldef) ldef--;
     970       30062 :     if (gc_needed(av,1))
     971             :     {
     972           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"ZM_hnf[2]. li=%ld",li);
     973           0 :       A = gerepilecopy(av, A);
     974             :     }
     975             :   }
     976             :   /* rank A = n - def */
     977       17308 :   if (remove) { GEN B = NULL; remove_0cols(def, &A, &B, remove); }
     978       17308 :   return gerepileupto(av0, ZM_copy(A));
     979             : }
     980             : 
     981             : GEN
     982       19588 : ZM_hnf(GEN x) { return lg(x) > 8? ZM_hnfall(x, NULL, 1): hnf_i(x, 1); }
     983             : 
     984             : /* u*z[1..k] mod p, in place */
     985             : static void
     986     2579137 : FpV_Fp_mul_part_ip(GEN z, GEN u, GEN p, long k)
     987             : {
     988             :   long i;
     989     2579137 :   if (is_pm1(u)) {
     990      205471 :     if (signe(u) > 0) {
     991      864975 :       for (i = 1; i <= k; i++)
     992      676416 :         if (signe(gel(z,i))) gel(z,i) = modii(gel(z,i), p);
     993             :     } else {
     994       43582 :       for (i = 1; i <= k; i++)
     995       26670 :         if (signe(gel(z,i))) gel(z,i) = modii(negi(gel(z,i)), p);
     996             :     }
     997             :   }
     998             :   else {
     999     9635504 :     for (i = 1; i <= k; i++)
    1000     7261838 :       if (signe(gel(z,i))) gel(z,i) = Fp_mul(u,gel(z,i), p);
    1001             :   }
    1002     2579137 : }
    1003             : static void
    1004     9474040 : FpV_red_part_ipvec(GEN z, GEN p, long k)
    1005             : {
    1006             :   long i;
    1007     9474040 :   for (i = 1; i <= k; i++) gel(z,i) = modii(gel(z,i), gel(p,i));
    1008     9474040 : }
    1009             : 
    1010             : /* return x * U, in echelon form (mod p^m), where (det(U),p) = 1.
    1011             :  * If early_abort is set, return NULL as soon as one pivot is 0 (mod p^m) */
    1012             : GEN
    1013       46945 : ZpM_echelon(GEN x, long early_abort, GEN p, GEN pm)
    1014             : {
    1015       46945 :   pari_sp av0 = avma, av;
    1016             :   long m, li, co, i, j, k, def, ldef;
    1017             : 
    1018       46945 :   co = lg(x); if (co == 1) return cgetg(1,t_MAT);
    1019       46945 :   li = lgcols(x);
    1020       46945 :   av = avma;
    1021       46945 :   x = RgM_shallowcopy(x);
    1022       46945 :   m = Z_pval(pm, p);
    1023             : 
    1024       46945 :   ldef = (li > co)? li - co: 0;
    1025      423594 :   for (def = co-1,i = li-1; i > ldef; i--)
    1026             :   {
    1027      377073 :     long vmin = LONG_MAX, kmin = 0;
    1028      377073 :     GEN umin = gen_0, pvmin, q;
    1029     3448364 :     for (k = 1; k <= def; k++)
    1030             :     {
    1031     3144436 :       GEN u = gcoeff(x,i,k);
    1032             :       long v;
    1033     3144436 :       if (!signe(u)) continue;
    1034     1517542 :       v = Z_pvalrem(u, p, &u);
    1035     1517542 :       if (v >= m) gcoeff(x,i,k) = gen_0;
    1036      981757 :       else if (v < vmin) {
    1037      365189 :         vmin = v; kmin = k; umin = u;
    1038      365189 :         if (!vmin) break;
    1039             :       }
    1040             :     }
    1041      377073 :     if (!kmin)
    1042             :     {
    1043      154327 :       if (early_abort) return NULL;
    1044      153903 :       gcoeff(x,i,def) = gen_0;
    1045      153903 :       ldef--;
    1046      153903 :       if (ldef < 0) ldef = 0;
    1047      153903 :       continue;
    1048             :     }
    1049      222746 :     if (kmin != def) swap(gel(x,def), gel(x,kmin));
    1050      222746 :     q = vmin? powiu(p, m-vmin): pm;
    1051             :     /* pivot has valuation vmin */
    1052      222746 :     umin = modii(umin, q);
    1053      222746 :     if (!equali1(umin))
    1054      174652 :       FpV_Fp_mul_part_ip(gel(x,def), Fp_inv(umin,q), pm, i-1);
    1055      222746 :     gcoeff(x, i, def) = pvmin = powiu(p, vmin);
    1056     2587401 :     for (j = def-1; j; j--)
    1057             :     { /* zero x[i, 1..def-1] using x[i,def] = pvmin */
    1058     2364655 :       GEN t, a = gcoeff(x,i,j) = modii(gcoeff(x,i,j), pm);
    1059     2364655 :       if (!signe(a)) continue;
    1060             : 
    1061     1165287 :       t = diviiexact(a, pvmin); togglesign(t);
    1062     1165287 :       ZC_lincomb1_inplace(gel(x,j), gel(x,def), t);
    1063     1165287 :       if (gc_needed(av,1))
    1064             :       {
    1065          14 :         if (DEBUGMEM>1) pari_warn(warnmem,"ZpM_echelon. i=%ld",i);
    1066          14 :         x = gerepilecopy(av, x); pvmin = gcoeff(x,i,def);
    1067             :       }
    1068             :     }
    1069      222746 :     def--;
    1070             :   }
    1071       46521 :   if (co > li)
    1072             :   {
    1073           0 :     x += co - li;
    1074           0 :     x[0] = evaltyp(t_MAT) | evallg(li);
    1075             :   }
    1076       46521 :   return gerepilecopy(av0, x);
    1077             : }
    1078             : GEN
    1079      766404 : zlm_echelon(GEN x, long early_abort, ulong p, ulong pm)
    1080             : {
    1081      766404 :   pari_sp av0 = avma;
    1082             :   long li, co, i, j, k, def, ldef;
    1083             :   ulong m;
    1084             : 
    1085      766404 :   co = lg(x); if (co == 1) return cgetg(1,t_MAT);
    1086      766404 :   li = lgcols(x);
    1087      766404 :   x = Flm_copy(x);
    1088      766404 :   m = u_lval(pm, p);
    1089             : 
    1090      766404 :   ldef = (li > co)? li - co: 0;
    1091     3923083 :   for (def = co-1,i = li-1; i > ldef; i--)
    1092             :   {
    1093     3212736 :     long vmin = LONG_MAX, kmin = 0;
    1094     3212736 :     ulong umin = 0, pvmin, q;
    1095    12432413 :     for (k = 1; k <= def; k++)
    1096             :     {
    1097    11062030 :       ulong u = ucoeff(x,i,k);
    1098             :       long v;
    1099    11062030 :       if (!u) continue;
    1100     5136121 :       v = u_lvalrem(u, p, &u);
    1101     5136121 :       if (v >= (long) m) ucoeff(x,i,k) = 0;
    1102     5136121 :       else if (v < vmin) {
    1103     3217828 :         vmin = v; kmin = k; umin = u;
    1104     3217828 :         if (!vmin) break;
    1105             :       }
    1106             :     }
    1107     3212736 :     if (!kmin)
    1108             :     {
    1109      234493 :       if (early_abort) return NULL;
    1110      178436 :       ucoeff(x,i,def) = 0;
    1111      178436 :       ldef--;
    1112      178436 :       if (ldef < 0) ldef = 0;
    1113      178436 :       continue;
    1114             :     }
    1115     2978243 :     if (kmin != def) swap(gel(x,def), gel(x,kmin));
    1116     2978243 :     q = vmin? upowuu(p, m-vmin): pm;
    1117             :     /* pivot has valuation vmin */
    1118     2978243 :     umin %= q;
    1119     2978243 :     if (umin != 1)
    1120     1539715 :       Flv_Fl_mul_part_inplace(gel(x,def), Fl_inv(umin,q), pm, i-1);
    1121     2978243 :     ucoeff(x, i, def) = pvmin = upowuu(p, vmin);
    1122    18335736 :     for (j = def-1; j; j--)
    1123             :     { /* zero x[i, 1..def-1] using x[i,def] = pvmin */
    1124    15357493 :       ulong t, a = ucoeff(x,i,j);
    1125    15357493 :       if (!a) continue;
    1126             : 
    1127     8549345 :       t = Fl_neg(a / pvmin, q);
    1128     8549345 :       Flc_lincomb1_inplace(gel(x,j), gel(x,def), t, pm);
    1129             :     }
    1130     2978243 :     def--;
    1131             :   }
    1132      710347 :   if (co > li)
    1133             :   {
    1134           0 :     x += co - li;
    1135           0 :     x[0] = evaltyp(t_MAT) | evallg(li);
    1136             :   }
    1137      710347 :   return gerepilecopy(av0, x);
    1138             : }
    1139             : 
    1140             : /* dm = multiple of diag element (usually detint(x))
    1141             :  * flag & hnf_MODID: reduce mod dm * matid [ otherwise as above ].
    1142             :  * flag & hnf_PART: don't reduce once diagonal is known; */
    1143             : 
    1144             : /* x a ZM, dm a t_INT */
    1145             : GEN
    1146     1978359 : ZM_hnfmodall_i(GEN x, GEN dm, long flag)
    1147             : {
    1148             :   pari_sp av;
    1149     1978359 :   const long center = (flag & hnf_CENTER);
    1150     1978359 :   long moddiag = (flag & hnf_MODID);
    1151             :   long li, co, i, j, k, def, ldef;
    1152             :   GEN a, b, p1, p2, u, dm2, LDM;
    1153             : 
    1154     1978359 :   co = lg(x); if (co == 1) return cgetg(1,t_MAT);
    1155     1978303 :   li = lgcols(x); if (li == 1) return cgetg(1,t_MAT);
    1156     1978296 :   if (typ(dm) == t_INT)
    1157             :   {
    1158     1963680 :     if (flag == hnf_MODID)
    1159             :     { /* if log p >> n^3, ispsp dominates hnf ! */
    1160     1847186 :       ulong C = itou_or_0(muliu(sqru(co-1), 2*(li-1)));
    1161     1847186 :       if ((lgefint(dm) == 3 || !C || (ulong)expi(dm) < C) && BPSW_psp(dm))
    1162      482011 :         return ZM_hnfmodprime(x, dm);
    1163             :     }
    1164     1481669 :     LDM= const_vecsmall(li-1, lgefint(dm));
    1165     1481669 :     dm2 = shifti(dm, -1);
    1166     1481669 :     dm = const_vec(li-1,dm);
    1167     1481669 :     dm2= const_vec(li-1,dm2);
    1168             :   }
    1169             :   else
    1170             :   {
    1171       14616 :     if (lg(dm) != li) pari_err_DIM("ZM_hnfmod");
    1172       14616 :     moddiag = 1;
    1173       14616 :     dm2 = cgetg(li, t_VEC);
    1174       14616 :     LDM = cgetg(li, t_VECSMALL);
    1175       46319 :     for (i=1; i<li; i++)
    1176             :     {
    1177       31703 :       gel(dm2,i) = shifti(gel(dm,i),-1);
    1178       31703 :       LDM[i] = lgefint(gel(dm,i));
    1179             :     }
    1180             :   }
    1181     1496285 :   av = avma;
    1182     1496285 :   x = RgM_shallowcopy(x);
    1183             : 
    1184     1496285 :   ldef = 0;
    1185     1496285 :   if (li > co)
    1186             :   {
    1187       25074 :     ldef = li - co;
    1188       25074 :     if (!moddiag)
    1189           7 :       pari_err_DOMAIN("ZM_hnfmod","nb lines",">", strtoGENstr("nb columns"), x);
    1190             :   }
    1191     6306187 :   for (def = co-1,i = li-1; i > ldef; i--,def--)
    1192             :   {
    1193     4809909 :     GEN d = gel(dm,i), d2 = gel(dm2,i);
    1194     4809909 :     gcoeff(x,i,def) = centermodii(gcoeff(x,i,def), d,d2);
    1195    27497063 :     for (j = def-1; j; j--)
    1196             :     {
    1197    22687154 :       gcoeff(x,i,j) = centermodii(gcoeff(x,i,j), d,d2);
    1198    22687154 :       a = gcoeff(x,i,j);
    1199    22687154 :       if (!signe(a)) continue;
    1200             : 
    1201    13854282 :       k = (j==1)? def: j-1;
    1202    13854282 :       gcoeff(x,i,k) = centermodii(gcoeff(x,i,k), d,d2);
    1203    13854282 :       ZC_elem(a,gcoeff(x,i,k), x,NULL, j,k);
    1204    13854282 :       p1 = gel(x,j);
    1205    13854282 :       p2 = gel(x,k);
    1206             :       /* prevent coeffs explosion: reduce mod dm when lg() > ldm */
    1207    83610387 :       for (k = 1; k < i; k++)
    1208             :       {
    1209    69756105 :         if (lgefint(gel(p1,k)) > LDM[k])
    1210    12754678 :           gel(p1,k) = centermodii(gel(p1,k), gel(dm,k),gel(dm2,k));
    1211    69756105 :         if (lgefint(gel(p2,k)) > LDM[k])
    1212     2036606 :           gel(p2,k) = centermodii(gel(p2,k), gel(dm,k),gel(dm2,k));
    1213             :       }
    1214             :     }
    1215     4809909 :     if (gc_needed(av,1))
    1216             :     {
    1217          12 :       if (DEBUGMEM>1) pari_warn(warnmem,"ZM_hnfmod[1]. i=%ld",i);
    1218          12 :       x = gerepilecopy(av, x);
    1219             :     }
    1220     4809909 :     if (moddiag && !signe(gcoeff(x,i,def)))
    1221             :     { /* missing pivot on line i, insert column */
    1222     1170840 :       GEN a = cgetg(co + 1, t_MAT);
    1223     1170840 :       for (k = 1; k <= def; k++) gel(a,k) = gel(x,k);
    1224     1170840 :       gel(a,k++) = Rg_col_ei(gel(dm,i), li-1, i);
    1225     1170840 :       for (     ; k <= co;  k++) gel(a,k) = gel(x,k-1);
    1226     1170840 :       ldef--; if (ldef < 0) ldef = 0;
    1227     1170840 :       co++; def++; x = a;
    1228             :     }
    1229             :   }
    1230     1496278 :   if (co < li)
    1231             :   { /* implies moddiag, add missing diag(dm) components */
    1232       24101 :     GEN a = cgetg(li+1, t_MAT);
    1233       24101 :     for (k = 1; k <= li-co; k++) gel(a,k) = Rg_col_ei(gel(dm,k), li-1, k);
    1234       24101 :     for (i = 1; i < co; i++) gel(a,k-1+i) = gel(x,i);
    1235       24101 :     gel(a,li) = zerocol(li-1); x = a;
    1236             :   }
    1237             :   else
    1238             :   {
    1239     1472177 :     x += co - li;
    1240     1472177 :     x[0] = evaltyp(t_MAT) | evallg(li); /* kill 0 columns */
    1241             :   }
    1242     1496278 :   if (moddiag)
    1243             :   { /* one column extra: an accumulator, discarded at the end */
    1244     1384810 :     if (lg(x) == li) x = shallowconcat(x, zerocol(li-1));
    1245             :     /* add up missing diag(dm) components */
    1246     5887090 :     for (i = li-1; i > 0; i--)
    1247             :     {
    1248     4502280 :       gcoeff(x, i, li) = gel(dm,i);
    1249    17806662 :       for (j = i; j > 0; j--)
    1250             :       {
    1251    13304382 :         GEN a = gcoeff(x, j, li);
    1252    13304382 :         if (!signe(a)) continue;
    1253     4737020 :         ZC_elem(a, gcoeff(x,j,j), x, NULL, li,j);
    1254     4737020 :         FpV_red_part_ipvec(gel(x,li), dm, j-1);
    1255     4737020 :         FpV_red_part_ipvec(gel(x,j),  dm, j-1);
    1256             :       }
    1257     4502280 :       if (gc_needed(av,1))
    1258             :       {
    1259           1 :         if (DEBUGMEM>1) pari_warn(warnmem,"ZM_hnfmod[2]. i=%ld", i);
    1260           1 :         x = gerepilecopy(av, x);
    1261             :       }
    1262             :     }
    1263             :   }
    1264             :   else
    1265             :   {
    1266      111468 :     b = gel(dm,1);
    1267      458822 :     for (i = li-1; i > 0; i--)
    1268             :     {
    1269      347354 :       GEN d = bezout(gcoeff(x,i,i),b, &u,NULL);
    1270      347354 :       gcoeff(x,i,i) = d;
    1271      347354 :       FpV_Fp_mul_part_ip(gel(x,i), u, b, i-1);
    1272      347354 :       if (i > 1) b = diviiexact(b,d);
    1273             :     }
    1274             :   }
    1275     1496278 :   x[0] = evaltyp(t_MAT) | evallg(li); /* kill 0 columns / discard accumulator */
    1276     1496278 :   if (flag & hnf_PART) return x;
    1277             : 
    1278     1496278 :   if (!moddiag)
    1279             :   { /* compute optimal value for dm */
    1280      111468 :     b = cgetg(li, t_VEC); gel(b,1) = gcoeff(x,1,1);
    1281      111468 :     for (i = 2; i < li; i++) gel(b,i) = mulii(gel(b,i-1), gcoeff(x,i,i));
    1282      111468 :     dm = b;
    1283             :   }
    1284             : 
    1285     6345912 :   for (i = li-1; i > 0; i--)
    1286             :   {
    1287     4849634 :     GEN diag = gcoeff(x,i,i);
    1288     4849634 :     if (signe(diag) < 0) { gel(x,i) = ZC_neg(gel(x,i)); diag = gcoeff(x,i,i); }
    1289    14417416 :     for (j = i+1; j < li; j++)
    1290             :     {
    1291     9567782 :       b = gcoeff(x,i,j);
    1292     9567782 :       b = center? diviiround(b,diag): truedivii(b, diag);
    1293     9567782 :       if (!signe(b)) continue;
    1294     5511668 :       togglesign(b);
    1295     5511668 :       ZC_lincomb1_inplace(gel(x,j), gel(x,i),b);
    1296     5511668 :       p1 = gel(x,j);
    1297    21075073 :       for (k=1; k<i; k++)
    1298    15563405 :         if (lgefint(gel(p1,k)) > LDM[k]) gel(p1,k) = remii(gel(p1,k), gel(dm,i));
    1299             :     }
    1300     4849634 :     if (gc_needed(av,1))
    1301             :     {
    1302          14 :       if (DEBUGMEM>1) pari_warn(warnmem,"ZM_hnfmod[3]. i=%ld", i);
    1303          14 :       gerepileall(av, 2, &x, &dm); diag = gcoeff(x,i,i);
    1304             :     }
    1305             :   }
    1306     1496278 :   return x;
    1307             : }
    1308             : GEN
    1309     1974005 : ZM_hnfmodall(GEN x, GEN dm, long flag)
    1310             : {
    1311     1974005 :   pari_sp av = avma;
    1312     1974005 :   return gerepilecopy(av, ZM_hnfmodall_i(x, dm, flag));
    1313             : }
    1314             : GEN
    1315      111461 : ZM_hnfmod(GEN x, GEN d) { return ZM_hnfmodall(x,d,0); }
    1316             : GEN
    1317     1657832 : ZM_hnfmodid(GEN x, GEN d)
    1318     1657832 : { return ZM_hnfmodall(x,d,hnf_MODID); }
    1319             : /* return the column echelon form of x with 1's as pivots,
    1320             :  * P contains the row indices containing the pivots in increasing order */
    1321             : static GEN
    1322     1386155 : FpM_echelon(GEN x, GEN *pP, GEN p)
    1323             : {
    1324             :   pari_sp av;
    1325             :   long iP, li, co, i, j, k, def, ldef;
    1326             :   GEN P;
    1327             : 
    1328     1386155 :   co = lg(x); if (co == 1) { *pP = cgetg(1,t_VECSMALL); return cgetg(1,t_MAT); }
    1329     1386155 :   li = lgcols(x);
    1330     1386155 :   iP = 1;
    1331     1386155 :   *pP = P = cgetg(li, t_VECSMALL);
    1332     1386155 :   av = avma;
    1333     1386155 :   x = FpM_red(x, p);
    1334             : 
    1335     1386155 :   ldef = (li > co)? li - co: 0;
    1336     5576908 :   for (def = co-1,i = li-1; i > ldef; i--)
    1337             :   {
    1338     4190753 :     GEN u = NULL;
    1339     7473119 :     for (k = def; k; k--)
    1340             :     {
    1341     5834182 :       u = gcoeff(x,i,k);
    1342     5834182 :       if (signe(u)) break;
    1343             :     }
    1344     4190753 :     if (!k)
    1345             :     {
    1346     1638937 :       if (--ldef < 0) ldef = 0;
    1347     1638937 :       continue;
    1348             :     }
    1349     2551816 :     P[iP++] = i;
    1350     2551816 :     if (k != def) swap(gel(x,def), gel(x,k));
    1351     2551816 :     if (!equali1(u))
    1352     2057131 :       FpV_Fp_mul_part_ip(gel(x,def), Fp_inv(u,p), p, i-1);
    1353     2551816 :     gcoeff(x, i, def) = gen_1;
    1354     9308782 :     for (j = def-1; j; j--)
    1355             :     { /* zero x[i, 1..def-1] using x[i,def] = 1*/
    1356     6756966 :       GEN xj = gel(x,j), u = gel(xj,i);
    1357     6756966 :       if (!signe(u)) continue;
    1358             : 
    1359     3955275 :       ZC_lincomb1_inplace(xj, gel(x,def), negi(u));
    1360     3955275 :       for (k = 1; k < i; k++) gel(xj,k) = modii(gel(xj,k), p);
    1361             :     }
    1362     2551816 :     if (gc_needed(av,2))
    1363             :     {
    1364           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"FpM_echelon. i=%ld",i);
    1365           0 :       x = gerepilecopy(av, x);
    1366             :     }
    1367     2551816 :     def--;
    1368             :   }
    1369             :   /* rank = iP - 1 */
    1370     1386155 :   setlg(P, iP); vecsmall_sort(P);
    1371     1386155 :   if (co > iP) x += co - iP;
    1372     1386155 :   x[0] = evaltyp(t_MAT) | evallg(iP);
    1373     1386155 :   return x;
    1374             : }
    1375             : /* given x square of maximal rank with 1 or p on diagonal from hnfmodid
    1376             :  * (=> a column containing p has its other entries at 0 ), return the HNF */
    1377             : static GEN
    1378     1374189 : FpM_hnfend(pari_sp av, GEN x, GEN p)
    1379             : {
    1380     1374189 :   long i, l = lgcols(x);
    1381     5551999 :   for (i = l-1; i > 0; i--)
    1382             :   {
    1383     4177810 :     GEN diag = gcoeff(x,i,i);
    1384             :     long j;
    1385     4177810 :     if (is_pm1(diag))
    1386     6062149 :       for (j = i+1; j < l; j++)
    1387             :       {
    1388     3533608 :         GEN xj = gel(x,j), b = gel(xj,i);
    1389             :         long k;
    1390     3533608 :         if (!signe(b)) continue;
    1391     2423183 :         ZC_lincomb1_inplace(xj, gel(x,i), negi(b));
    1392    11315030 :         for (k=1; k<i; k++)
    1393     8891847 :           if (lgefint(gel(xj,k)) > 3) gel(xj,k) = remii(gel(xj,k), p);
    1394             :       }
    1395             :     else
    1396     1649269 :       for (j = i+1; j < l; j++) gcoeff(x,i,j) = modii(gcoeff(x,i,j), p);
    1397     4177810 :     if (gc_needed(av,2))
    1398             :     {
    1399           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"FpM_hnfend. i=%ld",i);
    1400           0 :       x = gerepilecopy(av, x);
    1401             :     }
    1402             :   }
    1403     1374189 :   return x;
    1404             : }
    1405             : GEN
    1406     1386155 : ZM_hnfmodprime(GEN x, GEN p)
    1407             : {
    1408     1386155 :   pari_sp av = avma;
    1409             :   GEN P, y;
    1410             :   long l, lP, i;
    1411     1386155 :   if (lg(x) == 1) return cgetg(1, t_MAT);
    1412     1386155 :   l = lgcols(x);
    1413     1386155 :   x = FpM_echelon(x, &P, p);
    1414     1386155 :   lP = lg(P); /* rank = lP-1 */
    1415     1386155 :   if (lP == l) { avma = av; return matid(l-1); }
    1416     1374189 :   y = scalarmat_shallow(p, l-1);
    1417     1374189 :   for (i = 1; i < lP; i++) gel(y,P[i]) = gel(x,i);
    1418     1374189 :   return gerepilecopy(av, FpM_hnfend(av,y,p));
    1419             : }
    1420             : 
    1421             : static GEN
    1422      204047 : allhnfmod(GEN x, GEN dm, int flag)
    1423             : {
    1424      204047 :   if (typ(x)!=t_MAT) pari_err_TYPE("allhnfmod",x);
    1425      204047 :   RgM_check_ZM(x, "allhnfmod");
    1426      204047 :   if (isintzero(dm)) return ZM_hnf(x);
    1427      204047 :   return ZM_hnfmodall(x, dm, flag);
    1428             : }
    1429             : GEN
    1430          14 : hnfmod(GEN x, GEN d)
    1431             : {
    1432          14 :   if (typ(d) != t_INT) pari_err_TYPE("mathnfmod",d);
    1433          14 :   return allhnfmod(x, d, 0);
    1434             : }
    1435             : GEN
    1436      204033 : hnfmodid(GEN x, GEN d)
    1437             : {
    1438      204033 :   switch(typ(d))
    1439             :   {
    1440      204019 :     case t_INT: break;
    1441             :     case t_VEC: case t_COL:
    1442          14 :       if (RgV_is_ZV(d)) break;
    1443           0 :     default: pari_err_TYPE("mathnfmodid",d);
    1444             :   }
    1445      204033 :   return allhnfmod(x, d, hnf_MODID);
    1446             : }
    1447             : 
    1448             : /* M a ZM in HNF. Normalize with *centered* residues */
    1449             : GEN
    1450        1267 : ZM_hnfcenter(GEN M)
    1451             : {
    1452        1267 :   long i, j, k, N = lg(M)-1;
    1453        1267 :   pari_sp av = avma;
    1454             : 
    1455        5271 :   for (j=N-1; j>0; j--) /* skip last line */
    1456             :   {
    1457        4004 :     GEN Mj = gel(M,j), a = gel(Mj,j);
    1458       23135 :     for (k = j+1; k <= N; k++)
    1459             :     {
    1460       19131 :       GEN Mk = gel(M,k), q = diviiround(gel(Mk,j), a);
    1461       19131 :       long s = signe(q);
    1462       19131 :       if (!s) continue;
    1463       12376 :       if (is_pm1(q))
    1464             :       {
    1465        2247 :         if (s < 0)
    1466         357 :           for (i = 1; i <= j; i++) gel(Mk,i) = addii(gel(Mk,i), gel(Mj,i));
    1467             :         else
    1468        1890 :           for (i = 1; i <= j; i++) gel(Mk,i) = subii(gel(Mk,i), gel(Mj,i));
    1469             :       }
    1470             :       else
    1471       10129 :         for (i = 1; i <= j; i++) gel(Mk,i) = subii(gel(Mk,i), mulii(q,gel(Mj,i)));
    1472       12376 :       if (gc_needed(av,1))
    1473             :       {
    1474           0 :         if (DEBUGMEM>1) pari_warn(warnmem,"ZM_hnfcenter, j = %ld",j);
    1475           0 :         M = gerepilecopy(av, M);
    1476             :       }
    1477             :     }
    1478             :   }
    1479        1267 :   return M;
    1480             : }
    1481             : 
    1482             : /***********************************************************************/
    1483             : /*                                                                     */
    1484             : /*                 HNFLLL (Havas, Majewski, Mathews)                   */
    1485             : /*                                                                     */
    1486             : /***********************************************************************/
    1487             : 
    1488             : static void
    1489      505817 : Minus(long j, GEN lambda)
    1490             : {
    1491      505817 :   long k, n = lg(lambda);
    1492             : 
    1493      505817 :   for (k=1  ; k<j; k++) togglesign_safe(&gcoeff(lambda,k,j));
    1494      505817 :   for (k=j+1; k<n; k++) togglesign_safe(&gcoeff(lambda,j,k));
    1495      505817 : }
    1496             : 
    1497             : /* index of first non-zero entry */
    1498             : static long
    1499    36207301 : findi(GEN M)
    1500             : {
    1501    36207301 :   long i, n = lg(M);
    1502   242377464 :   for (i=1; i<n; i++)
    1503   221805252 :     if (signe(gel(M,i))) return i;
    1504    20572212 :   return 0;
    1505             : }
    1506             : 
    1507             : static long
    1508    36207301 : findi_normalize(GEN Aj, GEN B, long j, GEN lambda)
    1509             : {
    1510    36207301 :   long r = findi(Aj);
    1511    36207301 :   if (r && signe(gel(Aj,r)) < 0)
    1512             :   {
    1513      505817 :     ZV_togglesign(Aj); if (B) ZV_togglesign(gel(B,j));
    1514      505817 :     Minus(j,lambda);
    1515             :   }
    1516    36207301 :   return r;
    1517             : }
    1518             : 
    1519             : static void
    1520    18103196 : reduce2(GEN A, GEN B, long k, long j, long *row0, long *row1, GEN lambda, GEN D)
    1521             : {
    1522             :   GEN q;
    1523             :   long i;
    1524             : 
    1525    18103196 :   *row0 = findi_normalize(gel(A,j), B,j,lambda);
    1526    18103196 :   *row1 = findi_normalize(gel(A,k), B,k,lambda);
    1527    18103196 :   if (*row0)
    1528     5317453 :     q = truedivii(gcoeff(A,*row0,k), gcoeff(A,*row0,j));
    1529    12785743 :   else if (abscmpii(shifti(gcoeff(lambda,j,k), 1), gel(D,j)) > 0)
    1530     2608276 :     q = diviiround(gcoeff(lambda,j,k), gel(D,j));
    1531             :   else
    1532    28280663 :     return;
    1533             : 
    1534     7925729 :   if (signe(q))
    1535             :   {
    1536     4728960 :     GEN Lk = gel(lambda,k), Lj = gel(lambda,j);
    1537     4728960 :     togglesign_safe(&q);
    1538     4728960 :     if (*row0) ZC_lincomb1_inplace(gel(A,k),gel(A,j),q);
    1539     4728960 :     if (B) ZC_lincomb1_inplace(gel(B,k),gel(B,j),q);
    1540     4728960 :     gel(Lk,j) = addii(gel(Lk,j), mulii(q,gel(D,j)));
    1541     4728960 :     if (is_pm1(q))
    1542             :     {
    1543     1480263 :       if (signe(q) > 0)
    1544             :       {
    1545     5129356 :         for (i=1; i<j; i++)
    1546     4565062 :           if (signe(gel(Lj,i))) gel(Lk,i) = addii(gel(Lk,i), gel(Lj,i));
    1547             :       }
    1548             :       else
    1549             :       {
    1550     6217254 :         for (i=1; i<j; i++)
    1551     5301285 :           if (signe(gel(Lj,i))) gel(Lk,i) = subii(gel(Lk,i), gel(Lj,i));
    1552             :       }
    1553             :     }
    1554             :     else
    1555             :     {
    1556    18990223 :       for (i=1; i<j; i++)
    1557    15741526 :         if (signe(gel(Lj,i))) gel(Lk,i) = addii(gel(Lk,i), mulii(q,gel(Lj,i)));
    1558             :     }
    1559             :   }
    1560             : }
    1561             : 
    1562             : static void
    1563     2445462 : hnfswap(GEN A, GEN B, long k, GEN lambda, GEN D)
    1564             : {
    1565     2445462 :   GEN t, p1, p2, Lk = gel(lambda,k);
    1566     2445462 :   long i,j,n = lg(A);
    1567             : 
    1568     2445462 :   swap(gel(A,k), gel(A,k-1));
    1569     2445462 :   if (B) swap(gel(B,k), gel(B,k-1));
    1570     2445462 :   for (j=k-2; j; j--) swap(gcoeff(lambda,j,k-1), gel(Lk,j));
    1571    31145260 :   for (i=k+1; i<n; i++)
    1572             :   {
    1573    28699798 :     GEN Li = gel(lambda,i);
    1574    28699798 :     p1 = mulii(gel(Li,k-1), gel(D,k));
    1575    28699798 :     p2 = mulii(gel(Li,k), gel(Lk,k-1));
    1576    28699798 :     t = subii(p1,p2);
    1577             : 
    1578    28699798 :     p1 = mulii(gel(Li,k), gel(D,k-2));
    1579    28699798 :     p2 = mulii(gel(Li,k-1), gel(Lk,k-1));
    1580    28699798 :     gel(Li,k-1) = diviiexact(addii(p1,p2), gel(D,k-1));
    1581    28699798 :     gel(Li,k) = diviiexact(t, gel(D,k-1));
    1582             :   }
    1583     2445462 :   p1 = mulii(gel(D,k-2), gel(D,k));
    1584     2445462 :   p2 = sqri(gel(Lk,k-1));
    1585     2445462 :   gel(D,k-1) = diviiexact(addii(p1,p2), gel(D,k-1));
    1586     2445462 : }
    1587             : 
    1588             : /* reverse row order in matrix A, IN PLACE */
    1589             : static GEN
    1590      101086 : reverse_rows(GEN A)
    1591             : {
    1592      101086 :   long i, j, h, n = lg(A);
    1593      101086 :   if (n == 1) return A;
    1594      101030 :   h = lgcols(A);
    1595      873212 :   for (j=1; j<n; j++)
    1596             :   {
    1597      772182 :     GEN c = gel(A,j);
    1598             :     /* start at (h-1) >>1 : if h = 2i even, no need to swap c[i] and itself */
    1599      772182 :     for (i=(h-1)>>1; i; i--) swap(gel(c,i), gel(c,h-i));
    1600             :   }
    1601      101030 :   return A;
    1602             : }
    1603             : 
    1604             : GEN
    1605       50543 : ZM_hnflll(GEN A, GEN *ptB, int remove)
    1606             : {
    1607       50543 :   pari_sp av = avma;
    1608             : #ifdef HNFLLL_QUALITY
    1609             :   const long m1 = 1, n1 = 1; /* alpha = m1/n1. Maybe 3/4 here ? */
    1610             : #endif
    1611             :   long n, k, kmax;
    1612             :   GEN B, lambda, D;
    1613             : 
    1614       50543 :   n = lg(A);
    1615       50543 :   A = reverse_rows(ZM_copy(A)); /* ZM_copy for in place findi_normalize() */
    1616       50543 :   B = ptB? matid(n-1): NULL;
    1617       50543 :   D = const_vec(n, gen_1) + 1;
    1618       50543 :   lambda = zeromatcopy(n-1,n-1);
    1619       50543 :   k = kmax = 2;
    1620     4346496 :   while (k < n)
    1621             :   {
    1622             :     long row0, row1;
    1623             :     int do_swap;
    1624     4245410 :     reduce2(A,B,k,k-1,&row0,&row1,lambda,D);
    1625     4245410 :     if (row0)
    1626     2678553 :       do_swap = (!row1 || row0 <= row1);
    1627     1566857 :     else if (!row1)
    1628             :     { /* row0 == row1 == 0 */
    1629     1361326 :       pari_sp av1 = avma;
    1630     1361326 :       GEN z = addii(mulii(gel(D,k-2),gel(D,k)), sqri(gcoeff(lambda,k-1,k)));
    1631             : #ifdef HNFLLL_QUALITY
    1632             :       do_swap = (cmpii(mului(n1,z), mului(m1,sqri(gel(D,k-1)))) < 0);
    1633             : #else /* assume m1 = n1 = 1 */
    1634     1361326 :       do_swap = (cmpii(z, sqri(gel(D,k-1))) < 0);
    1635             : #endif
    1636     1361326 :       avma = av1;
    1637             :     }
    1638             :     else
    1639      205531 :       do_swap = 0;
    1640     4245410 :     if (do_swap)
    1641             :     {
    1642     2384578 :       hnfswap(A,B,k,lambda,D);
    1643     2384578 :       if (k > 2) k--;
    1644             :     }
    1645             :     else
    1646             :     {
    1647             :       long i;
    1648    15718618 :       for (i=k-2; i; i--)
    1649             :       {
    1650             :         long row0, row1;
    1651    13857786 :         reduce2(A,B,k,i,&row0,&row1,lambda,D);
    1652    13857786 :         if (gc_needed(av,3))
    1653             :         {
    1654          37 :           GEN b = D-1;
    1655          37 :           if (DEBUGMEM>1) pari_warn(warnmem,"hnflll (reducing), kmax = %ld",kmax);
    1656          37 :           gerepileall(av, B? 4: 3, &A, &lambda, &b, &B);
    1657          37 :           D = b+1;
    1658             :         }
    1659             :       }
    1660     1860832 :       if (++k > kmax) kmax = k;
    1661             :     }
    1662     4245410 :     if (gc_needed(av,3))
    1663             :     {
    1664          31 :       GEN b = D-1;
    1665          31 :       if (DEBUGMEM>1) pari_warn(warnmem,"hnflll, kmax = %ld / %ld",kmax,n-1);
    1666          31 :       gerepileall(av, B? 4: 3, &A, &lambda, &b, &B);
    1667          31 :       D = b+1;
    1668             :     }
    1669             :   }
    1670             :   /* handle trivial case: return negative diag coefficient otherwise */
    1671       50543 :   if (n == 2) (void)findi_normalize(gel(A,1), B,1,lambda);
    1672       50543 :   A = reverse_rows(A);
    1673       50543 :   if (remove)
    1674             :   {
    1675             :     long i;
    1676       59200 :     for (i = 1; i < n; i++)
    1677       53999 :       if (!ZV_equal0(gel(A,i))) break;
    1678        9042 :     remove_0cols(i-1, &A, &B, remove);
    1679             :   }
    1680       50543 :   gerepileall(av, B? 2: 1, &A, &B);
    1681       50543 :   if (B) *ptB = B;
    1682       50543 :   return A;
    1683             : }
    1684             : 
    1685             : GEN
    1686           7 : hnflll(GEN x)
    1687             : {
    1688           7 :   GEN z = cgetg(3, t_VEC);
    1689           7 :   gel(z,1) = ZM_hnflll(x, &gel(z,2), 1);
    1690           7 :   return z;
    1691             : }
    1692             : 
    1693             : /* Variation on HNFLLL: Extended GCD */
    1694             : 
    1695             : static void
    1696      368158 : reduce1(GEN A, GEN B, long k, long j, GEN lambda, GEN D)
    1697             : {
    1698             :   GEN q;
    1699             :   long i;
    1700             : 
    1701      368158 :   if (signe(gel(A,j)))
    1702       55239 :     q = diviiround(gel(A,k),gel(A,j));
    1703      312919 :   else if (abscmpii(shifti(gcoeff(lambda,j,k), 1), gel(D,j)) > 0)
    1704        8353 :     q = diviiround(gcoeff(lambda,j,k), gel(D,j));
    1705             :   else
    1706      672724 :     return;
    1707             : 
    1708       63592 :   if (signe(q))
    1709             :   {
    1710       49503 :     GEN Lk = gel(lambda,k), Lj = gel(lambda,j);
    1711       49503 :     togglesign_safe(&q);
    1712       49503 :     gel(A,k) = addii(gel(A,k), mulii(q,gel(A,j)));
    1713       49503 :     ZC_lincomb1_inplace(gel(B,k),gel(B,j),q);
    1714       49503 :     gel(Lk,j) = addii(gel(Lk,j), mulii(q,gel(D,j)));
    1715       73498 :     for (i=1; i<j; i++)
    1716       23995 :       if (signe(gel(Lj,i))) gel(Lk,i) = addii(gel(Lk,i), mulii(q,gel(Lj,i)));
    1717             :   }
    1718             : }
    1719             : 
    1720             : static GEN
    1721       50932 : ZV_gcdext_i(GEN A)
    1722             : {
    1723             : #ifdef HNFLLL_QUALITY
    1724             :   const long m1 = 1, n1 = 1; /* alpha = m1/n1. Maybe 3/4 here ? */
    1725             : #endif
    1726       50932 :   long k, n = lg(A);
    1727             :   GEN B, lambda, D;
    1728             : 
    1729       50932 :   if (n == 1) retmkvec2(gen_1, cgetg(1,t_MAT));
    1730       50932 :   A = leafcopy(A);
    1731       50932 :   B = matid(n-1);
    1732       50932 :   lambda = zeromatcopy(n-1,n-1);
    1733       50932 :   D = const_vec(n, gen_1) + 1;
    1734       50932 :   k = 2;
    1735      291168 :   while (k < n)
    1736             :   {
    1737             :     int do_swap;
    1738             : 
    1739      189304 :     reduce1(A,B,k,k-1,lambda,D);
    1740      189304 :     if (signe(gel(A,k-1))) do_swap = 1;
    1741      134065 :     else if (!signe(gel(A,k)))
    1742             :     {
    1743       60836 :       pari_sp av1 = avma;
    1744       60836 :       GEN z = addii(mulii(gel(D,k-2),gel(D,k)), sqri(gcoeff(lambda,k-1,k)));
    1745             : #ifdef HNFLLL_QUALITY
    1746             :       do_swap = (cmpii(mului(n1,z), mului(m1,sqri(gel(D,k-1)))) < 0);
    1747             : #else /* assume m1 = n1 = 1 */
    1748       60836 :       do_swap = (cmpii(z, sqri(gel(D,k-1))) < 0);
    1749             : #endif
    1750       60836 :       avma = av1;
    1751             :     }
    1752       73229 :     else do_swap = 0;
    1753             : 
    1754      189304 :     if (do_swap)
    1755             :     {
    1756       60884 :       hnfswap(A,B,k,lambda,D);
    1757       60884 :       if (k > 2) k--;
    1758             :     }
    1759             :     else
    1760             :     {
    1761             :       long i;
    1762      128420 :       for (i=k-2; i; i--) reduce1(A,B,k,i,lambda,D);
    1763      128420 :       k++;
    1764             :     }
    1765             :   }
    1766       50932 :   if (signe(gel(A,n-1)) < 0)
    1767             :   {
    1768        6682 :     gel(A,n-1) = negi(gel(A,n-1));
    1769        6682 :     ZV_togglesign(gel(B,n-1));
    1770             :   }
    1771       50932 :   return mkvec2(gel(A,n-1), B);
    1772             : }
    1773             : GEN
    1774       50918 : ZV_extgcd(GEN A)
    1775             : {
    1776       50918 :   pari_sp av = avma;
    1777       50918 :   return gerepilecopy(av, ZV_gcdext_i(A));
    1778             : }
    1779             : /* as ZV_extgcd, transforming the gcd into a t_MAT, for mathnf0 */
    1780             : static GEN
    1781          21 : ZV_hnfgcdext(GEN A)
    1782             : {
    1783          21 :   pari_sp av = avma;
    1784             :   GEN z;
    1785          21 :   if (lg(A) == 1) retmkvec2(cgetg(1,t_MAT),cgetg(1,t_MAT));
    1786          14 :   z = ZV_gcdext_i(A);
    1787          14 :   gel(z,1) = mkmat(mkcol(gel(z,1)));
    1788          14 :   return gerepilecopy(av, z);
    1789             : }
    1790             : 
    1791             : /* HNF with permutation. */
    1792             : GEN
    1793         238 : ZM_hnfperm(GEN A, GEN *ptU, GEN *ptperm)
    1794             : {
    1795             :   GEN U, c, l, perm, d, p, q, b;
    1796         238 :   pari_sp av = avma, av1;
    1797             :   long r, t, i, j, j1, k, m, n;
    1798             : 
    1799         238 :   n = lg(A)-1;
    1800         238 :   if (!n)
    1801             :   {
    1802           7 :     if (ptU) *ptU = cgetg(1,t_MAT);
    1803           7 :     if (ptperm) *ptperm = cgetg(1,t_VEC);
    1804           7 :     return cgetg(1, t_MAT);
    1805             :   }
    1806         231 :   m = nbrows(A);
    1807         231 :   c = zero_zv(m);
    1808         231 :   l = zero_zv(n);
    1809         231 :   perm = cgetg(m+1, t_VECSMALL);
    1810         231 :   av1 = avma;
    1811         231 :   A = RgM_shallowcopy(A);
    1812         231 :   U = ptU? matid(n): NULL;
    1813             :   /* U base change matrix : A0*U = A all along */
    1814        1148 :   for (r=0, k=1; k <= n; k++)
    1815             :   {
    1816        2849 :     for (j=1; j<k; j++)
    1817             :     {
    1818        1932 :       if (!l[j]) continue;
    1819        1799 :       t=l[j]; b=gcoeff(A,t,k);
    1820        1799 :       if (!signe(b)) continue;
    1821             : 
    1822         525 :       ZC_elem(b,gcoeff(A,t,j), A,U,k,j);
    1823         525 :       d = gcoeff(A,t,j);
    1824         525 :       if (signe(d) < 0)
    1825             :       {
    1826           0 :         ZV_neg_inplace(gel(A,j));
    1827           0 :         if (U) ZV_togglesign(gel(U,j));
    1828           0 :         d = gcoeff(A,t,j);
    1829             :       }
    1830        1253 :       for (j1=1; j1<j; j1++)
    1831             :       {
    1832         728 :         if (!l[j1]) continue;
    1833         714 :         q = truedivii(gcoeff(A,t,j1),d);
    1834         714 :         if (!signe(q)) continue;
    1835             : 
    1836         329 :         togglesign(q);
    1837         329 :         ZC_lincomb1_inplace(gel(A,j1), gel(A,j), q);
    1838         329 :         if (U) ZC_lincomb1_inplace(gel(U,j1), gel(U,j), q);
    1839             :       }
    1840             :     }
    1841         917 :     t = m; while (t && (c[t] || !signe(gcoeff(A,t,k)))) t--;
    1842         917 :     if (t)
    1843             :     {
    1844         819 :       p = gcoeff(A,t,k);
    1845       13846 :       for (i=t-1; i; i--)
    1846             :       {
    1847       13027 :         q = gcoeff(A,i,k);
    1848       13027 :         if (signe(q) && abscmpii(p,q) > 0) { p = q; t = i; }
    1849             :       }
    1850         819 :       perm[++r] = l[k] = t; c[t] = k;
    1851         819 :       if (signe(p) < 0)
    1852             :       {
    1853          97 :         ZV_neg_inplace(gel(A,k));
    1854          97 :         if (U) ZV_togglesign(gel(U,k));
    1855          97 :         p = gcoeff(A,t,k);
    1856             :       }
    1857             :       /* p > 0 */
    1858        2408 :       for (j=1; j<k; j++)
    1859             :       {
    1860        1589 :         if (!l[j]) continue;
    1861        1568 :         q = truedivii(gcoeff(A,t,j),p);
    1862        1568 :         if (!signe(q)) continue;
    1863             : 
    1864         244 :         togglesign(q);
    1865         244 :         ZC_lincomb1_inplace(gel(A,j), gel(A,k), q);
    1866         244 :         if (U) ZC_lincomb1_inplace(gel(U,j), gel(U,k), q);
    1867             :       }
    1868             :     }
    1869         917 :     if (gc_needed(av1,1))
    1870             :     {
    1871           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"hnfperm, k=%ld",k);
    1872           0 :       gerepileall(av1, U? 2: 1, &A, &U);
    1873             :     }
    1874             :   }
    1875         231 :   if (r < m)
    1876             :   {
    1877        3227 :     for (i=1,k=r; i<=m; i++)
    1878        3052 :       if (!c[i]) perm[++k] = i;
    1879             :   }
    1880             : 
    1881             : /* We have A0*U=A, U in Gl(n,Z)
    1882             :  * basis for Im(A):  columns of A s.t l[j]>0 (r   cols)
    1883             :  * basis for Ker(A): columns of U s.t l[j]=0 (n-r cols) */
    1884         231 :   p = cgetg(r+1,t_MAT);
    1885         231 :   for (i=1; i<=m/2; i++) lswap(perm[i], perm[m+1-i]);
    1886         231 :   if (U)
    1887             :   {
    1888          70 :     GEN u = cgetg(n+1,t_MAT);
    1889         322 :     for (t=1,k=r,j=1; j<=n; j++)
    1890         252 :       if (l[j])
    1891             :       {
    1892         154 :         u[k + n-r] = U[j];
    1893         154 :         gel(p,k--) = vecpermute(gel(A,j), perm);
    1894             :       }
    1895             :       else
    1896          98 :         u[t++] = U[j];
    1897          70 :     *ptU = u;
    1898          70 :     if (ptperm) *ptperm = perm;
    1899          70 :     gerepileall(av, ptperm? 3: 2, &p, ptU, ptperm);
    1900             :   }
    1901             :   else
    1902             :   {
    1903         826 :     for (k=r,j=1; j<=n; j++)
    1904         665 :       if (l[j]) gel(p,k--) = vecpermute(gel(A,j), perm);
    1905         161 :     if (ptperm) *ptperm = perm;
    1906         161 :     gerepileall(av, ptperm? 2: 1, &p, ptperm);
    1907             :   }
    1908         231 :   return p;
    1909             : }
    1910             : 
    1911             : GEN
    1912         154 : ZM_hnf_knapsack(GEN x)
    1913             : {
    1914         154 :   GEN t, perm, H = ZM_hnfperm(x,NULL,&perm);
    1915         154 :   long i,j, l = lg(H), h = lgcols(H);
    1916        1743 :   for (i=1; i<h; i++)
    1917             :   {
    1918        1673 :     int fl = 0;
    1919       10542 :     for (j=1; j<l; j++)
    1920             :     {
    1921        8953 :       t = gcoeff(H,i,j);
    1922        8953 :       if (signe(t))
    1923             :       {
    1924        1701 :         if (!is_pm1(t) || fl) return NULL;
    1925        1617 :         fl = 1;
    1926             :       }
    1927             :     }
    1928             :   }
    1929          70 :   return rowpermute(H, perm_inv(perm));
    1930             : }
    1931             : 
    1932             : GEN
    1933          14 : hnfperm(GEN A)
    1934             : {
    1935          14 :   GEN y = cgetg(4, t_VEC);
    1936          14 :   gel(y,1) = ZM_hnfperm(A, &gel(y,2), &gel(y,3));
    1937          14 :   return y;
    1938             : }
    1939             : 
    1940             : /* Hermite Normal Form, with base change matrix if ptB != NULL.
    1941             :  * If 'remove' = 1, remove 0 columns (do NOT update *ptB accordingly)
    1942             :  * If 'remove' = 2, remove 0 columns and update *ptB accordingly */
    1943             : GEN
    1944      159347 : ZM_hnfall_i(GEN A, GEN *ptB, long remove)
    1945             : {
    1946             :   pari_sp av;
    1947             :   long m, n, r, i, j, k, li;
    1948             :   GEN B, c, h, a;
    1949             : 
    1950      159347 :   RgM_dimensions(A, &m,&n);
    1951      159347 :   if (!n)
    1952             :   {
    1953           7 :     if (ptB) *ptB = cgetg(1,t_MAT);
    1954           7 :     return cgetg(1,t_MAT);
    1955             :   }
    1956      159340 :   c = zero_zv(m);
    1957      159340 :   h = const_vecsmall(n, m);
    1958      159340 :   av = avma;
    1959      159340 :   A = RgM_shallowcopy(A);
    1960      159340 :   B = ptB? matid(n): NULL;
    1961      159340 :   r = n+1;
    1962      564183 :   for (li=m; li; li--)
    1963             :   {
    1964      997599 :     for (j=1; j<r; j++)
    1965             :     {
    1966     1420201 :       for (i=h[j]; i>li; i--)
    1967             :       {
    1968      422763 :         a = gcoeff(A,i,j);
    1969      422763 :         k = c[i];
    1970             :         /* zero a = Aij  using  Aik */
    1971      422763 :         if (signe(a)) ZC_elem(a,gcoeff(A,i,k), A,B,j,k);
    1972      422763 :         ZM_reduce(A,B, i,k); /* ensure reduced entries */
    1973      422763 :         if (gc_needed(av,1))
    1974             :         {
    1975           0 :           if (DEBUGMEM>1)
    1976           0 :             pari_warn(warnmem,"ZM_hnfall[1], li = %ld, j = %ld", li, j);
    1977           0 :           gerepileall(av, B? 2: 1, &A, &B);
    1978             :         }
    1979             :       }
    1980      997438 :       if (signe( gcoeff(A,li,j) )) break;
    1981      592756 :       h[j] = li-1;
    1982             :     }
    1983      404843 :     if (j == r) continue;
    1984      404682 :     r--;
    1985      404682 :     if (j < r) /* A[j] != 0 */
    1986             :     {
    1987      328940 :       swap(gel(A,j), gel(A,r));
    1988      328940 :       if (B) swap(gel(B,j), gel(B,r));
    1989      328940 :       h[j] = h[r]; h[r] = li; c[li] = r;
    1990             :     }
    1991      404682 :     if (signe(gcoeff(A,li,r)) < 0)
    1992             :     {
    1993      117570 :       ZV_neg_inplace(gel(A,r));
    1994      117570 :       if (B) ZV_togglesign(gel(B,r));
    1995             :     }
    1996      404682 :     ZM_reduce(A,B, li,r);
    1997      404682 :     if (gc_needed(av,1))
    1998             :     {
    1999           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"ZM_hnfall[2], li = %ld", li);
    2000           0 :       gerepileall(av, B? 2: 1, &A, &B);
    2001             :     }
    2002             :   }
    2003             : 
    2004      159340 :   if (DEBUGLEVEL>5) err_printf("\nhnfall, final phase: ");
    2005      159340 :   r--; /* first r cols are in the image the n-r (independent) last ones */
    2006      852406 :   for (j=1; j<=r; j++)
    2007     1836358 :     for (i=h[j]; i; i--)
    2008             :     {
    2009     1143292 :       a = gcoeff(A,i,j);
    2010     1143292 :       k = c[i];
    2011     1143292 :       if (signe(a)) ZC_elem(a,gcoeff(A,i,k), A,B, j,k);
    2012     1143292 :       ZM_reduce(A,B, i,k); /* ensure reduced entries, even if a = 0 */
    2013     1143292 :       if (gc_needed(av,1))
    2014             :       {
    2015           0 :         if (DEBUGMEM>1) pari_warn(warnmem,"ZM_hnfall[3], j = %ld", j);
    2016           0 :         gerepileall(av, B? 2: 1, &A, &B);
    2017             :       }
    2018             :     }
    2019      159340 :   if (DEBUGLEVEL>5) err_printf("\n");
    2020      159340 :   if (remove) remove_0cols(r, &A, &B, remove);
    2021      159340 :   if (ptB) *ptB = B;
    2022      159340 :   return A;
    2023             : }
    2024             : GEN
    2025        2812 : ZM_hnfall(GEN A, GEN *ptB, long remove)
    2026             : {
    2027        2812 :   pari_sp av = avma;
    2028        2812 :   A = ZM_hnfall_i(A, ptB, remove);
    2029        2812 :   gerepileall(av, ptB? 2: 1, &A, ptB);
    2030        2812 :   return A;
    2031             : }
    2032             : 
    2033             : GEN
    2034          14 : hnfall(GEN x)
    2035             : {
    2036          14 :   GEN z = cgetg(3, t_VEC);
    2037          14 :   gel(z,1) = ZM_hnfall(x, (GEN*)(z+2), 1);
    2038          14 :   return z;
    2039             : }
    2040             : GEN
    2041          14 : hnf(GEN x) { return mathnf0(x,0); }
    2042             : 
    2043             : /* C = A^(-1)t where A and C are integral, A is upper triangular, t t_INT */
    2044             : GEN
    2045       20797 : hnf_invscale(GEN A, GEN t)
    2046             : {
    2047       20797 :   long n = lg(A)-1, i,j,k;
    2048       20797 :   GEN m, c = cgetg(n+1,t_MAT);
    2049             : 
    2050       20797 :   if (!n) return c;
    2051      114982 :   for (k=1; k<=n; k++)
    2052             :   { /* cf hnf_divscale with B = id, thus b = e_k */
    2053       94185 :     GEN u = cgetg(n+1, t_COL);
    2054       94185 :     pari_sp av = avma;
    2055       94185 :     gel(c,k) = u;
    2056       94185 :     gel(u,n) = k == n? gerepileuptoint(av, diviiexact(t, gcoeff(A,n,n))): gen_0;
    2057      621453 :     for (i=n-1; i>0; i--)
    2058             :     {
    2059      527268 :       av = avma; m = i == k? t: gen_0;
    2060      527268 :       for (j=i+1; j<=n; j++) m = subii(m, mulii(gcoeff(A,i,j),gel(u,j)));
    2061      527268 :       gel(u,i) = gerepileuptoint(av, diviiexact(m, gcoeff(A,i,i)));
    2062             :     }
    2063             :   }
    2064       20797 :   return c;
    2065             : }
    2066             : 
    2067             : /* C = A^(-1)(tB) where A, B, C are integral, A is upper triangular, t t_INT */
    2068             : GEN
    2069      180490 : hnf_divscale(GEN A, GEN B, GEN t)
    2070             : {
    2071      180490 :   long n = lg(A)-1, i,j,k;
    2072      180490 :   GEN m, c = cgetg(n+1,t_MAT);
    2073             : 
    2074      180490 :   if (!n) return c;
    2075      707293 :   for (k=1; k<=n; k++)
    2076             :   {
    2077      526803 :     GEN u = cgetg(n+1, t_COL), b = gel(B,k);
    2078      526803 :     pari_sp av = avma;
    2079      526803 :     gel(c,k) = u; m = mulii(gel(b,n),t);
    2080      526803 :     gel(u,n) = gerepileuptoint(av, diviiexact(m, gcoeff(A,n,n)));
    2081     2257787 :     for (i=n-1; i>0; i--)
    2082             :     {
    2083     1730984 :       av = avma; m = mulii(gel(b,i),t);
    2084     1730984 :       for (j=i+1; j<=n; j++) m = subii(m, mulii(gcoeff(A,i,j),gel(u,j)));
    2085     1730984 :       gel(u,i) = gerepileuptoint(av, diviiexact(m, gcoeff(A,i,i)));
    2086             :     }
    2087             :   }
    2088      180490 :   return c;
    2089             : }
    2090             : 
    2091             : /* A, B integral upper HNF. A^(-1) B integral ? */
    2092             : int
    2093          91 : hnfdivide(GEN A, GEN B)
    2094             : {
    2095          91 :   pari_sp av = avma;
    2096          91 :   long n = lg(A)-1, i,j,k;
    2097             :   GEN u, b, m, r;
    2098             : 
    2099          91 :   if (!n) return 1;
    2100          91 :   if (lg(B)-1 != n) pari_err_DIM("hnfdivide");
    2101          91 :   u = cgetg(n+1, t_COL);
    2102         413 :   for (k=1; k<=n; k++)
    2103             :   {
    2104         322 :     b = gel(B,k);
    2105         322 :     m = gel(b,k);
    2106         322 :     gel(u,k) = dvmdii(m, gcoeff(A,k,k), &r);
    2107         322 :     if (r != gen_0) { avma = av; return 0; }
    2108         798 :     for (i=k-1; i>0; i--)
    2109             :     {
    2110         476 :       m = gel(b,i);
    2111         476 :       for (j=i+1; j<=k; j++) m = subii(m, mulii(gcoeff(A,i,j),gel(u,j)));
    2112         476 :       m = dvmdii(m, gcoeff(A,i,i), &r);
    2113         476 :       if (r != gen_0) { avma = av; return 0; }
    2114         476 :       gel(u,i) = m;
    2115             :     }
    2116             :   }
    2117          91 :   avma = av; return 1;
    2118             : }
    2119             : 
    2120             : /* A upper HNF, b integral vector. Return A^(-1) b if integral,
    2121             :  * NULL otherwise. Assume #A[,1] = #b. */
    2122             : GEN
    2123      113211 : hnf_invimage(GEN A, GEN b)
    2124             : {
    2125      113211 :   pari_sp av = avma;
    2126      113211 :   long n = lg(A)-1, m, i, k;
    2127             :   GEN u, r;
    2128             : 
    2129      113211 :   if (!n) return lg(b)==1? cgetg(1,t_COL):NULL;
    2130      113169 :   m = nbrows(A); /* m >= n */
    2131      113169 :   u = cgetg(n+1, t_COL);
    2132      348313 :   for (i = n, k = m; k > 0; k--)
    2133             :   {
    2134      348313 :     pari_sp av2 = avma;
    2135             :     long j;
    2136      348313 :     GEN t = gel(b,k), Aki = gcoeff(A,k,i);
    2137      348313 :     if (typ(t) != t_INT) pari_err_TYPE("hnf_invimage",t);
    2138      348313 :     for (j=i+1; j<=n; j++) t = subii(t, mulii(gcoeff(A,k,j),gel(u,j)));
    2139      348313 :     if (!signe(Aki))
    2140             :     {
    2141           0 :       if (signe(t)) { avma = av;return NULL; }
    2142           0 :       avma = av2; gel(u,i) = gen_0; continue;
    2143             :     }
    2144      348313 :     t = dvmdii(t, Aki, &r);
    2145      348313 :     if (r != gen_0) { avma = av; return NULL; }
    2146      315308 :     gel(u,i) = gerepileuptoint(av2, t);
    2147      315308 :     if (--i == 0) break;
    2148             :   }
    2149             :   /* If there is a solution, it must be u. Check remaining equations */
    2150      160328 :   for (; k > 0; k--)
    2151             :   {
    2152       80164 :     pari_sp av2 = avma;
    2153             :     long j;
    2154       80164 :     GEN t = gel(b,k);
    2155       80164 :     if (typ(t) != t_INT) pari_err_TYPE("hnf_invimage",t);
    2156       80164 :     for (j=1; j<=n; j++) t = subii(t, mulii(gcoeff(A,k,j),gel(u,j)));
    2157       80164 :     if (signe(t)) { avma = av;return NULL; }
    2158       80164 :     avma = av2;
    2159             :   }
    2160       80164 :   return u;
    2161             : }
    2162             : 
    2163             : /* A upper HNF, B integral matrix or column. Return A^(-1) B if integral,
    2164             :  * NULL otherwise */
    2165             : GEN
    2166      106876 : hnf_solve(GEN A, GEN B)
    2167             : {
    2168             :   pari_sp av;
    2169             :   long i, l;
    2170             :   GEN C;
    2171             : 
    2172      106876 :   if (typ(B) == t_COL) return hnf_invimage(A, B);
    2173       40271 :   av = avma; C = cgetg_copy(B, &l);
    2174       63399 :   for (i = 1; i < l; i++) {
    2175       42189 :     GEN c = hnf_invimage(A, gel(B,i));
    2176       42189 :     if (!c) { avma = av; return NULL; }
    2177       23128 :     gel(C,i) = c;
    2178             :   }
    2179       21210 :   return C;
    2180             : }
    2181             : 
    2182             : /***************************************************************/
    2183             : /**                                                           **/
    2184             : /**               SMITH NORMAL FORM REDUCTION                 **/
    2185             : /**                                                           **/
    2186             : /***************************************************************/
    2187             : 
    2188             : static GEN
    2189           0 : trivsmith(long all)
    2190             : {
    2191             :   GEN z;
    2192           0 :   if (!all) return cgetg(1,t_VEC);
    2193           0 :   z=cgetg(4,t_VEC);
    2194           0 :   gel(z,1) = cgetg(1,t_MAT);
    2195           0 :   gel(z,2) = cgetg(1,t_MAT);
    2196           0 :   gel(z,3) = cgetg(1,t_MAT); return z;
    2197             : }
    2198             : 
    2199             : static void
    2200      132769 : snf_pile(pari_sp av, GEN *x, GEN *U, GEN *V)
    2201             : {
    2202             :   GEN *gptr[3];
    2203      132769 :   int c = 1; gptr[0]=x;
    2204      132769 :   if (*U) gptr[c++] = U;
    2205      132769 :   if (*V) gptr[c++] = V;
    2206      132769 :   gerepilemany(av,gptr,c);
    2207      132769 : }
    2208             : 
    2209             : static GEN
    2210       77247 : bezout_step(GEN *pa, GEN *pb, GEN *pu, GEN *pv)
    2211             : {
    2212       77247 :   GEN a = *pa, b = *pb, d;
    2213       77247 :   if (absequalii(a,b))
    2214             :   {
    2215       29281 :     long sa = signe(a), sb = signe(b);
    2216       29281 :     *pv = gen_0;
    2217       29281 :     if (sb == sa) {
    2218       29085 :       *pa = *pb = gen_1;
    2219       29085 :       if (sa > 0) {*pu=gen_1; return a;} else {*pu=gen_m1; return absi(a);}
    2220             :     }
    2221         196 :     if (sa > 0) { *pa = *pu = gen_1; *pb = gen_m1; return a; }
    2222           7 :     *pa = *pu = gen_m1; *pb = gen_1; return b;
    2223             :   }
    2224       47966 :   d = bezout(a,b, pu,pv);
    2225       47966 :   *pa = diviiexact(a, d);
    2226       47966 :   *pb = diviiexact(b, d); return d;
    2227             : }
    2228             : 
    2229             : static int
    2230      359254 : negcmpii(void *E, GEN x, GEN y) { (void)E; return -cmpii(x,y); }
    2231             : 
    2232             : /* does b = x[i,i] divide all entries in x[1..i-1,1..i-1] ? If so, return 0;
    2233             :  * else return the index of a problematic row */
    2234             : static long
    2235      156814 : ZM_snf_no_divide(GEN x, long i)
    2236             : {
    2237      156814 :   GEN b = gcoeff(x,i,i);
    2238             :   long j, k;
    2239             : 
    2240      156814 :   if (!signe(b))
    2241             :   { /* impossible in the current implementation : x square of maximal rank */
    2242           0 :     for (k = 1; k < i; k++)
    2243           0 :       for (j = 1; j < i; j++)
    2244           0 :         if (signe(gcoeff(x,k,j))) return k;
    2245           0 :     return 0;
    2246             :   }
    2247      156814 :   if (is_pm1(b)) return 0;
    2248      285652 :   for (k=1; k<i; k++)
    2249             :   {
    2250     1051612 :     for (j=1; j<i; j++)
    2251      844305 :       if (signe(remii(gcoeff(x,k,j),b))) return k;
    2252             :   }
    2253       75615 :   return 0;
    2254             : }
    2255             : 
    2256             : /* Return the SNF D of matrix X. If ptU/ptV non-NULL set them to U/V
    2257             :  * to that D = UXV */
    2258             : GEN
    2259       75881 : ZM_snfall_i(GEN x, GEN *ptU, GEN *ptV, int return_vec)
    2260             : {
    2261       75881 :   pari_sp av0 = avma, av;
    2262             :   long i, j, k, m0, m, n0, n;
    2263       75881 :   GEN p1, u, v, U, V, V0, mdet, ys, perm = NULL;
    2264             : 
    2265       75881 :   n0 = n = lg(x)-1;
    2266       75881 :   if (!n) {
    2267        3592 :     if (ptU) *ptU = cgetg(1,t_MAT);
    2268        3592 :     if (ptV) *ptV = cgetg(1,t_MAT);
    2269        3592 :     return cgetg(1, return_vec? t_VEC: t_MAT);
    2270             :   }
    2271       72289 :   av = avma;
    2272       72289 :   m0 = m = nbrows(x);
    2273             : 
    2274       72289 :   U = ptU? gen_1: NULL; /* TRANSPOSE of row transform matrix [act on columns]*/
    2275       72289 :   V = ptV? gen_1: NULL;
    2276       72289 :   V0 = NULL;
    2277       72289 :   x = RgM_shallowcopy(x);
    2278       72289 :   if (m == n && ZM_ishnf(x))
    2279             :   {
    2280       64498 :     mdet = ZM_det_triangular(x);
    2281       64498 :     if (V) *ptV = matid(n);
    2282             :   }
    2283             :   else
    2284             :   {
    2285        7791 :     mdet = ZM_detmult(x);
    2286        7791 :     if (signe(mdet))
    2287             :     {
    2288        7770 :       if (!V)
    2289          77 :         p1 = ZM_hnfmod(x,mdet);
    2290             :       else
    2291             :       {
    2292        7693 :         if (m == n)
    2293             :         {
    2294        7672 :           p1 = ZM_hnfmod(x,mdet);
    2295        7672 :           *ptV = RgM_solve(x,p1);
    2296             :         }
    2297             :         else
    2298          21 :           p1 = ZM_hnfperm(x, ptV, ptU? &perm: NULL);
    2299             :       }
    2300        7770 :       mdet = ZM_det_triangular(p1);
    2301             :     }
    2302             :     else
    2303          21 :       p1 = ZM_hnfperm(x, ptV, ptU? &perm: NULL);
    2304        7791 :     x = p1;
    2305             :   }
    2306       72289 :   n = lg(x)-1;
    2307       72289 :   if (V)
    2308             :   {
    2309       68446 :     V = *ptV;
    2310       68446 :     if (n != n0)
    2311             :     {
    2312          21 :       V0 = vecslice(V, 1, n0 - n); /* kernel */
    2313          21 :       V  = vecslice(V, n0-n+1, n0);
    2314          21 :       av = avma;
    2315             :     }
    2316             :   }
    2317             :   /* independent columns */
    2318       72289 :   if (!signe(mdet))
    2319             :   {
    2320          21 :     if (n)
    2321             :     {
    2322          21 :       x = ZM_snfall_i(shallowtrans(x), ptV, ptU, return_vec); /* swap ptV,ptU */
    2323          21 :       if (typ(x) == t_MAT && n != m) x = shallowtrans(x);
    2324          21 :       if (V) V = ZM_mul(V, shallowtrans(*ptV));
    2325          21 :       if (U) U = *ptU; /* TRANSPOSE */
    2326             :     }
    2327             :     else /* 0 matrix */
    2328             :     {
    2329           0 :       x = cgetg(1,t_MAT);
    2330           0 :       if (V) V = cgetg(1, t_MAT);
    2331           0 :       if (U) U = matid(m);
    2332             :     }
    2333          21 :     goto THEEND;
    2334             :   }
    2335       72268 :   if (U) U = matid(n);
    2336             : 
    2337             :   /* square, maximal rank n */
    2338       72268 :   p1 = gen_indexsort(RgM_diagonal_shallow(x), NULL, &negcmpii);
    2339       72268 :   ys = cgetg(n+1,t_MAT);
    2340       72268 :   for (j=1; j<=n; j++) gel(ys,j) = vecpermute(gel(x, p1[j]), p1);
    2341       72268 :   x = ys;
    2342       72268 :   if (U) U = vecpermute(U, p1);
    2343       72268 :   if (V) V = vecpermute(V, p1);
    2344             : 
    2345       72268 :   p1 = ZM_hnfmod(x, mdet);
    2346       72268 :   if (V) V = ZM_mul(V, RgM_solve(x,p1));
    2347       72268 :   x = p1;
    2348             : 
    2349       72268 :   if (DEBUGLEVEL>7) err_printf("starting SNF loop");
    2350      452704 :   for (i=n; i>1; i--)
    2351             :   {
    2352      154084 :     if (DEBUGLEVEL>7) err_printf("\ni = %ld: ",i);
    2353             :     for(;;)
    2354             :     {
    2355      211269 :       int c = 0;
    2356             :       GEN a, b;
    2357      927921 :       for (j=i-1; j>=1; j--)
    2358             :       {
    2359      716652 :         b = gcoeff(x,i,j); if (!signe(b)) continue;
    2360        8822 :         a = gcoeff(x,i,i);
    2361        8822 :         ZC_elem(b, a, x,V, j,i);
    2362        8822 :         if (gc_needed(av,1))
    2363             :         {
    2364           0 :           if (DEBUGMEM>1) pari_warn(warnmem,"[1]: ZM_snfall i = %ld", i);
    2365           0 :           snf_pile(av, &x,&U,&V);
    2366             :         }
    2367             :       }
    2368      211269 :       if (DEBUGLEVEL>7) err_printf("; ");
    2369      927921 :       for (j=i-1; j>=1; j--)
    2370             :       {
    2371             :         GEN d;
    2372      716652 :         b = gcoeff(x,j,i); if (!signe(b)) continue;
    2373       77247 :         a = gcoeff(x,i,i);
    2374       77247 :         d = bezout_step(&a, &b, &u, &v);
    2375      451257 :         for (k = 1; k < i; k++)
    2376             :         {
    2377      374010 :           GEN t = addii(mulii(u,gcoeff(x,i,k)),mulii(v,gcoeff(x,j,k)));
    2378      748020 :           gcoeff(x,j,k) = subii(mulii(a,gcoeff(x,j,k)),
    2379      374010 :                                 mulii(b,gcoeff(x,i,k)));
    2380      374010 :           gcoeff(x,i,k) = t;
    2381             :         }
    2382       77247 :         gcoeff(x,j,i) = gen_0;
    2383       77247 :         gcoeff(x,i,i) = d;
    2384       77247 :         if (U) update(u,v,a,b,(GEN*)(U+i),(GEN*)(U+j));
    2385       77247 :         if (gc_needed(av,1))
    2386             :         {
    2387           0 :           if (DEBUGMEM>1) pari_warn(warnmem,"[2]: ZM_snfall, i = %ld", i);
    2388           0 :           snf_pile(av, &x,&U,&V);
    2389             :         }
    2390       77247 :         c = 1;
    2391             :       }
    2392      211269 :       if (!c)
    2393             :       {
    2394      156814 :         k = ZM_snf_no_divide(x, i);
    2395      156814 :         if (!k) break;
    2396             : 
    2397             :         /* x[k,j] != 0 mod b */
    2398        9219 :         for (j=1; j<=i; j++)
    2399        6489 :           gcoeff(x,i,j) = addii(gcoeff(x,i,j),gcoeff(x,k,j));
    2400        2730 :         if (U) gel(U,i) = gadd(gel(U,i),gel(U,k));
    2401             :       }
    2402       57185 :       if (gc_needed(av,1))
    2403             :       {
    2404           0 :         if (DEBUGMEM>1) pari_warn(warnmem,"[3]: ZM_snfall");
    2405           0 :         snf_pile(av, &x,&U,&V);
    2406             :       }
    2407       57185 :     }
    2408             :   }
    2409       72268 :   if (DEBUGLEVEL>7) err_printf("\n");
    2410      298620 :   for (k=1; k<=n; k++)
    2411      226352 :     if (signe(gcoeff(x,k,k)) < 0)
    2412             :     {
    2413         112 :       if (V) ZV_togglesign(gel(V,k));
    2414         112 :       togglesign(gcoeff(x,k,k));
    2415             :     }
    2416             : THEEND:
    2417       72289 :   if (return_vec)
    2418             :   {
    2419       67275 :     long l = lg(x)-1;
    2420       67275 :     if (typ(x) == t_MAT) x = RgM_diagonal_shallow(x);
    2421       67275 :     if (m0 > l) x = shallowconcat(zerovec(m0-l), x);
    2422             :   }
    2423             : 
    2424       72289 :   if (V0)
    2425             :   {
    2426          21 :     if (!return_vec) x = shallowconcat(zeromat(m,n0-n), x);
    2427          21 :     if (V) V = shallowconcat(V0, V);
    2428             :   }
    2429       72289 :   if (U)
    2430             :   {
    2431       55918 :     U = shallowtrans(U);
    2432       55918 :     if (perm) U = vecpermute(U, perm_inv(perm));
    2433             :   }
    2434       72289 :   snf_pile(av0, &x,&U,&V);
    2435       72289 :   if (ptU) *ptU = U;
    2436       72289 :   if (ptV) *ptV = V;
    2437       72289 :   return x;
    2438             : }
    2439             : GEN
    2440        8143 : ZM_snfall(GEN x, GEN *U, GEN *V) { return ZM_snfall_i(x, U, V, 0); }
    2441             : GEN
    2442         217 : ZM_snf(GEN x) { return ZM_snfall_i(x, NULL,NULL, 1); }
    2443             : 
    2444             : GEN
    2445         322 : smith(GEN x) {
    2446         322 :   if (typ(x)!=t_MAT) pari_err_TYPE("smith",x);
    2447         322 :   RgM_check_ZM(x, "smith");
    2448         322 :   return ZM_snfall_i(x, NULL,NULL, 1);
    2449             : }
    2450             : GEN
    2451          21 : smithall(GEN x)
    2452             : {
    2453          21 :   GEN z = cgetg(4, t_VEC);
    2454          21 :   if (typ(x)!=t_MAT) pari_err_TYPE("smithall",x);
    2455          21 :   RgM_check_ZM(x, "smithall");
    2456          21 :   gel(z,3) = ZM_snfall_i(x, (GEN*)(z+1),(GEN*)(z+2), 0);
    2457          21 :   return z;
    2458             : }
    2459             : 
    2460             : void
    2461      124529 : ZM_snfclean(GEN d, GEN u, GEN v)
    2462             : {
    2463      124529 :   long i, c, l = lg(d);
    2464             : 
    2465      124529 :   if (typ(d) == t_VEC)
    2466      124529 :     for (c=1; c<l; c++) { GEN t = gel(d,c); if (is_pm1(t)) break; }
    2467             :   else
    2468             :   {
    2469           0 :     for (c=1; c<l; c++) { GEN t = gcoeff(d,c,c); if (is_pm1(t)) break; }
    2470           0 :     if (c < l) for (i = 1; i < c; i++) setlg(gel(d,i), c);
    2471             :   }
    2472      124529 :   setlg(d, c);
    2473      124529 :   if (u) for (i=1; i<l; i++) setlg(gel(u,i), c);
    2474      124529 :   if (v) setlg(v, c);
    2475      124529 : }
    2476             : 
    2477             : /* Assume z was computed by [g]smithall(). Remove the 1s on the diagonal */
    2478             : GEN
    2479         539 : smithclean(GEN z)
    2480             : {
    2481             :   long i, j, h, l, c, d;
    2482             :   GEN U, V, y, D, t;
    2483             : 
    2484         539 :   if (typ(z) != t_VEC) pari_err_TYPE("smithclean",z);
    2485         539 :   l = lg(z); if (l == 1) return cgetg(1,t_VEC);
    2486         532 :   U = gel(z,1);
    2487         532 :   if (l != 4 || typ(U) != t_MAT)
    2488             :   { /* assume z = vector of elementary divisors */
    2489        1470 :     for (c=1; c<l; c++)
    2490        1309 :       if (gequal1(gel(z,c))) break;
    2491         511 :     return gcopy_lg(z, c);
    2492             :   }
    2493          21 :   V = gel(z,2);
    2494          21 :   D = gel(z,3);
    2495          21 :   l = lg(D);
    2496          21 :   if (l == 1) return gcopy(z);
    2497          21 :   h = lgcols(D);
    2498          21 :   if (h > l)
    2499             :   { /* D = vconcat(zero matrix, diagonal matrix) */
    2500          21 :     for (c=1+h-l, d=1; c<h; c++,d++)
    2501          21 :       if (gequal1(gcoeff(D,c,d))) break;
    2502             :   }
    2503           7 :   else if (h < l)
    2504             :   { /* D = concat(zero matrix, diagonal matrix) */
    2505           7 :     for (c=1, d=1+l-h; d<l; c++,d++)
    2506           7 :       if (gequal1(gcoeff(D,c,d))) break;
    2507             :   }
    2508             :   else
    2509             :   { /* D diagonal */
    2510           0 :     for (c=1; c<l; c++)
    2511           0 :       if (gequal1(gcoeff(D,c,c))) break;
    2512           0 :     d = c;
    2513             :   }
    2514             :   /* U was (h-1)x(h-1), V was (l-1)x(l-1), D was (h-1)x(l-1) */
    2515          21 :   y = cgetg(4,t_VEC);
    2516             :   /* truncate U to (c-1) x (h-1) */
    2517          21 :   gel(y,1) = t = cgetg(h,t_MAT);
    2518          21 :   for (j=1; j<h; j++) gel(t,j) = gcopy_lg(gel(U,j), c);
    2519             :   /* truncate V to (l-1) x (d-1) */
    2520          21 :   gel(y,2) = gcopy_lg(V, d);
    2521          21 :   gel(y,3) = t = zeromatcopy(c-1, d-1);
    2522             :   /* truncate D to a (c-1) x (d-1) matrix */
    2523          21 :   if (d > 1)
    2524             :   {
    2525          14 :     if (h > l)
    2526             :     {
    2527          14 :       for (i=1+h-l, j=1; i<c; i++,j++)
    2528           7 :         gcoeff(t,i,j) = gcopy(gcoeff(D,i,j));
    2529             :     }
    2530           7 :     else if (h < l)
    2531             :     {
    2532           7 :       for (i=1, j=1+l-h; j<d; i++,j++)
    2533           0 :         gcoeff(t,i,j) = gcopy(gcoeff(D,i,j));
    2534             :     }
    2535             :     else
    2536             :     {
    2537           0 :       for (j=1; j<d; j++)
    2538           0 :         gcoeff(t,j,j) = gcopy(gcoeff(D,j,j));
    2539             :     }
    2540             :   }
    2541          21 :   return y;
    2542             : }
    2543             : 
    2544             : /* does b = x[i,i] divide all entries in x[1..i-1,1..i-1] ? If so, return 0;
    2545             :  * else return the index of a problematic row */
    2546             : static long
    2547         112 : gsnf_no_divide(GEN x, long i, long vx)
    2548             : {
    2549         112 :   GEN b = gcoeff(x,i,i);
    2550             :   long j, k;
    2551             : 
    2552         112 :   if (gequal0(b))
    2553             :   {
    2554          14 :     for (k = 1; k < i; k++)
    2555          14 :       for (j = 1; j < i; j++)
    2556          14 :         if (!gequal0(gcoeff(x,k,j))) return k;
    2557           0 :     return 0;
    2558             :   }
    2559             : 
    2560          98 :   if (!is_RgX(b,vx) || degpol(b)<=0) return 0;
    2561         133 :   for (k = 1; k < i; k++)
    2562         231 :     for (j = 1; j < i; j++)
    2563             :     {
    2564         154 :       GEN z = gcoeff(x,k,j), r;
    2565         154 :       if (!is_RgX(z,vx)) z = scalarpol(z, vx);
    2566         154 :       r = RgX_rem(z, b);
    2567         154 :       if (signe(r) && (! isinexactreal(r) ||
    2568           0 :              gexpo(r) > 16 + gexpo(b) - prec2nbits(gprecision(r)))
    2569           7 :          ) return k;
    2570             :     }
    2571          49 :   return 0;
    2572             : }
    2573             : 
    2574             : /* Hermite Normal Form, with base change matrix if ptB != NULL.
    2575             :  * If 'remove' = 1, remove 0 columns (do NOT update *ptB accordingly)
    2576             :  * If 'remove' = 2, remove 0 columns and update *ptB accordingly */
    2577             : GEN
    2578          35 : RgM_hnfall(GEN A, GEN *pB, long remove)
    2579             : {
    2580             :   pari_sp av;
    2581             :   long li, j, k, m, n, def, ldef;
    2582             :   GEN B;
    2583          35 :   long vx = gvar(A);
    2584             : 
    2585          35 :   n = lg(A)-1;
    2586          35 :   if (vx==NO_VARIABLE || !n)
    2587             :   {
    2588           0 :     RgM_check_ZM(A, "mathnf0");
    2589           0 :     return ZM_hnfall(A, pB, remove);
    2590             :   }
    2591          35 :   m = nbrows(A);
    2592          35 :   av = avma;
    2593          35 :   A = RgM_shallowcopy(A);
    2594          35 :   B = pB? matid(n): NULL;
    2595          35 :   def = n; ldef = (m>n)? m-n: 0;
    2596          77 :   for (li=m; li>ldef; li--)
    2597             :   {
    2598             :     GEN d, T;
    2599          63 :     for (j=def-1; j; j--)
    2600             :     {
    2601          21 :       GEN a = gcoeff(A,li,j);
    2602          21 :       if (gequal0(a)) continue;
    2603             : 
    2604           7 :       k = (j==1)? def: j-1;
    2605           7 :       RgC_elem(a,gcoeff(A,li,k), A,B, j,k, li, vx);
    2606             :     }
    2607          42 :     T = normalize_as_RgX(gcoeff(A,li,def), vx, &d);
    2608          42 :     if (gequal0(T))
    2609           0 :     { if (ldef) ldef--; }
    2610             :     else
    2611             :     {
    2612          42 :       gcoeff(A,li,def) = T;
    2613          42 :       if (B && !gequal1(d)) gel(B, def) = RgC_Rg_div(gel(B, def), d);
    2614          42 :       RgM_reduce(A, B, li, def, vx);
    2615          42 :       def--;
    2616             :     }
    2617          42 :     if (gc_needed(av,1))
    2618             :     {
    2619           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"ghnfall");
    2620           0 :       gerepileall(av, B? 2: 1, &A, &B);
    2621             :     }
    2622             :   }
    2623             :   /* rank A = n - def */
    2624          35 :   if (remove) remove_0cols(def, &A, &B, remove);
    2625          35 :   gerepileall(av, B? 2: 1, &A, &B);
    2626          35 :   if (B) *pB = B;
    2627          35 :   return A;
    2628             : }
    2629             : 
    2630             : static GEN
    2631          42 : gsmithall_i(GEN x,long all)
    2632             : {
    2633             :   pari_sp av;
    2634             :   long i, j, k, n;
    2635             :   GEN z, u, v, U, V;
    2636          42 :   long vx = gvar(x);
    2637          42 :   if (typ(x)!=t_MAT) pari_err_TYPE("gsmithall",x);
    2638          42 :   if (vx==NO_VARIABLE) return all? smithall(x): smith(x);
    2639          35 :   n = lg(x)-1;
    2640          35 :   if (!n) return trivsmith(all);
    2641          35 :   if (lgcols(x) != n+1) pari_err_DIM("gsmithall");
    2642          35 :   av = avma;
    2643          35 :   x = RgM_shallowcopy(x);
    2644          35 :   if (all) { U = matid(n); V = matid(n); }
    2645         252 :   for (i=n; i>=2; i--)
    2646             :   {
    2647             :     for(;;)
    2648             :     {
    2649             :       GEN a, b, d;
    2650         189 :       int c = 0;
    2651         546 :       for (j=i-1; j>=1; j--)
    2652             :       {
    2653         357 :         b = gcoeff(x,i,j); if (gequal0(b)) continue;
    2654         140 :         a = gcoeff(x,i,i);
    2655         140 :         d = gbezout_step(&b, &a, &v, &u, vx);
    2656         469 :         for (k = 1; k < i; k++)
    2657             :         {
    2658         329 :           GEN t = gadd(gmul(u,gcoeff(x,k,i)),gmul(v,gcoeff(x,k,j)));
    2659         329 :           gcoeff(x,k,j) = gsub(gmul(a,gcoeff(x,k,j)),gmul(b,gcoeff(x,k,i)));
    2660         329 :           gcoeff(x,k,i) = t;
    2661             :         }
    2662         140 :         gcoeff(x,i,j) = gen_0;
    2663         140 :         gcoeff(x,i,i) = d;
    2664         140 :         if (all) update(u,v,a,b,(GEN*)(V+i),(GEN*)(V+j));
    2665             :       }
    2666         546 :       for (j=i-1; j>=1; j--)
    2667             :       {
    2668         357 :         b = gcoeff(x,j,i); if (gequal0(b)) continue;
    2669         119 :         a = gcoeff(x,i,i);
    2670         119 :         d = gbezout_step(&b, &a, &v, &u, vx);
    2671         406 :         for (k = 1; k < i; k++)
    2672             :         {
    2673         287 :           GEN t = gadd(gmul(u,gcoeff(x,i,k)),gmul(v,gcoeff(x,j,k)));
    2674         287 :           gcoeff(x,j,k) = gsub(gmul(a,gcoeff(x,j,k)),gmul(b,gcoeff(x,i,k)));
    2675         287 :           gcoeff(x,i,k) = t;
    2676             :         }
    2677         119 :         gcoeff(x,j,i) = gen_0;
    2678         119 :         gcoeff(x,i,i) = d;
    2679         119 :         if (all) update(u,v,a,b,(GEN*)(U+i),(GEN*)(U+j));
    2680         119 :         c = 1;
    2681             :       }
    2682         189 :       if (!c)
    2683             :       {
    2684         112 :         k = gsnf_no_divide(x, i, vx);
    2685         112 :         if (!k) break;
    2686             : 
    2687          70 :         for (j=1; j<=i; j++)
    2688          49 :           gcoeff(x,i,j) = gadd(gcoeff(x,i,j),gcoeff(x,k,j));
    2689          21 :         if (all) gel(U,i) = gadd(gel(U,i),gel(U,k));
    2690             :       }
    2691          98 :       if (gc_needed(av,1))
    2692             :       {
    2693           0 :         if (DEBUGMEM>1) pari_warn(warnmem,"gsmithall");
    2694           0 :         gerepileall(av, all? 3: 1, &x, &U, &V);
    2695             :       }
    2696          98 :     }
    2697             :   }
    2698         161 :   for (k=1; k<=n; k++)
    2699             :   {
    2700         126 :     GEN d, T = normalize_as_RgX(gcoeff(x,k,k), vx, &d);
    2701         126 :     if (gequal0(T)) continue;
    2702         112 :     if (all && !gequal1(d)) gel(V,k) = RgC_Rg_div(gel(V,k), d);
    2703         112 :     gcoeff(x,k,k) = T;
    2704             :   }
    2705          35 :   z = all? mkvec3(shallowtrans(U), V, x): RgM_diagonal_shallow(x);
    2706          35 :   return gerepilecopy(av, z);
    2707             : }
    2708             : 
    2709             : GEN
    2710         378 : matsnf0(GEN x,long flag)
    2711             : {
    2712         378 :   pari_sp av = avma;
    2713         378 :   if (flag > 7) pari_err_FLAG("matsnf");
    2714         378 :   if (typ(x) == t_VEC && flag & 4) return smithclean(x);
    2715         378 :   if (flag & 2) x = flag&1 ? gsmithall(x): gsmith(x);
    2716         336 :   else          x = flag&1 ?  smithall(x):  smith(x);
    2717         378 :   if (flag & 4) x = gerepileupto(av, smithclean(x));
    2718         378 :   return x;
    2719             : }
    2720             : 
    2721             : GEN
    2722          42 : gsmith(GEN x) { return gsmithall_i(x,0); }
    2723             : 
    2724             : GEN
    2725           0 : gsmithall(GEN x) { return gsmithall_i(x,1); }
    2726             : 
    2727             : /* H is a relation matrix, either in HNF or a t_VEC (diagonal HNF) */
    2728             : static GEN
    2729      124529 : snf_group(GEN H, GEN D, GEN *newU, GEN *newUi)
    2730             : {
    2731             :   long i, j, l;
    2732             : 
    2733      124529 :   ZM_snfclean(D, newU? *newU: NULL, newUi? *newUi: NULL);
    2734      124529 :   l = lg(D);
    2735      124529 :   if (newU) {
    2736      108276 :     GEN U = *newU;
    2737      346185 :     for (i = 1; i < l; i++)
    2738             :     {
    2739      237909 :       GEN d = gel(D,i), d2 = shifti(d, 1);
    2740     1282141 :       for (j = 1; j < lg(U); j++)
    2741     1044232 :         gcoeff(U,i,j) = centermodii(gcoeff(U,i,j), d, d2);
    2742             :     }
    2743      108276 :     *newU = U;
    2744             :   }
    2745      124529 :   if (newUi && l > 1)
    2746             :   { /* UHV=D -> U^-1 = (HV)D^-1 -> U^-1 = H(VD^-1 mod 1) mod H */
    2747             :     /* Ui = ZM_inv(U, NULL); setlg(Ui, l); */
    2748      120454 :     GEN V = *newUi, Ui;
    2749      120454 :     int Hvec = (typ(H) == t_VEC);
    2750      120454 :     for (i = 1; i < l; i++) gel(V,i) = FpC_red(gel(V,i), gel(D,i));
    2751      120454 :     if (!Hvec)
    2752             :     {
    2753       60170 :       if (ZM_isdiagonal(H)) { H = RgM_diagonal_shallow(H); Hvec = 1; }
    2754             :     }
    2755      120454 :     Ui = Hvec? ZM_diag_mul(H, V): ZM_mul(H, V);
    2756      120454 :     for (i = 1; i < l; i++) gel(Ui,i) = ZC_Z_divexact(gel(Ui,i), gel(D,i));
    2757      120454 :     if (Hvec)
    2758       84843 :     { for (i = 1; i < l; i++) gel(Ui,i) = vecmodii(gel(Ui,i), H); }
    2759             :     else
    2760       35611 :       Ui = ZM_hnfrem(Ui, H);
    2761      120454 :     *newUi = Ui;
    2762             :   }
    2763      124529 :   return D;
    2764             : }
    2765             : /* H relation matrix among row of generators g in HNF.  Let URV = D its SNF,
    2766             :  * newU R newV = newD its clean SNF (no 1 in Dnew). Return the diagonal of
    2767             :  * newD, newU and newUi such that  1/U = (newUi, ?).
    2768             :  * Rationale: let (G,0) = g Ui be the new generators then
    2769             :  * 0 = G U R --> G D = 0,  g = G newU  and  G = g newUi */
    2770             : GEN
    2771       64049 : ZM_snf_group(GEN H, GEN *newU, GEN *newUi)
    2772             : {
    2773       64049 :   GEN D = ZM_snfall_i(H, newU, newUi, 1);
    2774       64049 :   return snf_group(H, D, newU, newUi);
    2775             : }
    2776             : 
    2777             : /* D a ZV: SNF for matdiagonal(D). Faster because we only ensure elementary
    2778             :  * divisors condition: d[n] | ... | d[1] and need not convert D to matrix form*/
    2779             : GEN
    2780       60480 : ZV_snfall(GEN D, GEN *pU, GEN *pV)
    2781             : {
    2782       60480 :   pari_sp av = avma;
    2783       60480 :   long j, n = lg(D)-1;
    2784       60480 :   GEN U = pU? matid(n): NULL;
    2785       60480 :   GEN V = pV? matid(n): NULL;
    2786             :   GEN p;
    2787             : 
    2788       60480 :   D = leafcopy(D);
    2789      194005 :   for (j = n; j > 0; j--)
    2790             :   {
    2791      133525 :     GEN b = gel(D,j);
    2792      133525 :     if (signe(b) < 0)
    2793             :     {
    2794           0 :       gel(D,j) = absi(b);
    2795           0 :       if (V) ZV_togglesign(gel(V,j));
    2796             :     }
    2797             :   }
    2798             :   /* entries are non-negative integers */
    2799       60480 :   p = gen_indexsort(D, NULL, &negcmpii);
    2800       60480 :   D = vecpermute(D, p);
    2801       60480 :   if (U) U = vecpermute(U, p);
    2802       60480 :   if (V) V = vecpermute(V, p);
    2803             :   /* entries are sorted by decreasing value */
    2804      194005 :   for (j = n; j > 0; j--)
    2805             :   {
    2806      133525 :     GEN b = gel(D,j);
    2807             :     long i;
    2808      241038 :     for (i = j-1; i > 0; i--)
    2809             :     { /* invariant: a >= b. If au+bv = d is a Bezout relation, A=a/d and B=b/d
    2810             :        * we have [B,-A;u,v]*diag(a,b)*[1-u*A,1; -u*A,1]] = diag(Ab, d) */
    2811      108185 :       GEN a = gel(D,i), u,v, d = bezout(a,b, &u,&v), A, Wi, Wj;
    2812      108185 :       if (equalii(d,b)) continue;
    2813       13608 :       A = diviiexact(a,d);
    2814       13608 :       if (V)
    2815             :       {
    2816       13566 :         GEN t = mulii(u,A);
    2817       13566 :         Wi = ZC_lincomb(subui(1,t), negi(t), gel(V,i), gel(V,j));
    2818       13566 :         Wj = ZC_add(gel(V,i), gel(V,j));
    2819       13566 :         gel(V,i) = Wi;
    2820       13566 :         gel(V,j) = Wj;
    2821             :       }
    2822       13608 :       if (U)
    2823             :       {
    2824       13608 :         GEN B = diviiexact(b,d);
    2825       13608 :         Wi = ZC_lincomb(B, negi(A), gel(U,i), gel(U,j));
    2826       13608 :         Wj = ZC_lincomb(u, v, gel(U,i), gel(U,j));
    2827       13608 :         gel(U,i) = Wi;
    2828       13608 :         gel(U,j) = Wj;
    2829             :       }
    2830       13608 :       gel(D,i) = mulii(A,b); /* lcm(a,b) */
    2831       13608 :       gel(D,j) = d; /* gcd(a,b) */
    2832       13608 :       b = gel(D,j); if (equali1(b)) break;
    2833             :     }
    2834             :   }
    2835       60480 :   snf_pile(av, &D,&U,&V);
    2836       60480 :   if (U) *pU = shallowtrans(U);
    2837       60480 :   if (V) *pV = V;
    2838       60480 :   return D;
    2839             : }
    2840             : GEN
    2841       60480 : ZV_snf_group(GEN d, GEN *newU, GEN *newUi)
    2842             : {
    2843       60480 :   GEN D = ZV_snfall(d, newU, newUi);
    2844       60480 :   return snf_group(d, D, newU, newUi);
    2845             : }
    2846             : 
    2847             : /* D a vector of elementary divisors. Truncate (setlg) to leave out trivial
    2848             :  * entries (= 1) */
    2849             : void
    2850         224 : ZV_snf_trunc(GEN D)
    2851             : {
    2852         224 :   long i, l = lg(D);
    2853         490 :   for (i = 1; i < l; i++)
    2854         371 :     if (is_pm1(gel(D,i))) { setlg(D,i); break; }
    2855         224 : }

Generated by: LCOV version 1.11