Code coverage tests

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

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

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

LCOV - code coverage report
Current view: top level - basemath - hnf_snf.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.8.0 lcov report (development 16912-212c0f0) Lines: 1338 1485 90.1 %
Date: 2014-10-20 Functions: 71 78 91.0 %
Legend: Lines: hit not hit | Branches: + taken - not taken # not executed Branches: 1105 1385 79.8 %

           Branch data     Line data    Source code
       1                 :            : /* Copyright (C) 2000  The PARI group.
       2                 :            : 
       3                 :            : This file is part of the PARI/GP package.
       4                 :            : 
       5                 :            : PARI/GP is free software; you can redistribute it and/or modify it under the
       6                 :            : terms of the GNU General Public License as published by the Free Software
       7                 :            : Foundation. It is distributed in the hope that it will be useful, but WITHOUT
       8                 :            : ANY WARRANTY WHATSOEVER.
       9                 :            : 
      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                 :          7 : hnfallgen(GEN x)
      24                 :            : {
      25                 :          7 :   GEN z = cgetg(3, t_VEC);
      26                 :          7 :   gel(z,1) = RgM_hnfall(x, (GEN*)(z+2), 1);
      27                 :          7 :   return z;
      28                 :            : }
      29                 :            : GEN
      30                 :      11081 : mathnf0(GEN x, long flag)
      31                 :            : {
      32      [ +  +  - ]:      11081 :   switch(typ(x))
      33                 :            :   {
      34                 :            :     case t_VEC:
      35         [ +  + ]:      11018 :       if (RgV_is_ZV(x))
      36      [ +  +  + ]:      10997 :         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                 :      10976 :             return ZV_hnfgcdext(x);
      44                 :            :         }
      45                 :         28 :       x = gtomat(x); break;
      46                 :         63 :     case t_MAT: break;
      47                 :          0 :     default: pari_err_TYPE("mathnf0",x);
      48                 :            :   }
      49                 :            : 
      50   [ +  +  -  -  :         91 :   switch(flag)
                +  +  - ]
      51                 :            :   {
      52         [ +  + ]:         49 :     case 0: return RgM_is_ZM(x)? ZM_hnf(x): RgM_hnfall(x,NULL,1);
      53         [ +  + ]:         21 :     case 1: return RgM_is_ZM(x)? hnfall(x): hnfallgen(x);
      54                 :          0 :     case 2: return RgM_hnfall(x, NULL, 1);
      55                 :          0 :     case 3: return hnfallgen(x);
      56                 :          7 :     case 4: RgM_check_ZM(x, "mathnf0"); return hnflll(x);
      57                 :         14 :     case 5: RgM_check_ZM(x, "mathnf0"); return hnfperm(x);
      58                 :          0 :     default: pari_err_FLAG("mathnf");
      59                 :            :   }
      60                 :      11081 :   return NULL; /* not reached */
      61                 :            : }
      62                 :            : 
      63                 :            : /*******************************************************************/
      64                 :            : /*                                                                 */
      65                 :            : /*                SPECIAL HNF (FOR INTERNAL USE !!!)               */
      66                 :            : /*                                                                 */
      67                 :            : /*******************************************************************/
      68                 :            : static int
      69                 :    2703998 : count(GEN mat, long row, long len, long *firstnonzero)
      70                 :            : {
      71                 :    2703998 :   long j, n = 0;
      72                 :            : 
      73         [ +  + ]:  189423663 :   for (j=1; j<=len; j++)
      74                 :            :   {
      75                 :  187316112 :     long p = mael(mat,j,row);
      76         [ +  + ]:  187316112 :     if (p)
      77                 :            :     {
      78         [ +  + ]:    7359399 :       if (labs(p)!=1) return -1;
      79                 :    6762952 :       n++; *firstnonzero=j;
      80                 :            :     }
      81                 :            :   }
      82                 :    2703998 :   return n;
      83                 :            : }
      84                 :            : 
      85                 :            : static int
      86                 :      69021 : count2(GEN mat, long row, long len)
      87                 :            : {
      88                 :            :   long j;
      89         [ +  + ]:     799489 :   for (j=len; j; j--)
      90         [ +  + ]:     773263 :     if (labs(mael(mat,j,row)) == 1) return j;
      91                 :      69021 :   return 0;
      92                 :            : }
      93                 :            : 
      94                 :            : static GEN
      95                 :      15523 : hnffinal(GEN matgen,GEN perm,GEN* ptdep,GEN* ptB,GEN* ptC)
      96                 :            : {
      97                 :            :   GEN p1,p2,U,H,Hnew,Bnew,Cnew,diagH1;
      98                 :      15523 :   GEN B = *ptB, C = *ptC, dep = *ptdep, depnew;
      99                 :            :   pari_sp av;
     100                 :            :   long i,j,k,s,i1,j1,zc;
     101                 :      15523 :   long co = lg(C);
     102                 :      15523 :   long col = lg(matgen)-1;
     103                 :            :   long lnz, nlze, lig;
     104                 :            : 
     105         [ +  + ]:      15523 :   if (col == 0) return matgen;
     106                 :      15383 :   lnz = nbrows(matgen); /* was called lnz-1 - nr in hnfspec */
     107                 :      15383 :   nlze = nbrows(dep);
     108                 :      15383 :   lig = nlze + lnz;
     109                 :            :   /* H: lnz x lnz [disregarding initial 0 cols], U: col x col */
     110                 :      15383 :   H = ZM_hnflll(matgen, &U, 0);
     111                 :      15383 :   H += (lg(H)-1 - lnz); H[0] = evaltyp(t_MAT) | evallg(lnz+1);
     112                 :            :   /* Only keep the part above the H (above the 0s is 0 since the dep rows
     113                 :            :    * are dependent from the ones in matgen) */
     114                 :      15383 :   zc = col - lnz; /* # of 0 columns, correspond to units */
     115         [ +  + ]:      15383 :   if (nlze) { dep = ZM_mul(dep,U); dep += zc; }
     116                 :            : 
     117                 :      15383 :   diagH1 = new_chunk(lnz+1); /* diagH1[i] = 0 iff H[i,i] != 1 (set later) */
     118                 :            : 
     119                 :      15383 :   av = avma;
     120                 :      15383 :   Cnew = cgetg(co, typ(C));
     121                 :      15383 :   setlg(C, col+1); p1 = gmul(C,U);
     122         [ +  + ]:     166935 :   for (j=1; j<=col; j++) gel(Cnew,j) = gel(p1,j);
     123         [ +  + ]:     985023 :   for (   ; j<co ; j++)  gel(Cnew,j) = gel(C,j);
     124                 :            : 
     125                 :            :   /* Clean up B using new H */
     126         [ +  + ]:     103399 :   for (s=0,i=lnz; i; i--)
     127                 :            :   {
     128                 :      88016 :     GEN Di = gel(dep,i), Hi = gel(H,i);
     129                 :      88016 :     GEN h = gel(Hi,i); /* H[i,i] */
     130         [ +  + ]:      88016 :     if ( (diagH1[i] = is_pm1(h)) ) { h = NULL; s++; }
     131         [ +  + ]:    6576607 :     for (j=col+1; j<co; j++)
     132                 :            :     {
     133                 :    6488591 :       GEN z = gel(B,j-col);
     134                 :    6488591 :       p1 = gel(z,i+nlze);
     135         [ +  + ]:    6488591 :       if (h) p1 = truedivii(p1,h);
     136         [ +  + ]:    6488591 :       if (!signe(p1)) continue;
     137         [ +  + ]:   50565175 :       for (k=1; k<=nlze; k++) gel(z,k) = subii(gel(z,k), mulii(p1, gel(Di,k)));
     138         [ +  + ]:   43146618 :       for (   ; k<=lig;  k++) gel(z,k) = subii(gel(z,k), mulii(p1, gel(Hi,k-nlze)));
     139                 :    3962234 :       gel(Cnew,j) = gsub(gel(Cnew,j), gmul(p1, gel(Cnew,i+zc)));
     140                 :            :     }
     141         [ +  + ]:      88016 :     if (gc_needed(av,3))
     142                 :            :     {
     143         [ -  + ]:         37 :       if(DEBUGMEM>1) pari_warn(warnmem,"hnffinal, i = %ld",i);
     144                 :         37 :       gerepileall(av, 2, &Cnew, &B);
     145                 :            :     }
     146                 :            :   }
     147                 :      15383 :   p1 = cgetg(lnz+1,t_VEC); p2 = perm + nlze;
     148         [ +  + ]:     103399 :   for (i1=0, j1=lnz-s, i=1; i<=lnz; i++) /* push the 1 rows down */
     149         [ +  + ]:      88016 :     if (diagH1[i])
     150                 :      46241 :       gel(p1,++j1) = gel(p2,i);
     151                 :            :     else
     152                 :      41775 :       gel(p2,++i1) = gel(p2,i);
     153         [ +  + ]:      61624 :   for (i=i1+1; i<=lnz; i++) gel(p2,i) = gel(p1,i);
     154                 :            : 
     155                 :            :   /* s = # extra redundant generators taken from H
     156                 :            :    *          zc  col-s  co   zc = col - lnz
     157                 :            :    *       [ 0 |dep |     ]    i = nlze + lnz - s = lig - s
     158                 :            :    *  nlze [--------|  B' ]
     159                 :            :    *       [ 0 | H' |     ]    H' = H minus the s rows with a 1 on diagonal
     160                 :            :    *     i [--------|-----] lig-s           (= "1-rows")
     161                 :            :    *       [   0    | Id  ]
     162                 :            :    *       [        |     ] li */
     163                 :      15383 :   lig -= s; col -= s; lnz -= s;
     164                 :      15383 :   Hnew = cgetg(lnz+1,t_MAT);
     165                 :      15383 :   depnew = cgetg(lnz+1,t_MAT); /* only used if nlze > 0 */
     166                 :      15383 :   Bnew = cgetg(co-col,t_MAT);
     167                 :      15383 :   C = shallowcopy(Cnew);
     168         [ +  + ]:     103399 :   for (j=1,i1=j1=0; j<=lnz+s; j++)
     169                 :            :   {
     170                 :      88016 :     GEN z = gel(H,j);
     171         [ +  + ]:      88016 :     if (diagH1[j])
     172                 :            :     { /* hit exactly s times */
     173                 :      46241 :       i1++; C[i1+col] = Cnew[j+zc];
     174                 :      46241 :       p1 = cgetg(lig+1,t_COL); gel(Bnew,i1) = p1;
     175         [ +  + ]:     604708 :       for (i=1; i<=nlze; i++) gel(p1,i) = gcoeff(dep,i,j);
     176                 :      46241 :       p1 += nlze;
     177                 :            :     }
     178                 :            :     else
     179                 :            :     {
     180                 :      41775 :       j1++; C[j1+zc] = Cnew[j+zc];
     181                 :      41775 :       p1 = cgetg(lnz+1,t_COL); gel(Hnew,j1) = p1;
     182                 :      41775 :       depnew[j1] = dep[j];
     183                 :            :     }
     184         [ +  + ]:     541238 :     for (i=k=1; k<=lnz; i++)
     185         [ +  + ]:     453222 :       if (!diagH1[i]) p1[k++] = z[i];
     186                 :            :   }
     187         [ +  + ]:     985023 :   for (j=s+1; j<co-col; j++)
     188                 :            :   {
     189                 :     969640 :     GEN z = gel(B,j-s);
     190                 :     969640 :     p1 = cgetg(lig+1,t_COL); gel(Bnew,j) = p1;
     191         [ +  + ]:   10415237 :     for (i=1; i<=nlze; i++) gel(p1,i) = gel(z,i);
     192                 :     969640 :     z += nlze; p1 += nlze;
     193         [ +  + ]:    4764577 :     for (i=k=1; k<=lnz; i++)
     194         [ +  + ]:    3794937 :       if (!diagH1[i]) gel(p1,k++) = gel(z,i);
     195                 :            :   }
     196                 :      15383 :   *ptdep = depnew;
     197                 :      15383 :   *ptC = C;
     198                 :      15523 :   *ptB = Bnew; return Hnew;
     199                 :            : }
     200                 :            : 
     201                 :            : /* for debugging */
     202                 :            : static void
     203                 :          0 : p_mat(GEN mat, GEN perm, long k)
     204                 :            : {
     205                 :          0 :   pari_sp av = avma;
     206                 :          0 :   perm = vecslice(perm, k+1, lg(perm)-1);
     207                 :          0 :   err_printf("Permutation: %Ps\n",perm);
     208         [ #  # ]:          0 :   if (DEBUGLEVEL > 6)
     209                 :          0 :     err_printf("matgen = %Ps\n", zm_to_ZM( rowpermute(mat, perm) ));
     210                 :          0 :   avma = av;
     211                 :          0 : }
     212                 :            : 
     213                 :            : static GEN
     214                 :     405027 : col_dup(long l, GEN col)
     215                 :            : {
     216                 :     405027 :   GEN c = new_chunk(l);
     217                 :     405027 :   memcpy(c,col,l * sizeof(long)); return c;
     218                 :            : }
     219                 :            : 
     220                 :            : /* HNF reduce a relation matrix (column operations + row permutation)
     221                 :            : ** Input:
     222                 :            : **   mat = (li-1) x (co-1) matrix of long
     223                 :            : **   C   = r x (co-1) matrix of GEN
     224                 :            : **   perm= permutation vector (length li-1), indexing the rows of mat: easier
     225                 :            : **     to maintain perm than to copy rows. For columns we can do it directly
     226                 :            : **     using e.g. swap(mat[i], mat[j])
     227                 :            : **   k0 = integer. The k0 first lines of mat are dense, the others are sparse.
     228                 :            : ** Output: cf ASCII art in the function body
     229                 :            : **
     230                 :            : ** row permutations applied to perm
     231                 :            : ** column operations applied to C. IN PLACE
     232                 :            : **/
     233                 :            : GEN
     234                 :       6404 : hnfspec_i(GEN mat0, GEN perm, GEN* ptdep, GEN* ptB, GEN* ptC, long k0)
     235                 :            : {
     236                 :            :   pari_sp av;
     237                 :            :   long co, n, s, nlze, lnz, nr, i, j, k, lk0, col, lig, *p;
     238                 :            :   GEN mat;
     239                 :            :   GEN p1, p2, matb, matbnew, vmax, matt, T, extramat, B, C, H, dep, permpro;
     240                 :       6404 :   const long li = lg(perm); /* = lgcols(mat0) */
     241                 :       6404 :   const long CO = lg(mat0);
     242                 :            : 
     243                 :       6404 :   n = 0; /* -Wall */
     244                 :            : 
     245                 :       6404 :   C = *ptC; co = CO;
     246 [ +  + ][ -  + ]:       6404 :   if (co > 300 && co > 1.5 * li)
     247                 :            :   { /* treat the rest at the end */
     248                 :          0 :     co = (long)(1.2 * li);
     249                 :          0 :     setlg(C, co);
     250                 :            :   }
     251                 :            : 
     252         [ -  + ]:       6404 :   if (DEBUGLEVEL>5)
     253                 :            :   {
     254                 :          0 :     err_printf("Entering hnfspec\n");
     255                 :          0 :     p_mat(mat0,perm,0);
     256                 :            :   }
     257                 :       6404 :   matt = cgetg(co, t_MAT); /* dense part of mat (top) */
     258                 :       6404 :   mat = cgetg(co, t_MAT);
     259         [ +  + ]:     411431 :   for (j = 1; j < co; j++)
     260                 :            :   {
     261                 :     405027 :     GEN matj = col_dup(li, gel(mat0,j));
     262                 :     405027 :     p1 = cgetg(k0+1,t_COL); gel(matt,j) = p1; gel(mat,j) = matj;
     263         [ +  + ]:    2938991 :     for (i=1; i<=k0; i++) gel(p1,i) = stoi(matj[perm[i]]);
     264                 :            :   }
     265                 :       6404 :   av = avma;
     266                 :            : 
     267                 :       6404 :   i = lig = li-1; col = co-1; lk0 = k0;
     268 [ +  + ][ +  - ]:       6404 :   T = (k0 || (lg(C) > 1 && lgcols(C) > 1))? matid(col): NULL;
                 [ +  + ]
     269                 :            :   /* Look for lines with a single non-0 entry, equal to 1 in absolute value */
     270 [ +  + ][ +  + ]:    2069558 :   while (i > lk0 && col)
     271      [ +  +  + ]:    2063154 :     switch( count(mat,perm[i],col,&n) )
     272                 :            :     {
     273                 :            :       case 0: /* move zero lines between k0+1 and lk0 */
     274                 :      15106 :         lk0++; lswap(perm[i], perm[lk0]);
     275                 :      15106 :         i = lig; continue;
     276                 :            : 
     277                 :            :       case 1: /* move trivial generator between lig+1 and li */
     278                 :      99954 :         lswap(perm[i], perm[lig]);
     279         [ +  + ]:      99954 :         if (T) swap(gel(T,n), gel(T,col));
     280                 :      99954 :         swap(gel(mat,n), gel(mat,col)); p = gel(mat,col);
     281         [ +  + ]:      99954 :         if (p[perm[lig]] < 0) /* = -1 */
     282                 :            :         { /* convert relation -g = 0 to g = 0 */
     283         [ +  + ]:    5713645 :           for (i=lk0+1; i<lig; i++) p[perm[i]] = -p[perm[i]];
     284         [ +  - ]:      86121 :           if (T)
     285                 :            :           {
     286                 :      86121 :             p1 = gel(T,col);
     287                 :      86121 :             for (i=1; ; i++) /* T is a permuted identity: single non-0 entry */
     288         [ +  + ]:    5480013 :               if (signe(gel(p1,i))) { togglesign_safe(&gel(p1,i)); break; }
     289                 :    5393892 :           }
     290                 :            :         }
     291                 :      99954 :         lig--; col--; i = lig; continue;
     292                 :            : 
     293                 :    1948094 :       default: i--;
     294                 :            :     }
     295         [ -  + ]:       6404 :   if (DEBUGLEVEL>5) { err_printf("    after phase1:\n"); p_mat(mat,perm,0); }
     296                 :            : 
     297                 :            : #define absmax(s,z) {long _z; _z = labs(z); if (_z > s) s = _z;}
     298                 :            :   /* Get rid of all lines containing only 0 and +/- 1, keeping track of column
     299                 :            :    * operations in T. Leave the rows 1..lk0 alone [up to k0, coefficient
     300                 :            :    * explosion, between k0+1 and lk0, row is 0] */
     301                 :       6404 :   s = 0;
     302 [ +  + ][ +  + ]:     186986 :   while (lig > lk0 && col && s < (long)(HIGHBIT>>1))
                 [ +  - ]
     303                 :            :   {
     304         [ +  + ]:     646884 :     for (i=lig; i>lk0; i--)
     305         [ +  + ]:     640844 :       if (count(mat,perm[i],col,&n) > 0) break;
     306         [ +  + ]:     186622 :     if (i == lk0) break;
     307                 :            : 
     308                 :            :     /* only 0, +/- 1 entries, at least 2 of them non-zero */
     309                 :     180582 :     lswap(perm[i], perm[lig]);
     310                 :     180582 :     swap(gel(mat,n), gel(mat,col)); p = gel(mat,col);
     311         [ +  - ]:     180582 :     if (T) swap(gel(T,n), gel(T,col));
     312         [ +  + ]:     180582 :     if (p[perm[lig]] < 0)
     313                 :            :     {
     314         [ +  + ]:    4703978 :       for (i=lk0+1; i<=lig; i++) p[perm[i]] = -p[perm[i]];
     315         [ +  - ]:     130945 :       if (T) ZV_togglesign(gel(T,col));
     316                 :            :     }
     317         [ +  + ]:    8753353 :     for (j=1; j<col; j++)
     318                 :            :     {
     319                 :    8572771 :       GEN matj = gel(mat,j);
     320                 :            :       long t;
     321         [ +  + ]:    8572771 :       if (! (t = matj[perm[lig]]) ) continue;
     322         [ +  + ]:     528623 :       if (t == 1) {
     323 [ +  + ][ +  + ]:    8302148 :         for (i=lk0+1; i<=lig; i++) absmax(s, matj[perm[i]] -= p[perm[i]]);
     324                 :            :       }
     325                 :            :       else { /* t = -1 */
     326 [ +  + ][ +  + ]:    8788987 :         for (i=lk0+1; i<=lig; i++) absmax(s, matj[perm[i]] += p[perm[i]]);
     327                 :            :       }
     328         [ +  - ]:     528623 :       if (T) ZC_lincomb1_inplace(gel(T,j), gel(T,col), stoi(-t));
     329                 :            :     }
     330                 :     180582 :     lig--; col--;
     331         [ -  + ]:     180582 :     if (gc_needed(av,3))
     332                 :            :     {
     333         [ #  # ]:          0 :       if(DEBUGMEM>1) pari_warn(warnmem,"hnfspec[1]");
     334         [ #  # ]:          0 :       if (T) T = gerepilecopy(av, T); else avma = av;
     335                 :            :     }
     336                 :            :   }
     337                 :            :   /* As above with lines containing a +/- 1 (no other assumption).
     338                 :            :    * Stop when single precision becomes dangerous */
     339                 :       6404 :   vmax = cgetg(co,t_VECSMALL);
     340         [ +  + ]:     130895 :   for (j=1; j<=col; j++)
     341                 :            :   {
     342                 :     124491 :     GEN matj = gel(mat,j);
     343 [ +  + ][ +  + ]:    1462534 :     for (s=0, i=lk0+1; i<=lig; i++) absmax(s, matj[i]);
     344                 :     124491 :     vmax[j] = s;
     345                 :            :   }
     346 [ +  + ][ +  + ]:      49199 :   while (lig > lk0 && col)
     347                 :            :   {
     348         [ +  + ]:      72639 :     for (i=lig; i>lk0; i--)
     349         [ +  + ]:      69021 :       if ( (n = count2(mat,perm[i],col)) ) break;
     350         [ +  + ]:      46413 :     if (i == lk0) break;
     351                 :            : 
     352                 :      42795 :     lswap(vmax[n], vmax[col]);
     353                 :      42795 :     lswap(perm[i], perm[lig]);
     354                 :      42795 :     swap(gel(mat,n), gel(mat,col)); p = gel(mat,col);
     355         [ +  - ]:      42795 :     if (T) swap(gel(T,n), gel(T,col));
     356         [ +  + ]:      42795 :     if (p[perm[lig]] < 0)
     357                 :            :     {
     358         [ +  + ]:     182203 :       for (i=lk0+1; i<=lig; i++) p[perm[i]] = -p[perm[i]];
     359         [ +  - ]:      21738 :       if (T) ZV_togglesign(gel(T,col));
     360                 :            :     }
     361         [ +  + ]:     895532 :     for (j=1; j<col; j++)
     362                 :            :     {
     363                 :     852737 :       GEN matj = gel(mat,j);
     364                 :            :       long t;
     365         [ +  + ]:     852737 :       if (! (t = matj[perm[lig]]) ) continue;
     366 [ +  + ][ -  + ]:     497933 :       if (vmax[col] && (ulong)labs(t) >= (HIGHBIT-vmax[j]) / vmax[col])
     367                 :          0 :         goto END2;
     368                 :            : 
     369 [ +  + ][ +  + ]:    5251373 :       for (s=0, i=lk0+1; i<=lig; i++) absmax(s, matj[perm[i]] -= t*p[perm[i]]);
     370                 :     497933 :       vmax[j] = s;
     371         [ +  - ]:     497933 :       if (T) ZC_lincomb1_inplace(gel(T,j), gel(T,col), stoi(-t));
     372                 :            :     }
     373                 :      42795 :     lig--; col--;
     374         [ -  + ]:      42795 :     if (gc_needed(av,3))
     375                 :            :     {
     376         [ #  # ]:          0 :       if(DEBUGMEM>1) pari_warn(warnmem,"hnfspec[2]");
     377         [ #  # ]:          0 :       gerepileall(av, T? 2: 1, &vmax, &T);
     378                 :            :     }
     379                 :            :   }
     380                 :            : 
     381                 :            : END2: /* clean up mat: remove everything to the right of the 1s on diagonal */
     382                 :            :   /* go multiprecision first */
     383                 :       6404 :   matb = cgetg(co,t_MAT); /* bottom part (complement of matt) */
     384         [ +  + ]:     411431 :   for (j=1; j<co; j++)
     385                 :            :   {
     386                 :     405027 :     GEN matj = gel(mat,j);
     387                 :     405027 :     p1 = cgetg(li-k0,t_COL); gel(matb,j) = p1;
     388                 :     405027 :     p1 -= k0;
     389         [ +  + ]:   32296786 :     for (i=k0+1; i<li; i++) gel(p1,i) = stoi(matj[perm[i]]);
     390                 :            :   }
     391         [ -  + ]:       6404 :   if (DEBUGLEVEL>5)
     392                 :            :   {
     393                 :          0 :     err_printf("    after phase2:\n");
     394                 :          0 :     p_mat(mat,perm,lk0);
     395                 :            :   }
     396         [ +  + ]:     323513 :   for (i=li-2; i>lig; i--)
     397                 :            :   {
     398                 :     317109 :     long h, i0 = i - k0, k = i + co-li;
     399                 :     317109 :     GEN Bk = gel(matb,k);
     400         [ +  + ]:   12490104 :     for (j=k+1; j<co; j++)
     401                 :            :     {
     402                 :   12172995 :       GEN Bj = gel(matb,j), v = gel(Bj,i0);
     403         [ +  + ]:   12172995 :       s = signe(v); if (!s) continue;
     404                 :            : 
     405                 :    2507584 :       gel(Bj,i0) = gen_0;
     406         [ +  + ]:    2507584 :       if (is_pm1(v))
     407                 :            :       {
     408         [ +  + ]:    1401587 :         if (s > 0) /* v = 1 */
     409         [ +  + ]:   15219626 :         { for (h=1; h<i0; h++) gel(Bj,h) = subii(gel(Bj,h), gel(Bk,h)); }
     410                 :            :         else /* v = -1 */
     411         [ +  + ]:   14678295 :         { for (h=1; h<i0; h++) gel(Bj,h) = addii(gel(Bj,h), gel(Bk,h)); }
     412                 :            :       }
     413                 :            :       else {
     414         [ +  + ]:   13102363 :         for (h=1; h<i0; h++) gel(Bj,h) = subii(gel(Bj,h), mulii(v,gel(Bk,h)));
     415                 :            :       }
     416         [ +  - ]:    2507584 :       if (T) ZC_lincomb1_inplace(gel(T,j), gel(T,k), negi(v));
     417         [ +  + ]:    2507584 :       if (gc_needed(av,3))
     418                 :            :       {
     419         [ -  + ]:         26 :         if(DEBUGMEM>1) pari_warn(warnmem,"hnfspec[3], (i,j) = %ld,%ld", i,j);
     420         [ +  + ]:       9744 :         for (h=1; h<co; h++) setlg(matb[h], i0+1); /* bottom can be forgotten */
     421         [ +  - ]:         26 :         gerepileall(av, T? 2: 1, &matb, &T);
     422                 :         26 :         Bk = gel(matb,k);
     423                 :            :       }
     424                 :            :     }
     425                 :            :   }
     426         [ +  + ]:     411431 :   for (j=1; j<co; j++) setlg(matb[j], lig-k0+1); /* bottom can be forgotten */
     427         [ +  + ]:       6404 :   gerepileall(av, T? 2: 1, &matb, &T);
     428         [ -  + ]:       6404 :   if (DEBUGLEVEL>5) err_printf("    matb cleaned up (using Id block)\n");
     429                 :            : 
     430                 :       6404 :   nlze = lk0 - k0;  /* # of 0 rows */
     431                 :       6404 :   lnz = lig-nlze+1; /* 1 + # of non-0 rows (!= 0...0 1 0 ... 0) */
     432         [ +  + ]:       6404 :   if (T) matt = ZM_mul(matt,T); /* update top rows */
     433                 :       6404 :   extramat = cgetg(col+1,t_MAT); /* = new C minus the 0 rows */
     434         [ +  + ]:      88100 :   for (j=1; j<=col; j++)
     435                 :            :   {
     436                 :      81696 :     GEN z = gel(matt,j);
     437                 :      81696 :     GEN t = (gel(matb,j)) + nlze - k0;
     438                 :      81696 :     p2=cgetg(lnz,t_COL); gel(extramat,j) = p2;
     439         [ +  + ]:     570573 :     for (i=1; i<=k0; i++) gel(p2,i) = gel(z,i); /* top k0 rows */
     440         [ +  + ]:     248765 :     for (   ; i<lnz; i++) gel(p2,i) = gel(t,i); /* other non-0 rows */
     441                 :            :   }
     442         [ +  + ]:       6404 :   if (!col) {
     443                 :        140 :     permpro = identity_perm(lnz);
     444                 :        140 :     nr = lnz;
     445                 :            :   }
     446                 :            :   else
     447                 :       6264 :     permpro = ZM_imagecomplspec(extramat, &nr);
     448                 :            :   /* lnz = lg(permpro) */
     449         [ +  + ]:       6404 :   if (nlze)
     450                 :            :   { /* put the nlze 0 rows (trivial generators) at the top */
     451                 :       2961 :     p1 = new_chunk(lk0+1);
     452         [ +  + ]:      18067 :     for (i=1; i<=nlze; i++) p1[i] = perm[i + k0];
     453         [ +  + ]:      22092 :     for (   ; i<=lk0; i++)  p1[i] = perm[i - nlze];
     454         [ +  + ]:      37198 :     for (i=1; i<=lk0; i++)  perm[i] = p1[i];
     455                 :            :   }
     456                 :            :   /* sort other rows according to permpro (nr redundant generators first) */
     457                 :       6404 :   p1 = new_chunk(lnz); p2 = perm + nlze;
     458         [ +  + ]:      50026 :   for (i=1; i<lnz; i++) p1[i] = p2[permpro[i]];
     459         [ +  + ]:      50026 :   for (i=1; i<lnz; i++) p2[i] = p1[i];
     460                 :            :   /* perm indexes the rows of mat
     461                 :            :    *   |_0__|__redund__|__dense__|__too big__|_____done______|
     462                 :            :    *   0  nlze                              lig             li
     463                 :            :    *         \___nr___/ \___k0__/
     464                 :            :    *         \____________lnz ______________/
     465                 :            :    *
     466                 :            :    *               col   co
     467                 :            :    *       [dep     |     ]
     468                 :            :    *    i0 [--------|  B  ] (i0 = nlze + nr)
     469                 :            :    *       [matbnew |     ] matbnew has maximal rank = lnz-1 - nr
     470                 :            :    * mat = [--------|-----] lig
     471                 :            :    *       [   0    | Id  ]
     472                 :            :    *       [        |     ] li */
     473                 :            : 
     474                 :       6404 :   matbnew = cgetg(col+1,t_MAT); /* dense+toobig, maximal rank. For hnffinal */
     475                 :       6404 :   dep    = cgetg(col+1,t_MAT); /* rows dependent from the ones in matbnew */
     476         [ +  + ]:      88100 :   for (j=1; j<=col; j++)
     477                 :            :   {
     478                 :      81696 :     GEN z = gel(extramat,j);
     479                 :      81696 :     p1 = cgetg(nlze+nr+1,t_COL); gel(dep,j) = p1;
     480                 :      81696 :     p2 = cgetg(lnz-nr,t_COL); gel(matbnew,j) = p2;
     481         [ +  + ]:     281294 :     for (i=1; i<=nlze; i++) gel(p1,i) = gen_0;
     482         [ +  + ]:     107960 :     p1 += nlze; for (i=1; i<=nr; i++) p1[i] = z[permpro[i]];
     483         [ +  + ]:     711378 :     p2 -= nr;   for (   ; i<lnz; i++) p2[i] = z[permpro[i]];
     484                 :            :   }
     485                 :            : 
     486                 :            :   /* redundant generators in terms of the genuine generators
     487                 :            :    * (x_i) = - (g_i) B */
     488                 :       6404 :   B = cgetg(co-col,t_MAT);
     489         [ +  + ]:     329735 :   for (j=col+1; j<co; j++)
     490                 :            :   {
     491                 :     323331 :     GEN y = gel(matt,j);
     492                 :     323331 :     GEN z = gel(matb,j);
     493                 :     323331 :     p1=cgetg(lig+1,t_COL); gel(B,j-col) = p1;
     494         [ +  + ]:    1372890 :     for (i=1; i<=nlze; i++) gel(p1,i) = gel(z,i);
     495                 :     323331 :     p1 += nlze; z += nlze-k0;
     496         [ +  + ]:    3155568 :     for (k=1; k<lnz; k++)
     497                 :            :     {
     498                 :    2832237 :       i = permpro[k];
     499         [ +  + ]:    2832237 :       gel(p1,k) = (i <= k0)? gel(y,i): gel(z,i);
     500                 :            :     }
     501                 :            :   }
     502 [ +  + ][ +  + ]:       6404 :   if (T) C = typ(C)==t_MAT? RgM_mul(C,T): RgV_RgM_mul(C,T);
     503                 :       6404 :   gerepileall(av, 4, &matbnew, &B, &dep, &C);
     504                 :       6404 :   *ptdep = dep;
     505                 :       6404 :   *ptB = B;
     506                 :       6404 :   H = hnffinal(matbnew, perm, ptdep, ptB, &C);
     507         [ -  + ]:       6404 :   if (CO > co)
     508                 :            :   { /* treat the rest, N cols at a time (hnflll slow otherwise) */
     509                 :          0 :     const long N = 300;
     510                 :          0 :     long a, L = CO - co, l = minss(L, N); /* L columns to add */
     511                 :          0 :     GEN CC = *ptC, m0 = mat0;
     512                 :          0 :     setlg(CC, CO); /* restore */
     513                 :          0 :     CC += co-1;
     514                 :          0 :     m0 += co-1;
     515                 :          0 :     for (a = l;;)
     516                 :            :     {
     517                 :            :       GEN MAT, emb;
     518                 :          0 :       gerepileall(av, 4, &H,&C,ptB,ptdep);
     519                 :          0 :       MAT = cgetg(l + 1, t_MAT);
     520                 :          0 :       emb = cgetg(l + 1, typ(C));
     521         [ #  # ]:          0 :       for (j = 1 ; j <= l; j++)
     522                 :            :       {
     523                 :          0 :         gel(MAT,j) = gel(m0,j);
     524                 :          0 :         emb[j] = CC[j];
     525                 :            :       }
     526                 :          0 :       H = hnfadd_i(H, perm, ptdep, ptB, &C, MAT, emb);
     527         [ #  # ]:          0 :       if (a == L) break;
     528                 :          0 :       CC += l;
     529                 :          0 :       m0 += l;
     530         [ #  # ]:          0 :       a += l; if (a > L) { l = L - (a - l); a = L; }
     531                 :          0 :     }
     532                 :            :   }
     533                 :       6404 :   *ptC = C; return H;
     534                 :            : }
     535                 :            : 
     536                 :            : GEN
     537                 :        133 : hnfspec(GEN mat, GEN perm, GEN* ptdep, GEN* ptB, GEN* ptC, long k0)
     538                 :            : {
     539                 :        133 :   pari_sp av = avma;
     540                 :        133 :   GEN H = hnfspec_i(mat, perm, ptdep, ptB, ptC, k0);
     541                 :        133 :   gerepileall(av, 4, ptC, ptdep, ptB, &H); return H;
     542                 :            : }
     543                 :            : 
     544                 :            : /* HNF reduce x, apply same transforms to C */
     545                 :            : GEN
     546                 :        133 : mathnfspec(GEN x, GEN *ptperm, GEN *ptdep, GEN *ptB, GEN *ptC)
     547                 :            : {
     548                 :        133 :   long i,j,k,ly,lx = lg(x);
     549                 :            :   GEN z, perm;
     550         [ -  + ]:        133 :   if (lx == 1) return cgetg(1, t_MAT);
     551                 :        133 :   ly = lgcols(x);
     552                 :        133 :   *ptperm = perm = identity_perm(ly-1);
     553                 :        133 :   z = cgetg(lx,t_MAT);
     554         [ +  + ]:       1869 :   for (i=1; i<lx; i++)
     555                 :            :   {
     556                 :       1736 :     GEN C = cgetg(ly,t_COL), D = gel(x,i);
     557                 :       1736 :     gel(z,i) = C;
     558         [ +  + ]:      85064 :     for (j=1; j<ly; j++)
     559                 :            :     {
     560                 :      83328 :       GEN d = gel(D,j);
     561         [ -  + ]:      83328 :       if (is_bigint(d)) goto TOOLARGE;
     562                 :      83328 :       C[j] = itos(d);
     563                 :            :     }
     564                 :            :   }
     565                 :            :   /*  [ dep |     ]
     566                 :            :    *  [-----|  B  ]
     567                 :            :    *  [  H  |     ]
     568                 :            :    *  [-----|-----]
     569                 :            :    *  [  0  | Id  ] */
     570                 :        133 :   return hnfspec(z,perm, ptdep, ptB, ptC, 0);
     571                 :            : 
     572                 :            : TOOLARGE:
     573 [ #  # ][ #  # ]:          0 :   if (lg(*ptC) > 1 && lgcols(*ptC) > 1)
     574                 :          0 :     pari_err_IMPL("mathnfspec with large entries");
     575                 :          0 :   x = ZM_hnf(x); lx = lg(x); j = ly; k = 0;
     576         [ #  # ]:          0 :   for (i=1; i<ly; i++)
     577                 :            :   {
     578         [ #  # ]:          0 :     if (equali1(gcoeff(x,i,i + lx-ly)))
     579                 :          0 :       perm[--j] = i;
     580                 :            :     else
     581                 :          0 :       perm[++k] = i;
     582                 :            :   }
     583                 :          0 :   setlg(perm,k+1);
     584                 :          0 :   x = rowpermute(x, perm); /* upper part */
     585                 :          0 :   setlg(perm,ly);
     586                 :          0 :   *ptB = vecslice(x, j+lx-ly, lx-1);
     587                 :          0 :   setlg(x, j);
     588                 :          0 :   *ptdep = rowslice(x, 1, lx-ly);
     589                 :        133 :   return rowslice(x, lx-ly+1, k); /* H */
     590                 :            : }
     591                 :            : 
     592                 :            : /* add new relations to a matrix treated by hnfspec (extramat / extraC) */
     593                 :            : GEN
     594                 :       9119 : hnfadd_i(GEN H, GEN perm, GEN* ptdep, GEN* ptB, GEN* ptC, /* cf hnfspec */
     595                 :            :        GEN extramat,GEN extraC)
     596                 :            : {
     597                 :       9119 :   GEN matb, extratop, Cnew, permpro, B = *ptB, C = *ptC, dep = *ptdep;
     598                 :            :   long i, lH, lB, li, lig, co, col, nlze;
     599                 :            : 
     600         [ -  + ]:       9119 :   if (lg(extramat) == 1) return H;
     601                 :       9119 :   co = lg(C)-1;
     602                 :       9119 :   lH = lg(H)-1;
     603                 :       9119 :   lB = lg(B)-1;
     604                 :       9119 :   li = lg(perm)-1;
     605                 :       9119 :   lig = li - lB;
     606                 :       9119 :   col = co - lB;
     607                 :       9119 :   nlze = lig - lH;
     608                 :            : 
     609                 :            :  /*               col    co
     610                 :            :   *       [ 0 |dep |     ]
     611                 :            :   *  nlze [--------|  B  ]
     612                 :            :   *       [ 0 | H  |     ]
     613                 :            :   *       [--------|-----] lig
     614                 :            :   *       [   0    | Id  ]
     615                 :            :   *       [        |     ] li */
     616                 :       9119 :   extratop = zm_to_ZM( rowslicepermute(extramat, perm, 1, lig) );
     617         [ +  - ]:       9119 :   if (li != lig)
     618                 :            :   { /* zero out bottom part, using the Id block */
     619                 :       9119 :     GEN A = vecslice(C, col+1, co);
     620                 :       9119 :     GEN c = rowslicepermute(extramat, perm, lig+1, li);
     621         [ +  + ]:       9119 :     extraC   = gsub(extraC, typ(A)==t_MAT? RgM_zm_mul(A, c): RgV_zm_mul(A,c));
     622                 :       9119 :     extratop = ZM_sub(extratop, ZM_zm_mul(B, c));
     623                 :            :   }
     624                 :            : 
     625                 :       9119 :   extramat = shallowconcat(extratop, vconcat(dep, H));
     626                 :       9119 :   Cnew     = shallowconcat(extraC, vecslice(C, col-lH+1, co));
     627         [ -  + ]:       9119 :   if (DEBUGLEVEL>5) err_printf("    1st phase done\n");
     628                 :       9119 :   permpro = ZM_imagecomplspec(extramat, &nlze);
     629                 :       9119 :   extramat = rowpermute(extramat, permpro);
     630                 :       9119 :   *ptB     = rowpermute(B,        permpro);
     631                 :       9119 :   permpro = vecsmallpermute(perm, permpro);
     632         [ +  + ]:     102667 :   for (i=1; i<=lig; i++) perm[i] = permpro[i]; /* perm o= permpro */
     633                 :            : 
     634                 :       9119 :   *ptdep  = rowslice(extramat, 1, nlze);
     635                 :       9119 :   matb    = rowslice(extramat, nlze+1, lig);
     636         [ -  + ]:       9119 :   if (DEBUGLEVEL>5) err_printf("    2nd phase done\n");
     637                 :       9119 :   H = hnffinal(matb,perm,ptdep,ptB,&Cnew);
     638                 :       9119 :   *ptC = shallowconcat(vecslice(C, 1, col-lH), Cnew);
     639                 :       9119 :   return H;
     640                 :            : }
     641                 :            : 
     642                 :            : GEN
     643                 :          0 : hnfadd(GEN H, GEN perm, GEN* ptdep, GEN* ptB, GEN* ptC, /* cf hnfspec */
     644                 :            :        GEN extramat,GEN extraC)
     645                 :            : {
     646                 :          0 :   pari_sp av = avma;
     647                 :          0 :   H = hnfadd_i(H, perm, ptdep, ptB, ptC, ZM_to_zm(extramat), extraC);
     648                 :          0 :   gerepileall(av, 4, ptC, ptdep, ptB, &H); return H;
     649                 :            : }
     650                 :            : 
     651                 :            : /* zero aj = Aij (!= 0)  using  ak = Aik (maybe 0), via linear combination of
     652                 :            :  * A[j] and A[k] of determinant 1. If U != NULL, likewise update its columns */
     653                 :            : static void
     654                 :   20714496 : ZC_elem(GEN aj, GEN ak, GEN A, GEN U, long j, long k)
     655                 :            : {
     656                 :            :   GEN p1,u,v,d;
     657                 :            : 
     658         [ +  + ]:   20714496 :   if (!signe(ak)) {
     659                 :    4302392 :     swap(gel(A,j), gel(A,k));
     660         [ -  + ]:    4302392 :     if (U) swap(gel(U,j), gel(U,k));
     661                 :            :     return;
     662                 :            :   }
     663                 :   16412104 :   d = bezout(aj,ak,&u,&v);
     664                 :            :   /* frequent special case (u,v) = (1,0) or (0,1) */
     665         [ +  + ]:   16412104 :   if (!signe(u))
     666                 :            :   { /* ak | aj */
     667                 :   11504207 :     p1 = diviiexact(aj,ak); togglesign(p1);
     668                 :   11504207 :     ZC_lincomb1_inplace(gel(A,j), gel(A,k), p1);
     669         [ +  + ]:   11504207 :     if (U)
     670                 :     680180 :       ZC_lincomb1_inplace(gel(U,j), gel(U,k), p1);
     671                 :            :     return;
     672                 :            :   }
     673         [ +  + ]:    4907897 :   if (!signe(v))
     674                 :            :   { /* aj | ak */
     675                 :    3017737 :     p1 = diviiexact(ak,aj); togglesign(p1);
     676                 :    3017737 :     ZC_lincomb1_inplace(gel(A,k), gel(A,j), p1);
     677                 :    3017737 :     swap(gel(A,j), gel(A,k));
     678         [ +  + ]:    3017737 :     if (U) {
     679                 :     100983 :       ZC_lincomb1_inplace(gel(U,k), gel(U,j), p1);
     680                 :     100983 :       swap(gel(U,j), gel(U,k));
     681                 :            :     }
     682                 :            :     return;
     683                 :            :   }
     684                 :            : 
     685         [ +  + ]:    1890160 :   if (!is_pm1(d)) { aj = diviiexact(aj, d); ak = diviiexact(ak, d); }
     686                 :    1890160 :   p1 = gel(A,k); aj = negi(aj); /* NOT togglesign */
     687                 :    1890160 :   gel(A,k) = ZC_lincomb(u,v, gel(A,j),p1);
     688                 :    1890160 :   gel(A,j) = ZC_lincomb(aj,ak, p1,gel(A,j));
     689         [ +  + ]:    1890160 :   if (U)
     690                 :            :   {
     691                 :     135003 :     p1 = gel(U,k);
     692                 :     135003 :     gel(U,k) = ZC_lincomb(u,v, gel(U,j),p1);
     693                 :   20714496 :     gel(U,j) = ZC_lincomb(aj,ak, p1,gel(U,j));
     694                 :            :   }
     695                 :            : }
     696                 :            : 
     697                 :            : INLINE int
     698 [ +  + ][ +  - ]:        931 : is_RgX(GEN a, long v) { return typ(a) == t_POL && varn(a)==v; }
     699                 :            : /* set u,v such that au + bv = gcd(a,b), divide a,b by the gcd */
     700                 :            : static GEN
     701                 :        266 : gbezout_step(GEN *pa, GEN *pb, GEN *pu, GEN *pv, long vx)
     702                 :            : {
     703                 :        266 :   GEN a = *pa, b = *pb, d;
     704         [ -  + ]:        266 :   if (gequal0(a))
     705                 :            :   {
     706                 :          0 :     *pa = gen_0; *pu = gen_0;
     707                 :          0 :     *pb = gen_1; *pv = gen_1; return b;
     708                 :            :   }
     709         [ +  + ]:        266 :   a = is_RgX(a,vx)? RgX_renormalize(a): scalarpol(a, vx);
     710         [ +  + ]:        266 :   b = is_RgX(b,vx)? RgX_renormalize(b): scalarpol(b, vx);
     711                 :        266 :   d = RgX_extgcd(a,b, pu,pv);
     712         [ +  + ]:        266 :   if (degpol(d)) { a = RgX_div(a, d); b = RgX_div(b, d); }
     713 [ -  + ][ #  # ]:        154 :   else if (typ(gel(d,2)) == t_REAL && lg(gel(d,2)) <= 3)
     714                 :            : #if 1
     715                 :            :   { /* possible accuracy problem */
     716                 :          0 :     GEN D = RgX_gcd_simple(a,b);
     717         [ #  # ]:          0 :     if (degpol(D)) {
     718                 :          0 :       D = RgX_Rg_div(D, leading_term(D));
     719                 :          0 :       a = RgX_div(a, D);
     720                 :          0 :       b = RgX_div(b, D);
     721                 :          0 :       d = RgX_extgcd(a,b, pu,pv); /* retry now */
     722                 :          0 :       d = RgX_mul(d, D);
     723                 :            :     }
     724                 :            :   }
     725                 :            : #else
     726                 :            :   { /* less stable */
     727                 :            :     d = RgX_extgcd_simple(a,b, pu,pv);
     728                 :            :     if (degpol(d)) { a = RgX_div(a, d); b = RgX_div(b, d); }
     729                 :            :   }
     730                 :            : #endif
     731                 :        266 :   *pa = a;
     732                 :        266 :   *pb = b; return d;
     733                 :            : }
     734                 :            : static GEN
     735                 :     262752 : col_mul(GEN x, GEN c)
     736                 :            : {
     737         [ +  - ]:     262752 :   if (typ(x) == t_INT)
     738                 :            :   {
     739                 :     262752 :     long s = signe(x);
     740         [ +  + ]:     262752 :     if (!s) return NULL;
     741 [ +  + ][ +  + ]:     197190 :     if (is_pm1(x)) return (s > 0)? c: RgC_neg(c);
     742                 :            :   }
     743                 :     262752 :   return RgC_Rg_mul(c, x);
     744                 :            : }
     745                 :            : static void
     746                 :          0 : do_zero(GEN x)
     747                 :            : {
     748                 :          0 :   long i, lx = lg(x);
     749         [ #  # ]:          0 :   for (i=1; i<lx; i++) gel(x,i) = gen_0;
     750                 :          0 : }
     751                 :            : 
     752                 :            : /* (c1, c2) *= [u,-b; v,a] */
     753                 :            : static void
     754                 :      65688 : update(GEN u, GEN v, GEN a, GEN b, GEN *c1, GEN *c2)
     755                 :            : {
     756                 :            :   GEN p1,p2;
     757                 :            : 
     758                 :      65688 :   u = col_mul(u,*c1);
     759                 :      65688 :   v = col_mul(v,*c2);
     760 [ +  + ][ +  + ]:      65688 :   if (u) p1 = v? gadd(u,v): u;
     761                 :        140 :   else   p1 = v? v: NULL;
     762                 :            : 
     763                 :      65688 :   a = col_mul(a,*c2);
     764                 :      65688 :   b = col_mul(gneg_i(b),*c1);
     765 [ +  - ][ +  - ]:      65688 :   if (a) p2 = b? RgC_add(a,b): a;
     766                 :          0 :   else   p2 = b? b: NULL;
     767                 :            : 
     768         [ -  + ]:      65688 :   if (!p1) do_zero(*c1); else *c1 = p1;
     769         [ -  + ]:      65688 :   if (!p2) do_zero(*c2); else *c2 = p2;
     770                 :      65688 : }
     771                 :            : 
     772                 :            : /* zero aj = Aij (!= 0)  using  ak = Aik (maybe 0), via linear combination of
     773                 :            :  * A[j] and A[k] of determinant 1. If U != NULL, likewise update its columns */
     774                 :            : static void
     775                 :          7 : RgC_elem(GEN aj, GEN ak, GEN A, GEN V, long j, long k, long li, long vx)
     776                 :            : {
     777                 :          7 :   GEN u,v, d = gbezout_step(&aj, &ak, &u, &v, vx);
     778                 :            :   long l;
     779                 :            :   /* (A[,k], A[,j]) *= [v, -aj; u, ak ] */
     780         [ +  + ]:         14 :   for (l = 1; l < li; l++)
     781                 :            :   {
     782                 :          7 :     GEN t = gadd(gmul(u,gcoeff(A,l,j)), gmul(v,gcoeff(A,l,k)));
     783                 :          7 :     gcoeff(A,l,j) = gsub(gmul(ak,gcoeff(A,l,j)), gmul(aj,gcoeff(A,l,k)));
     784                 :          7 :     gcoeff(A,l,k) = t;
     785                 :            :   }
     786                 :          7 :   gcoeff(A,li,j) = gen_0;
     787                 :          7 :   gcoeff(A,li,k) = d;
     788         [ -  + ]:          7 :   if (V) update(v,u,ak,aj,(GEN*)(V+k),(GEN*)(V+j));
     789                 :          7 : }
     790                 :            : 
     791                 :            : /* reduce A[i,j] mod A[i,j0] for j=j0+1... via column operations */
     792                 :            : static void
     793                 :    1202872 : ZM_reduce(GEN A, GEN U, long i, long j0)
     794                 :            : {
     795                 :    1202872 :   long j, lA = lg(A);
     796                 :    1202872 :   GEN d = gcoeff(A,i,j0);
     797         [ +  + ]:    1202872 :   if (signe(d) < 0)
     798                 :            :   {
     799                 :       3907 :     ZV_neg_inplace(gel(A,j0));
     800         [ +  + ]:       3907 :     if (U) ZV_togglesign(gel(U,j0));
     801                 :       3907 :     d = gcoeff(A,i,j0);
     802                 :            :   }
     803         [ +  + ]:    3245971 :   for (j=j0+1; j<lA; j++)
     804                 :            :   {
     805                 :    2043099 :     GEN q = truedivii(gcoeff(A,i,j), d);
     806         [ +  + ]:    2043099 :     if (!signe(q)) continue;
     807                 :            : 
     808                 :     104878 :     togglesign(q);
     809                 :     104878 :     ZC_lincomb1_inplace(gel(A,j), gel(A,j0), q);
     810         [ +  + ]:     104878 :     if (U) ZC_lincomb1_inplace(gel(U,j), gel(U,j0), q);
     811                 :            :   }
     812                 :    1202872 : }
     813                 :            : 
     814                 :            : static GEN
     815                 :        168 : RgX_normalize_all(GEN T, GEN *pd)
     816                 :            : {
     817                 :        168 :   GEN d = leading_term(T);
     818 [ -  + ][ -  + ]:        168 :   while (gequal0(d) || ( typ(d) == t_REAL && lg(d) == 3
                 [ #  # ]
     819         [ #  # ]:          0 :                        && gexpo(T) - expo(d) > (long)BITS_IN_LONG)) {
     820                 :          0 :      T = normalizepol_lg(T, lg(T)-1);
     821         [ #  # ]:          0 :      if (!signe(T)) { *pd = gen_1; return T; }
     822                 :          0 :      d = leading_term(T);
     823                 :            :   }
     824                 :        168 :   *pd = d;
     825                 :        168 :   return RgX_Rg_div(T, d);
     826                 :            : }
     827                 :            : /* reduce A[i,j] mod A[i,j0] for j=j0+1... via column operations */
     828                 :            : static void
     829                 :         28 : RgM_reduce(GEN A, GEN U, long i, long j0, long vx)
     830                 :            : {
     831                 :         28 :   long j, lA = lg(A);
     832                 :         28 :   GEN d, T = gcoeff(A,i,j0);
     833         [ +  - ]:         28 :   if (is_RgX(T,vx)) {
     834                 :         28 :     T = RgX_normalize_all(T, &d);
     835         [ +  + ]:         28 :     if (degpol(T) == 0) { d = gel(T,2); T = gen_1; }
     836                 :            :   } else {
     837                 :          0 :     d = T; T = gen_1;
     838                 :            :   }
     839 [ +  + ][ -  + ]:         28 :   if (U && !gequal1(d)) gel(U,j0) = RgC_Rg_div(gel(U,j0), d);
     840                 :         28 :   gcoeff(A,i,j0) = T;
     841                 :            : 
     842         [ +  + ]:         35 :   for (j=j0+1; j<lA; j++)
     843                 :            :   {
     844                 :          7 :     GEN t = gcoeff(A,i,j), q;
     845         [ -  + ]:          7 :     if (gequal0(t)) continue;
     846         [ -  + ]:          7 :     if (T == gen_1)
     847                 :          0 :       q = t;
     848         [ +  - ]:          7 :     else if (is_RgX(t,vx))
     849                 :          7 :       q = RgX_div(t, T);
     850                 :          0 :     else continue;
     851                 :            : 
     852         [ +  - ]:          7 :     if (gequal0(q)) continue;
     853                 :          0 :     gel(A,j) = RgC_sub(gel(A,j), RgC_Rg_mul(gel(A,j0), q));
     854         [ #  # ]:          0 :     if (U) gel(U,j) = RgC_sub(gel(U,j), RgC_Rg_mul(gel(U,j0), q));
     855                 :            :   }
     856                 :         28 : }
     857                 :            : 
     858                 :            : /* A,B square integral in upper HNF, of the same dimension > 0. Return Au
     859                 :            :  * in Z^n (v in Z^n not computed), such that Au + Bv = [1, 0, ..., 0] */
     860                 :            : GEN
     861                 :     178795 : hnfmerge_get_1(GEN A, GEN B)
     862                 :            : {
     863                 :     178795 :   pari_sp av = avma;
     864                 :     178795 :   long j, k, c, l = lg(A), lb;
     865                 :     178795 :   GEN b, t, U = cgetg(l + 1, t_MAT), C = cgetg(l + 1, t_VEC);
     866                 :            : 
     867                 :     178795 :   t = NULL; /* -Wall */
     868                 :     178795 :   b = gcoeff(B,1,1); lb = lgefint(b);
     869         [ -  + ]:     178795 :   if (!signe(b)) {
     870         [ #  # ]:          0 :     if (!is_pm1(gcoeff(A,1,1))) return NULL;
     871                 :          0 :     return scalarcol_shallow(gen_1, l-1);
     872                 :            :   }
     873         [ +  + ]:     432763 :   for (j = 1; j < l; j++)
     874                 :            :   {
     875                 :     432756 :     c = j+1;
     876                 :     432756 :     gel(U,j) = col_ei(l-1, j);
     877                 :     432756 :     gel(U,c) = zerocol(l-1); /* dummy */
     878                 :     432756 :     gel(C,j) = vecslice(gel(A,j), 1,j);
     879                 :     432756 :     gel(C,c) = vecslice(gel(B,j), 1,j);
     880         [ +  + ]:    1277547 :     for (k = j; k > 0; k--)
     881                 :            :     {
     882                 :     844791 :       t = gcoeff(C,k,c);
     883         [ +  + ]:     844791 :       if (gequal0(t)) continue;
     884                 :     736571 :       setlg(C[c], k+1);
     885                 :     736571 :       ZC_elem(t, gcoeff(C,k,k), C, U, c, k);
     886         [ -  + ]:     736571 :       if (lgefint(gcoeff(C,k,k)) > lb) gel(C,k) = FpC_red(gel(C,k), b);
     887         [ +  + ]:     736571 :       if (j > 4)
     888                 :            :       {
     889                 :      52122 :         GEN u = gel(U,k);
     890                 :            :         long h;
     891         [ +  + ]:     639170 :         for (h=1; h<l; h++)
     892         [ +  + ]:     587048 :           if (lgefint(gel(u,h)) > lb) gel(u,h) = remii(gel(u,h), b);
     893                 :            :       }
     894                 :            :     }
     895         [ +  + ]:     432756 :     if (j == 1)
     896                 :     178795 :       t = gcoeff(C,1,1);
     897                 :            :     else
     898                 :            :     {
     899                 :            :       GEN u;
     900                 :     253961 :       t = bezout(gcoeff(C,1,1), b, &u, NULL); /* >= 0 */
     901 [ +  + ][ +  + ]:     253961 :       if (signe(u) && !equali1(u)) gel(U,1) = ZC_Z_mul(gel(U,1), u);
     902                 :     253961 :       gcoeff(C,1,1) = t;
     903                 :            :     }
     904         [ +  + ]:     432756 :     if (equali1(t)) break;
     905                 :            :   }
     906         [ +  + ]:     178795 :   if (j >= l) return NULL;
     907                 :     178795 :   return gerepileupto(av, ZM_ZC_mul(A,gel(U,1)));
     908                 :            : }
     909                 :            : 
     910                 :            : /* remove the first r columns */
     911                 :            : static void
     912                 :      17595 : remove_0cols(long r, GEN *pA, GEN *pB, long remove)
     913                 :            : {
     914                 :      17595 :   GEN A = *pA, B = *pB;
     915                 :      17595 :   long l = lg(A);
     916                 :      17595 :   A += r; A[0] = evaltyp(t_MAT) | evallg(l-r);
     917 [ +  + ][ +  + ]:      17595 :   if (B && remove == 2) { B += r; B[0] = A[0]; }
     918                 :      17595 :   *pA = A; *pB = B;
     919                 :      17595 : }
     920                 :            : 
     921                 :            : /* Inefficient compared to hnfall. 'remove' = throw away lin.dep columns */
     922                 :            : static GEN
     923                 :       6916 : hnf_i(GEN A, int remove)
     924                 :            : {
     925                 :       6916 :   pari_sp av0 = avma, av;
     926                 :            :   long s, n, m, j, k, li, def, ldef;
     927                 :            : 
     928                 :       6916 :   RgM_dimensions(A, &m, &n);
     929         [ +  + ]:       6916 :   if (!n) return cgetg(1,t_MAT);
     930                 :       6846 :   av = avma;
     931                 :       6846 :   A = RgM_shallowcopy(A);
     932         [ +  + ]:       6846 :   def = n; ldef = (m>n)? m-n: 0;
     933         [ +  + ]:      25522 :   for (li=m; li>ldef; li--)
     934                 :            :   {
     935         [ +  + ]:      47831 :     for (j=def-1; j; j--)
     936                 :            :     {
     937                 :      29155 :       GEN a = gcoeff(A,li,j);
     938         [ +  + ]:      29155 :       if (!signe(a)) continue;
     939                 :            : 
     940                 :            :       /* zero a = Aij  using  b = Aik */
     941         [ +  + ]:      13510 :       k = (j==1)? def: j-1;
     942                 :      13510 :       ZC_elem(a,gcoeff(A,li,k), A,NULL, j,k);
     943         [ -  + ]:      13510 :       if (gc_needed(av,1))
     944                 :            :       {
     945         [ #  # ]:          0 :         if (DEBUGMEM>1) pari_warn(warnmem,"ZM_hnf[1]. li=%ld",li);
     946                 :          0 :         A = gerepilecopy(av, A);
     947                 :            :       }
     948                 :            :     }
     949                 :      18676 :     s = signe(gcoeff(A,li,def));
     950         [ +  + ]:      18676 :     if (s)
     951                 :            :     {
     952         [ +  + ]:      18648 :       if (s < 0) ZV_neg_inplace(gel(A,def));
     953                 :      18648 :       ZM_reduce(A, NULL, li,def);
     954                 :      18648 :       def--;
     955                 :            :     }
     956                 :            :     else
     957         [ +  + ]:         28 :       if (ldef) ldef--;
     958         [ -  + ]:      18676 :     if (gc_needed(av,1))
     959                 :            :     {
     960         [ #  # ]:          0 :       if (DEBUGMEM>1) pari_warn(warnmem,"ZM_hnf[2]. li=%ld",li);
     961                 :          0 :       A = gerepilecopy(av, A);
     962                 :            :     }
     963                 :            :   }
     964                 :            :   /* rank A = n - def */
     965         [ +  - ]:       6846 :   if (remove) { GEN B = NULL; remove_0cols(def, &A, &B, remove); }
     966                 :       6916 :   return gerepileupto(av0, ZM_copy(A));
     967                 :            : }
     968                 :            : 
     969                 :            : GEN
     970         [ +  + ]:       8524 : ZM_hnf(GEN x) { return lg(x) > 8? ZM_hnfall(x, NULL, 1): hnf_i(x, 1); }
     971                 :            : 
     972                 :            : /* u*z[1..k] mod p, in place */
     973                 :            : static void
     974                 :     917927 : FpV_Fp_mul_part_ip(GEN z, GEN u, GEN p, long k)
     975                 :            : {
     976                 :            :   long i;
     977         [ +  + ]:     917927 :   if (is_pm1(u)) {
     978         [ +  + ]:     476885 :     if (signe(u) > 0) {
     979         [ +  + ]:    1564851 :       for (i = 1; i <= k; i++)
     980         [ +  + ]:    1122758 :         if (signe(gel(z,i))) gel(z,i) = modii(gel(z,i), p);
     981                 :            :     } else {
     982         [ +  + ]:     101214 :       for (i = 1; i <= k; i++)
     983         [ +  + ]:      66422 :         if (signe(gel(z,i))) gel(z,i) = modii(negi(gel(z,i)), p);
     984                 :            :     }
     985                 :            :   }
     986                 :            :   else {
     987         [ +  + ]:    2540626 :     for (i = 1; i <= k; i++)
     988         [ +  + ]:    2099584 :       if (signe(gel(z,i))) gel(z,i) = Fp_mul(u,gel(z,i), p);
     989                 :            :   }
     990                 :     917927 : }
     991                 :            : static void
     992                 :   14896362 : FpV_red_part_ipvec(GEN z, GEN p, long k)
     993                 :            : {
     994                 :            :   long i;
     995         [ +  + ]:   41780984 :   for (i = 1; i <= k; i++) gel(z,i) = modii(gel(z,i), gel(p,i));
     996                 :   14896362 : }
     997                 :            : 
     998                 :            : /* return x * U, in echelon form (mod p^m), where (det(U),p) = 1.
     999                 :            :  * If early_abort is set, return NULL as soon as one pivot is 0 (mod p^m) */
    1000                 :            : GEN
    1001                 :      44565 : ZpM_echelon(GEN x, long early_abort, GEN p, GEN pm)
    1002                 :            : {
    1003                 :      44565 :   pari_sp av0 = avma, av;
    1004                 :            :   long m, li, co, i, j, k, def, ldef;
    1005                 :            : 
    1006         [ -  + ]:      44565 :   co = lg(x); if (co == 1) return cgetg(1,t_MAT);
    1007                 :      44565 :   li = lgcols(x);
    1008                 :      44565 :   av = avma;
    1009                 :      44565 :   x = RgM_shallowcopy(x);
    1010                 :      44565 :   m = Z_pval(pm, p);
    1011                 :            : 
    1012         [ -  + ]:      44565 :   ldef = (li > co)? li - co: 0;
    1013         [ +  + ]:     408564 :   for (def = co-1,i = li-1; i > ldef; i--)
    1014                 :            :   {
    1015                 :     364424 :     long vmin = LONG_MAX, kmin = 0;
    1016                 :     364424 :     GEN umin = gen_0, pvmin, q;
    1017         [ +  + ]:    3389419 :     for (k = 1; k <= def; k++)
    1018                 :            :     {
    1019                 :    3093657 :       GEN u = gcoeff(x,i,k);
    1020                 :            :       long v;
    1021         [ +  + ]:    3093657 :       if (!signe(u)) continue;
    1022                 :    1493667 :       v = Z_pvalrem(u, p, &u);
    1023         [ +  + ]:    1493667 :       if (v >= m) gcoeff(x,i,k) = gen_0;
    1024         [ +  + ]:     966481 :       else if (v < vmin) {
    1025                 :     356163 :         vmin = v; kmin = k; umin = u;
    1026         [ +  + ]:    1493667 :         if (!vmin) break;
    1027                 :            :       }
    1028                 :            :     }
    1029         [ +  + ]:     364424 :     if (!kmin)
    1030                 :            :     {
    1031         [ +  + ]:     149295 :       if (early_abort) return NULL;
    1032                 :     148870 :       gcoeff(x,i,def) = gen_0;
    1033                 :     148870 :       ldef--;
    1034         [ +  - ]:     148870 :       if (ldef < 0) ldef = 0;
    1035                 :     148870 :       continue;
    1036                 :            :     }
    1037         [ +  + ]:     215129 :     if (kmin != def) swap(gel(x,def), gel(x,kmin));
    1038         [ +  + ]:     215129 :     q = vmin? powiu(p, m-vmin): pm;
    1039                 :            :     /* pivot has valuation vmin */
    1040                 :     215129 :     umin = modii(umin, q);
    1041         [ +  + ]:     215129 :     if (!equali1(umin))
    1042                 :     171099 :       FpV_Fp_mul_part_ip(gel(x,def), Fp_inv(umin,q), pm, i-1);
    1043                 :     215129 :     gcoeff(x, i, def) = pvmin = powiu(p, vmin);
    1044         [ +  + ]:    2552684 :     for (j = def-1; j; j--)
    1045                 :            :     { /* zero x[i, 1..def-1] using x[i,def] = pvmin */
    1046                 :    2337555 :       GEN t, a = gcoeff(x,i,j) = modii(gcoeff(x,i,j), pm);
    1047         [ +  + ]:    2337555 :       if (!signe(a)) continue;
    1048                 :            : 
    1049                 :    1144800 :       t = diviiexact(a, pvmin); togglesign(t);
    1050                 :    1144800 :       ZC_lincomb1_inplace(gel(x,j), gel(x,def), t);
    1051         [ +  + ]:    1144800 :       if (gc_needed(av,1))
    1052                 :            :       {
    1053         [ -  + ]:         14 :         if (DEBUGMEM>1) pari_warn(warnmem,"ZpM_echelon. i=%ld",i);
    1054                 :         14 :         x = gerepilecopy(av, x); pvmin = gcoeff(x,i,def);
    1055                 :            :       }
    1056                 :            :     }
    1057                 :     215129 :     def--;
    1058                 :            :   }
    1059         [ -  + ]:      44140 :   if (co > li)
    1060                 :            :   {
    1061                 :          0 :     x += co - li;
    1062                 :          0 :     x[0] = evaltyp(t_MAT) | evallg(li);
    1063                 :            :   }
    1064                 :      44565 :   return gerepilecopy(av0, x);
    1065                 :            : }
    1066                 :            : GEN
    1067                 :     182331 : zlm_echelon(GEN x, long early_abort, ulong p, ulong pm)
    1068                 :            : {
    1069                 :     182331 :   pari_sp av0 = avma;
    1070                 :            :   long li, co, i, j, k, def, ldef;
    1071                 :            :   ulong m;
    1072                 :            : 
    1073         [ -  + ]:     182331 :   co = lg(x); if (co == 1) return cgetg(1,t_MAT);
    1074                 :     182331 :   li = lgcols(x);
    1075                 :     182331 :   x = Flm_copy(x);
    1076                 :     182331 :   m = u_lval(pm, p);
    1077                 :            : 
    1078         [ -  + ]:     182331 :   ldef = (li > co)? li - co: 0;
    1079         [ +  + ]:    1558905 :   for (def = co-1,i = li-1; i > ldef; i--)
    1080                 :            :   {
    1081                 :    1402658 :     long vmin = LONG_MAX, kmin = 0;
    1082                 :    1402658 :     ulong umin = 0, pvmin, q;
    1083         [ +  + ]:    7998980 :     for (k = 1; k <= def; k++)
    1084                 :            :     {
    1085                 :    7277513 :       ulong u = ucoeff(x,i,k);
    1086                 :            :       long v;
    1087         [ +  + ]:    7277513 :       if (!u) continue;
    1088                 :    3295701 :       v = u_lvalrem(u, p, &u);
    1089         [ -  + ]:    3295701 :       if (v >= (long) m) ucoeff(x,i,k) = 0;
    1090         [ +  + ]:    3295701 :       else if (v < vmin) {
    1091                 :    1432662 :         vmin = v; kmin = k; umin = u;
    1092         [ +  + ]:    3295701 :         if (!vmin) break;
    1093                 :            :       }
    1094                 :            :     }
    1095         [ +  + ]:    1402658 :     if (!kmin)
    1096                 :            :     {
    1097         [ +  + ]:     200866 :       if (early_abort) return NULL;
    1098                 :     174782 :       ucoeff(x,i,def) = 0;
    1099                 :     174782 :       ldef--;
    1100         [ +  - ]:     174782 :       if (ldef < 0) ldef = 0;
    1101                 :     174782 :       continue;
    1102                 :            :     }
    1103         [ +  + ]:    1201792 :     if (kmin != def) swap(gel(x,def), gel(x,kmin));
    1104         [ +  + ]:    1201792 :     q = vmin? upowuu(p, m-vmin): pm;
    1105                 :            :     /* pivot has valuation vmin */
    1106                 :    1201792 :     umin %= q;
    1107         [ +  + ]:    1201792 :     if (umin != 1)
    1108                 :     976514 :       Flc_Fl_mul_part_inplace(gel(x,def), Fl_inv(umin,q), pm, i-1);
    1109                 :    1201792 :     ucoeff(x, i, def) = pvmin = upowuu(p, vmin);
    1110         [ +  + ]:   13287640 :     for (j = def-1; j; j--)
    1111                 :            :     { /* zero x[i, 1..def-1] using x[i,def] = pvmin */
    1112                 :   12085848 :       ulong t, a = ucoeff(x,i,j);
    1113         [ +  + ]:   12085848 :       if (!a) continue;
    1114                 :            : 
    1115                 :    7130072 :       t = Fl_neg(a / pvmin, q);
    1116                 :    7130072 :       Flc_lincomb1_inplace(gel(x,j), gel(x,def), t, pm);
    1117                 :            :     }
    1118                 :    1201792 :     def--;
    1119                 :            :   }
    1120         [ -  + ]:     156247 :   if (co > li)
    1121                 :            :   {
    1122                 :          0 :     x += co - li;
    1123                 :          0 :     x[0] = evaltyp(t_MAT) | evallg(li);
    1124                 :            :   }
    1125                 :     182331 :   return gerepilecopy(av0, x);
    1126                 :            : }
    1127                 :            : 
    1128                 :            : /* dm = multiple of diag element (usually detint(x))
    1129                 :            :  * flag & hnf_MODID: reduce mod dm * matid [ otherwise as above ].
    1130                 :            :  * flag & hnf_PART: don't reduce once diagonal is known; */
    1131                 :            : 
    1132                 :            : /* x a ZM, dm a t_INT */
    1133                 :            : GEN
    1134                 :    2567014 : ZM_hnfmodall_i(GEN x, GEN dm, long flag)
    1135                 :            : {
    1136                 :            :   pari_sp av;
    1137                 :    2567014 :   const long center = (flag & hnf_CENTER);
    1138                 :    2567014 :   long moddiag = (flag & hnf_MODID);
    1139                 :            :   long li, co, i, j, k, def, ldef;
    1140                 :            :   GEN a, b, p1, p2, u, dm2, LDM;
    1141                 :            : 
    1142         [ -  + ]:    2567014 :   co = lg(x); if (co == 1) return cgetg(1,t_MAT);
    1143         [ +  + ]:    2567014 :   li = lgcols(x); if (li == 1) return cgetg(1,t_MAT);
    1144         [ +  + ]:    2566951 :   if (typ(dm) == t_INT)
    1145                 :            :   {
    1146                 :    2552384 :     long ldm = lgefint(dm);
    1147                 :    2552384 :     dm2 = shifti(dm, -1);
    1148                 :    2552384 :     dm = const_vec(li-1,dm);
    1149                 :    2552384 :     dm2= const_vec(li-1,dm2);
    1150                 :    2552384 :     LDM= const_vecsmall(li-1,ldm);
    1151                 :            :   }
    1152                 :            :   else
    1153                 :            :   {
    1154         [ -  + ]:      14567 :     if (lg(dm) != li) pari_err_DIM("ZM_hnfmod");
    1155                 :      14567 :     moddiag = 1;
    1156                 :      14567 :     dm2 = cgetg(li, t_VEC);
    1157                 :      14567 :     LDM = cgetg(li, t_VECSMALL);
    1158         [ +  + ]:      42497 :     for (i=1; i<li; i++)
    1159                 :            :     {
    1160                 :      27930 :       gel(dm2,i) = shifti(gel(dm,i),-1);
    1161                 :      27930 :       LDM[i] = lgefint(gel(dm,i));
    1162                 :            :     }
    1163                 :            :   }
    1164                 :    2566951 :   av = avma;
    1165                 :    2566951 :   x = RgM_shallowcopy(x);
    1166                 :            : 
    1167                 :    2566951 :   ldef = 0;
    1168         [ +  + ]:    2566951 :   if (li > co)
    1169                 :            :   {
    1170                 :       1652 :     ldef = li - co;
    1171         [ +  + ]:       1652 :     if (!moddiag)
    1172                 :          7 :       pari_err_DOMAIN("ZM_hnfmod","nb lines",">", strtoGENstr("nb columns"), x);
    1173                 :            :   }
    1174         [ +  + ]:   10727082 :   for (def = co-1,i = li-1; i > ldef; i--,def--)
    1175                 :            :   {
    1176                 :    8160138 :     GEN d = gel(dm,i), d2 = gel(dm2,i);
    1177                 :    8160138 :     gcoeff(x,i,def) = centermodii(gcoeff(x,i,def), d,d2);
    1178         [ +  + ]:   32591814 :     for (j = def-1; j; j--)
    1179                 :            :     {
    1180                 :   24431676 :       gcoeff(x,i,j) = centermodii(gcoeff(x,i,j), d,d2);
    1181                 :   24431676 :       a = gcoeff(x,i,j);
    1182         [ +  + ]:   24431676 :       if (!signe(a)) continue;
    1183                 :            : 
    1184         [ +  + ]:   12195368 :       k = (j==1)? def: j-1;
    1185                 :   12195368 :       gcoeff(x,i,k) = centermodii(gcoeff(x,i,k), d,d2);
    1186                 :   12195368 :       ZC_elem(a,gcoeff(x,i,k), x,NULL, j,k);
    1187                 :   12195368 :       p1 = gel(x,j);
    1188                 :   12195368 :       p2 = gel(x,k);
    1189                 :            :       /* prevent coeffs explosion: reduce mod dm when lg() > ldm */
    1190         [ +  + ]:   79918338 :       for (k = 1; k < i; k++)
    1191                 :            :       {
    1192         [ +  + ]:   67722970 :         if (lgefint(gel(p1,k)) > LDM[k])
    1193                 :    6195024 :           gel(p1,k) = centermodii(gel(p1,k), gel(dm,k),gel(dm2,k));
    1194         [ +  + ]:   67722970 :         if (lgefint(gel(p2,k)) > LDM[k])
    1195                 :    1081007 :           gel(p2,k) = centermodii(gel(p2,k), gel(dm,k),gel(dm2,k));
    1196                 :            :       }
    1197         [ +  + ]:   12195368 :       if (gc_needed(av,1))
    1198                 :            :       {
    1199         [ -  + ]:         48 :         if (DEBUGMEM>1) pari_warn(warnmem,"ZM_hnfmod[1]. i=%ld",i);
    1200                 :         48 :         x = gerepilecopy(av, x);
    1201                 :            :       }
    1202                 :            :     }
    1203 [ +  + ][ +  + ]:    8160138 :     if (moddiag && !signe(gcoeff(x,i,def)))
    1204                 :            :     { /* missing pivot on line i, insert column */
    1205                 :    3085740 :       GEN a = cgetg(co + 1, t_MAT);
    1206         [ +  + ]:   13866889 :       for (k = 1; k <= def; k++) gel(a,k) = gel(x,k);
    1207                 :    3085740 :       gel(a,k++) = Rg_col_ei(gel(dm,i), li-1, i);
    1208         [ +  + ]:   11287019 :       for (     ; k <= co;  k++) gel(a,k) = gel(x,k-1);
    1209         [ +  + ]:    3085740 :       ldef--; if (ldef < 0) ldef = 0;
    1210                 :    3085740 :       co++; def++; x = a;
    1211                 :            :     }
    1212                 :            :   }
    1213         [ +  + ]:    2566944 :   if (co < li)
    1214                 :            :   { /* implies moddiag, add missing diag(dm) components */
    1215                 :       1057 :     GEN a = cgetg(li+1, t_MAT);
    1216         [ +  + ]:       2261 :     for (k = 1; k <= li-co; k++) gel(a,k) = Rg_col_ei(gel(dm,k), li-1, k);
    1217         [ +  + ]:       2541 :     for (i = 1; i < co; i++) gel(a,k-1+i) = gel(x,i);
    1218                 :       1057 :     gel(a,li) = zerocol(li-1); x = a;
    1219                 :            :   }
    1220                 :            :   else
    1221                 :            :   {
    1222                 :    2565887 :     x += co - li;
    1223                 :    2565887 :     x[0] = evaltyp(t_MAT) | evallg(li); /* kill 0 columns */
    1224                 :            :   }
    1225         [ +  + ]:    2566944 :   if (moddiag)
    1226                 :            :   { /* one column extra: an accumulator, discarded at the end */
    1227         [ +  + ]:    2335497 :     if (lg(x) == li) x = shallowconcat(x, zerocol(li-1));
    1228                 :            :     /* add up missing diag(dm) components */
    1229         [ +  + ]:    9750011 :     for (i = li-1; i > 0; i--)
    1230                 :            :     {
    1231                 :    7414514 :       gcoeff(x, i, li) = gel(dm,i);
    1232         [ +  + ]:   28096367 :       for (j = i; j > 0; j--)
    1233                 :            :       {
    1234                 :   20681853 :         GEN a = gcoeff(x, j, li);
    1235         [ +  + ]:   20681853 :         if (!signe(a)) continue;
    1236                 :    7448181 :         ZC_elem(a, gcoeff(x,j,j), x, NULL, li,j);
    1237                 :    7448181 :         FpV_red_part_ipvec(gel(x,li), dm, j-1);
    1238                 :    7448181 :         FpV_red_part_ipvec(gel(x,j),  dm, j-1);
    1239         [ +  + ]:    7448181 :         if (gc_needed(av,1))
    1240                 :            :         {
    1241         [ -  + ]:          4 :           if (DEBUGMEM>1) pari_warn(warnmem,"ZM_hnfmod[2]. i=%ld", i);
    1242                 :          4 :           x = gerepilecopy(av, x);
    1243                 :            :         }
    1244                 :            :       }
    1245                 :            :     }
    1246                 :            :   }
    1247                 :            :   else
    1248                 :            :   {
    1249                 :     231447 :     b = gel(dm,1);
    1250         [ +  + ]:     978275 :     for (i = li-1; i > 0; i--)
    1251                 :            :     {
    1252                 :     746828 :       GEN d = bezout(gcoeff(x,i,i),b, &u,NULL);
    1253                 :     746828 :       gcoeff(x,i,i) = d;
    1254                 :     746828 :       FpV_Fp_mul_part_ip(gel(x,i), u, b, i-1);
    1255         [ +  + ]:     746828 :       if (i > 1) b = diviiexact(b,d);
    1256                 :            :     }
    1257                 :            :   }
    1258                 :    2566944 :   x[0] = evaltyp(t_MAT) | evallg(li); /* kill 0 columns / discard accumulator */
    1259         [ -  + ]:    2566944 :   if (flag & hnf_PART) return x;
    1260                 :            : 
    1261         [ +  + ]:    2566944 :   if (!moddiag)
    1262                 :            :   { /* compute optimal value for dm */
    1263                 :     231447 :     b = cgetg(li, t_VEC); gel(b,1) = gcoeff(x,1,1);
    1264         [ +  + ]:     746828 :     for (i = 2; i < li; i++) gel(b,i) = mulii(gel(b,i-1), gcoeff(x,i,i));
    1265                 :     231447 :     dm = b;
    1266                 :            :   }
    1267                 :            : 
    1268         [ +  + ]:   10728286 :   for (i = li-1; i > 0; i--)
    1269                 :            :   {
    1270                 :    8161342 :     GEN diag = gcoeff(x,i,i);
    1271         [ +  + ]:    8161342 :     if (signe(diag) < 0) { gel(x,i) = ZC_neg(gel(x,i)); diag = gcoeff(x,i,i); }
    1272         [ +  + ]:   22708735 :     for (j = i+1; j < li; j++)
    1273                 :            :     {
    1274                 :   14547393 :       b = gcoeff(x,i,j);
    1275         [ +  + ]:   14547393 :       b = center? diviiround(b,diag): truedivii(b, diag);
    1276         [ +  + ]:   14547393 :       if (!signe(b)) continue;
    1277                 :    5660652 :       togglesign(b);
    1278                 :    5660652 :       ZC_lincomb1_inplace(gel(x,j), gel(x,i),b);
    1279                 :    5660652 :       p1 = gel(x,j);
    1280         [ +  + ]:   23497596 :       for (k=1; k<i; k++)
    1281         [ +  + ]:   17836944 :         if (lgefint(gel(p1,k)) > LDM[k]) gel(p1,k) = remii(gel(p1,k), gel(dm,i));
    1282         [ +  + ]:    5660652 :       if (gc_needed(av,1))
    1283                 :            :       {
    1284         [ -  + ]:         29 :         if (DEBUGMEM>1) pari_warn(warnmem,"ZM_hnfmod[3]. i=%ld", i);
    1285                 :         29 :         gerepileall(av, 2, &x, &dm); diag = gcoeff(x,i,i);
    1286                 :            :       }
    1287                 :            :     }
    1288                 :            :   }
    1289                 :    2567007 :   return x;
    1290                 :            : }
    1291                 :            : GEN
    1292                 :    2563563 : ZM_hnfmodall(GEN x, GEN dm, long flag)
    1293                 :            : {
    1294                 :    2563563 :   pari_sp av = avma;
    1295                 :    2563563 :   return gerepilecopy(av, ZM_hnfmodall_i(x, dm, flag));
    1296                 :            : }
    1297                 :            : GEN
    1298                 :     231440 : ZM_hnfmod(GEN x, GEN d) { return ZM_hnfmodall(x,d,0); }
    1299                 :            : GEN
    1300                 :    2331885 : ZM_hnfmodid(GEN x, GEN d) { return ZM_hnfmodall(x,d,hnf_MODID); }
    1301                 :            : 
    1302                 :            : static GEN
    1303                 :         84 : allhnfmod(GEN x, GEN dm, int flag)
    1304                 :            : {
    1305         [ -  + ]:         84 :   if (typ(x)!=t_MAT) pari_err_TYPE("allhnfmod",x);
    1306                 :         84 :   RgM_check_ZM(x, "allhnfmod");
    1307         [ -  + ]:         84 :   if (isintzero(dm)) return ZM_hnf(x);
    1308                 :         84 :   return ZM_hnfmodall(x, dm, flag);
    1309                 :            : }
    1310                 :            : GEN
    1311                 :         14 : hnfmod(GEN x, GEN d)
    1312                 :            : {
    1313         [ -  + ]:         14 :   if (typ(d) != t_INT) pari_err_TYPE("mathnfmod",d);
    1314                 :         14 :   return allhnfmod(x, d, 0);
    1315                 :            : }
    1316                 :            : GEN
    1317                 :         70 : hnfmodid(GEN x, GEN d)
    1318                 :            : {
    1319      [ +  +  - ]:         70 :   switch(typ(d))
    1320                 :            :   {
    1321                 :         56 :     case t_INT: break;
    1322                 :            :     case t_VEC: case t_COL:
    1323         [ +  - ]:         14 :       if (RgV_is_ZV(d)) break;
    1324                 :          0 :     default: pari_err_TYPE("mathnfmodid",d);
    1325                 :            :   }
    1326                 :         70 :   return allhnfmod(x, d, hnf_MODID);
    1327                 :            : }
    1328                 :            : 
    1329                 :            : /* M a ZM in HNF. Normalize with *centered* residues */
    1330                 :            : GEN
    1331                 :        616 : ZM_hnfcenter(GEN M)
    1332                 :            : {
    1333                 :        616 :   long i, j, k, N = lg(M)-1;
    1334                 :        616 :   pari_sp av = avma;
    1335                 :            : 
    1336         [ +  + ]:       2520 :   for (j=N-1; j>0; j--) /* skip last line */
    1337                 :            :   {
    1338                 :       1904 :     GEN Mj = gel(M,j), a = gel(Mj,j);
    1339         [ +  + ]:      12691 :     for (k = j+1; k <= N; k++)
    1340                 :            :     {
    1341                 :      10787 :       GEN Mk = gel(M,k), q = diviiround(gel(Mk,j), a);
    1342                 :      10787 :       long s = signe(q);
    1343         [ +  + ]:      10787 :       if (!s) continue;
    1344         [ +  + ]:       9079 :       if (is_pm1(q))
    1345                 :            :       {
    1346         [ +  + ]:       1120 :         if (s < 0)
    1347         [ +  + ]:        630 :           for (i = 1; i <= j; i++) gel(Mk,i) = addii(gel(Mk,i), gel(Mj,i));
    1348                 :            :         else
    1349         [ +  + ]:       3311 :           for (i = 1; i <= j; i++) gel(Mk,i) = subii(gel(Mk,i), gel(Mj,i));
    1350                 :            :       }
    1351                 :            :       else
    1352         [ +  + ]:      88865 :         for (i = 1; i <= j; i++) gel(Mk,i) = subii(gel(Mk,i), mulii(q,gel(Mj,i)));
    1353         [ -  + ]:       9079 :       if (gc_needed(av,1))
    1354                 :            :       {
    1355         [ #  # ]:          0 :         if (DEBUGMEM) pari_warn(warnmem,"ZM_hnfcenter, j = %ld",j);
    1356                 :          0 :         M = gerepilecopy(av, M);
    1357                 :            :       }
    1358                 :            :     }
    1359                 :            :   }
    1360                 :        616 :   return M;
    1361                 :            : }
    1362                 :            : 
    1363                 :            : /***********************************************************************/
    1364                 :            : /*                                                                     */
    1365                 :            : /*                 HNFLLL (Havas, Majewski, Mathews)                   */
    1366                 :            : /*                                                                     */
    1367                 :            : /***********************************************************************/
    1368                 :            : 
    1369                 :            : static void
    1370                 :     290843 : Minus(long j, GEN lambda)
    1371                 :            : {
    1372                 :     290843 :   long k, n = lg(lambda);
    1373                 :            : 
    1374         [ +  + ]:    2392282 :   for (k=1  ; k<j; k++) togglesign_safe(&gcoeff(lambda,k,j));
    1375         [ +  + ]:    3325858 :   for (k=j+1; k<n; k++) togglesign_safe(&gcoeff(lambda,j,k));
    1376                 :     290843 : }
    1377                 :            : 
    1378                 :            : /* index of first non-zero entry */
    1379                 :            : static long
    1380                 :   23523209 : findi(GEN M)
    1381                 :            : {
    1382                 :   23523209 :   long i, n = lg(M);
    1383         [ +  + ]:  174784236 :   for (i=1; i<n; i++)
    1384         [ +  + ]:  162629256 :     if (signe(gel(M,i))) return i;
    1385                 :   23523209 :   return 0;
    1386                 :            : }
    1387                 :            : 
    1388                 :            : static long
    1389                 :   23523209 : findi_normalize(GEN Aj, GEN B, long j, GEN lambda)
    1390                 :            : {
    1391                 :   23523209 :   long r = findi(Aj);
    1392 [ +  + ][ +  + ]:   23523209 :   if (r && signe(gel(Aj,r)) < 0)
    1393                 :            :   {
    1394         [ +  - ]:     290843 :     ZV_togglesign(Aj); if (B) ZV_togglesign(gel(B,j));
    1395                 :     290843 :     Minus(j,lambda);
    1396                 :            :   }
    1397                 :   23523209 :   return r;
    1398                 :            : }
    1399                 :            : 
    1400                 :            : static void
    1401                 :   11761496 : reduce2(GEN A, GEN B, long k, long j, long *row0, long *row1, GEN lambda, GEN D)
    1402                 :            : {
    1403                 :            :   GEN q;
    1404                 :            :   long i;
    1405                 :            : 
    1406                 :   11761496 :   *row0 = findi_normalize(gel(A,j), B,j,lambda);
    1407                 :   11761496 :   *row1 = findi_normalize(gel(A,k), B,k,lambda);
    1408         [ +  + ]:   11761496 :   if (*row0)
    1409                 :    4160771 :     q = truedivii(gcoeff(A,*row0,k), gcoeff(A,*row0,j));
    1410         [ +  + ]:    7600725 :   else if (absi_cmp(shifti(gcoeff(lambda,j,k), 1), gel(D,j)) > 0)
    1411                 :    1972732 :     q = diviiround(gcoeff(lambda,j,k), gel(D,j));
    1412                 :            :   else
    1413                 :   11761496 :     return;
    1414                 :            : 
    1415         [ +  + ]:    6133503 :   if (signe(q))
    1416                 :            :   {
    1417                 :    3563362 :     GEN Lk = gel(lambda,k), Lj = gel(lambda,j);
    1418                 :    3563362 :     togglesign_safe(&q);
    1419         [ +  + ]:    3563362 :     if (*row0) ZC_lincomb1_inplace(gel(A,k),gel(A,j),q);
    1420         [ +  - ]:    3563362 :     if (B) ZC_lincomb1_inplace(gel(B,k),gel(B,j),q);
    1421                 :    3563362 :     gel(Lk,j) = addii(gel(Lk,j), mulii(q,gel(D,j)));
    1422         [ +  + ]:    3563362 :     if (is_pm1(q))
    1423                 :            :     {
    1424         [ +  + ]:     941185 :       if (signe(q) > 0)
    1425                 :            :       {
    1426         [ +  + ]:    3241172 :         for (i=1; i<j; i++)
    1427         [ +  + ]:    2900676 :           if (signe(gel(Lj,i))) gel(Lk,i) = addii(gel(Lk,i), gel(Lj,i));
    1428                 :            :       }
    1429                 :            :       else
    1430                 :            :       {
    1431         [ +  + ]:    4151863 :         for (i=1; i<j; i++)
    1432         [ +  + ]:    3551174 :           if (signe(gel(Lj,i))) gel(Lk,i) = subii(gel(Lk,i), gel(Lj,i));
    1433                 :            :       }
    1434                 :            :     }
    1435                 :            :     else
    1436                 :            :     {
    1437         [ +  + ]:   24045450 :       for (i=1; i<j; i++)
    1438         [ +  + ]:   12283954 :         if (signe(gel(Lj,i))) gel(Lk,i) = addii(gel(Lk,i), mulii(q,gel(Lj,i)));
    1439                 :            :     }
    1440                 :            :   }
    1441                 :            : }
    1442                 :            : 
    1443                 :            : static void
    1444                 :    1655329 : hnfswap(GEN A, GEN B, long k, GEN lambda, GEN D)
    1445                 :            : {
    1446                 :    1655329 :   GEN t, p1, p2, Lk = gel(lambda,k);
    1447                 :    1655329 :   long i,j,n = lg(A);
    1448                 :            : 
    1449                 :    1655329 :   swap(gel(A,k), gel(A,k-1));
    1450         [ +  - ]:    1655329 :   if (B) swap(gel(B,k), gel(B,k-1));
    1451         [ +  + ]:   10509815 :   for (j=k-2; j; j--) swap(gcoeff(lambda,j,k-1), gel(Lk,j));
    1452         [ +  + ]:   24370622 :   for (i=k+1; i<n; i++)
    1453                 :            :   {
    1454                 :   22715293 :     GEN Li = gel(lambda,i);
    1455                 :   22715293 :     p1 = mulii(gel(Li,k-1), gel(D,k));
    1456                 :   22715293 :     p2 = mulii(gel(Li,k), gel(Lk,k-1));
    1457                 :   22715293 :     t = subii(p1,p2);
    1458                 :            : 
    1459                 :   22715293 :     p1 = mulii(gel(Li,k), gel(D,k-2));
    1460                 :   22715293 :     p2 = mulii(gel(Li,k-1), gel(Lk,k-1));
    1461                 :   22715293 :     gel(Li,k-1) = diviiexact(addii(p1,p2), gel(D,k-1));
    1462                 :   22715293 :     gel(Li,k) = diviiexact(t, gel(D,k-1));
    1463                 :            :   }
    1464                 :    1655329 :   p1 = mulii(gel(D,k-2), gel(D,k));
    1465                 :    1655329 :   p2 = sqri(gel(Lk,k-1));
    1466                 :    1655329 :   gel(D,k-1) = diviiexact(addii(p1,p2), gel(D,k-1));
    1467                 :    1655329 : }
    1468                 :            : 
    1469                 :            : /* reverse row order in matrix A, IN PLACE */
    1470                 :            : static GEN
    1471                 :      33634 : reverse_rows(GEN A)
    1472                 :            : {
    1473                 :      33634 :   long i, j, h, n = lg(A);
    1474         [ +  + ]:      33634 :   if (n == 1) return A;
    1475                 :      33620 :   h = lgcols(A);
    1476         [ +  + ]:     374180 :   for (j=1; j<n; j++)
    1477                 :            :   {
    1478                 :     340560 :     GEN c = gel(A,j);
    1479                 :            :     /* start at (h-1) >>1 : if h = 2i even, no need to swap c[i] and itself */
    1480         [ +  + ]:    1423540 :     for (i=(h-1)>>1; i; i--) swap(gel(c,i), gel(c,h-i));
    1481                 :            :   }
    1482                 :      33634 :   return A;
    1483                 :            : }
    1484                 :            : 
    1485                 :            : GEN
    1486                 :      16817 : ZM_hnflll(GEN A, GEN *ptB, int remove)
    1487                 :            : {
    1488                 :      16817 :   pari_sp av = avma;
    1489                 :            : #ifdef HNFLLL_QUALITY
    1490                 :            :   const long m1 = 1, n1 = 1; /* alpha = m1/n1. Maybe 3/4 here ? */
    1491                 :            : #endif
    1492                 :            :   long n, k, kmax;
    1493                 :            :   GEN B, lambda, D;
    1494                 :            : 
    1495                 :      16817 :   n = lg(A);
    1496                 :      16817 :   A = reverse_rows(ZM_copy(A)); /* ZM_copy for in place findi_normalize() */
    1497         [ +  - ]:      16817 :   B = ptB? matid(n-1): NULL;
    1498                 :      16817 :   D = const_vec(n, gen_1) + 1;
    1499                 :      16817 :   lambda = zeromatcopy(n-1,n-1);
    1500                 :      16817 :   k = kmax = 2;
    1501         [ +  + ]:    2796427 :   while (k < n)
    1502                 :            :   {
    1503                 :            :     long row0, row1;
    1504                 :            :     int do_swap;
    1505                 :    2779610 :     reduce2(A,B,k,k-1,&row0,&row1,lambda,D);
    1506         [ +  + ]:    2779610 :     if (row0)
    1507 [ +  + ][ +  + ]:    1859300 :       do_swap = (!row1 || row0 <= row1);
    1508         [ +  + ]:     920310 :     else if (!row1)
    1509                 :            :     { /* row0 == row1 == 0 */
    1510                 :     838578 :       pari_sp av1 = avma;
    1511                 :     838578 :       GEN z = addii(mulii(gel(D,k-2),gel(D,k)), sqri(gcoeff(lambda,k-1,k)));
    1512                 :            : #ifdef HNFLLL_QUALITY
    1513                 :            :       do_swap = (cmpii(mului(n1,z), mului(m1,sqri(gel(D,k-1)))) < 0);
    1514                 :            : #else /* assume m1 = n1 = 1 */
    1515                 :     838578 :       do_swap = (cmpii(z, sqri(gel(D,k-1))) < 0);
    1516                 :            : #endif
    1517                 :     838578 :       avma = av1;
    1518                 :            :     }
    1519                 :            :     else
    1520                 :      81732 :       do_swap = 0;
    1521         [ +  + ]:    2779610 :     if (do_swap)
    1522                 :            :     {
    1523                 :    1639973 :       hnfswap(A,B,k,lambda,D);
    1524         [ +  + ]:    1639973 :       if (k > 2) k--;
    1525                 :            :     }
    1526                 :            :     else
    1527                 :            :     {
    1528                 :            :       long i;
    1529         [ +  + ]:   10121523 :       for (i=k-2; i; i--)
    1530                 :            :       {
    1531                 :            :         long row0, row1;
    1532                 :    8981886 :         reduce2(A,B,k,i,&row0,&row1,lambda,D);
    1533         [ +  + ]:    8981886 :         if (gc_needed(av,3))
    1534                 :            :         {
    1535                 :         39 :           GEN b = D-1;
    1536         [ -  + ]:         39 :           if (DEBUGMEM) pari_warn(warnmem,"hnflll (reducing), kmax = %ld",kmax);
    1537         [ +  - ]:         39 :           gerepileall(av, B? 4: 3, &A, &lambda, &b, &B);
    1538                 :         39 :           D = b+1;
    1539                 :            :         }
    1540                 :            :       }
    1541         [ +  + ]:    1139637 :       if (++k > kmax) kmax = k;
    1542                 :            :     }
    1543         [ +  + ]:    2779610 :     if (gc_needed(av,3))
    1544                 :            :     {
    1545                 :         41 :       GEN b = D-1;
    1546         [ -  + ]:         41 :       if (DEBUGMEM) pari_warn(warnmem,"hnflll, kmax = %ld / %ld",kmax,n-1);
    1547         [ +  - ]:         41 :       gerepileall(av, B? 4: 3, &A, &lambda, &b, &B);
    1548                 :    2779610 :       D = b+1;
    1549                 :            :     }
    1550                 :            :   }
    1551                 :            :   /* handle trivial case: return negative diag coefficient otherwise */
    1552         [ +  + ]:      16817 :   if (n == 2) (void)findi_normalize(gel(A,1), B,1,lambda);
    1553                 :      16817 :   A = reverse_rows(A);
    1554         [ +  + ]:      16817 :   if (remove)
    1555                 :            :   {
    1556                 :            :     long i;
    1557         [ +  + ]:      18065 :     for (i = 1; i < n; i++)
    1558         [ +  + ]:      17834 :       if (!ZV_equal0(gel(A,i))) break;
    1559                 :       1434 :     remove_0cols(i-1, &A, &B, remove);
    1560                 :            :   }
    1561         [ +  - ]:      16817 :   gerepileall(av, B? 2: 1, &A, &B);
    1562         [ +  - ]:      16817 :   if (B) *ptB = B;
    1563                 :      16817 :   return A;
    1564                 :            : }
    1565                 :            : 
    1566                 :            : GEN
    1567                 :          7 : hnflll(GEN x)
    1568                 :            : {
    1569                 :          7 :   GEN z = cgetg(3, t_VEC);
    1570                 :          7 :   gel(z,1) = ZM_hnflll(x, &gel(z,2), 1);
    1571                 :          7 :   return z;
    1572                 :            : }
    1573                 :            : 
    1574                 :            : /* Variation on HNFLLL: Extended GCD */
    1575                 :            : 
    1576                 :            : static void
    1577                 :     248353 : reduce1(GEN A, GEN B, long k, long j, GEN lambda, GEN D)
    1578                 :            : {
    1579                 :            :   GEN q;
    1580                 :            :   long i;
    1581                 :            : 
    1582         [ +  + ]:     248353 :   if (signe(gel(A,j)))
    1583                 :      12028 :     q = diviiround(gel(A,k),gel(A,j));
    1584         [ +  + ]:     236325 :   else if (absi_cmp(shifti(gcoeff(lambda,j,k), 1), gel(D,j)) > 0)
    1585                 :       1598 :     q = diviiround(gcoeff(lambda,j,k), gel(D,j));
    1586                 :            :   else
    1587                 :     248353 :     return;
    1588                 :            : 
    1589         [ +  + ]:      13626 :   if (signe(q))
    1590                 :            :   {
    1591                 :       3947 :     GEN Lk = gel(lambda,k), Lj = gel(lambda,j);
    1592                 :       3947 :     togglesign_safe(&q);
    1593                 :       3947 :     gel(A,k) = addii(gel(A,k), mulii(q,gel(A,j)));
    1594                 :       3947 :     ZC_lincomb1_inplace(gel(B,k),gel(B,j),q);
    1595                 :       3947 :     gel(Lk,j) = addii(gel(Lk,j), mulii(q,gel(D,j)));
    1596         [ +  + ]:     261106 :     for (i=1; i<j; i++)
    1597         [ +  + ]:      12753 :       if (signe(gel(Lj,i))) gel(Lk,i) = addii(gel(Lk,i), mulii(q,gel(Lj,i)));
    1598                 :            :   }
    1599                 :            : }
    1600                 :            : 
    1601                 :            : static GEN
    1602                 :      10969 : ZV_gcdext_i(GEN A)
    1603                 :            : {
    1604                 :            : #ifdef HNFLLL_QUALITY
    1605                 :            :   const long m1 = 1, n1 = 1; /* alpha = m1/n1. Maybe 3/4 here ? */
    1606                 :            : #endif
    1607                 :      10969 :   long k, n = lg(A);
    1608                 :            :   GEN B, lambda, D;
    1609                 :            : 
    1610         [ -  + ]:      10969 :   if (n == 1) retmkvec2(gen_1, cgetg(1,t_MAT));
    1611                 :      10969 :   A = leafcopy(A);
    1612                 :      10969 :   B = matid(n-1);
    1613                 :      10969 :   lambda = zeromatcopy(n-1,n-1);
    1614                 :      10969 :   D = const_vec(n, gen_1) + 1;
    1615                 :      10969 :   k = 2;
    1616         [ +  + ]:      93390 :   while (k < n)
    1617                 :            :   {
    1618                 :            :     int do_swap;
    1619                 :            : 
    1620                 :      82421 :     reduce1(A,B,k,k-1,lambda,D);
    1621         [ +  + ]:      82421 :     if (signe(gel(A,k-1))) do_swap = 1;
    1622         [ +  + ]:      70393 :     else if (!signe(gel(A,k)))
    1623                 :            :     {
    1624                 :      47487 :       pari_sp av1 = avma;
    1625                 :      47487 :       GEN z = addii(mulii(gel(D,k-2),gel(D,k)), sqri(gcoeff(lambda,k-1,k)));
    1626                 :            : #ifdef HNFLLL_QUALITY
    1627                 :            :       do_swap = (cmpii(mului(n1,z), mului(m1,sqri(gel(D,k-1)))) < 0);
    1628                 :            : #else /* assume m1 = n1 = 1 */
    1629                 :      47487 :       do_swap = (cmpii(z, sqri(gel(D,k-1))) < 0);
    1630                 :            : #endif
    1631                 :      47487 :       avma = av1;
    1632                 :            :     }
    1633                 :      22906 :     else do_swap = 0;
    1634                 :            : 
    1635         [ +  + ]:      82421 :     if (do_swap)
    1636                 :            :     {
    1637                 :      15356 :       hnfswap(A,B,k,lambda,D);
    1638         [ +  + ]:      15356 :       if (k > 2) k--;
    1639                 :            :     }
    1640                 :            :     else
    1641                 :            :     {
    1642                 :            :       long i;
    1643         [ +  + ]:     232997 :       for (i=k-2; i; i--) reduce1(A,B,k,i,lambda,D);
    1644                 :      67065 :       k++;
    1645                 :            :     }
    1646                 :            :   }
    1647         [ +  + ]:      10969 :   if (signe(gel(A,n-1)) < 0)
    1648                 :            :   {
    1649                 :       5262 :     gel(A,n-1) = negi(gel(A,n-1));
    1650                 :       5262 :     ZV_togglesign(gel(B,n-1));
    1651                 :            :   }
    1652                 :      10969 :   return mkvec2(gel(A,n-1), B);
    1653                 :            : }
    1654                 :            : GEN
    1655                 :          0 : ZV_extgcd(GEN A)
    1656                 :            : {
    1657                 :          0 :   pari_sp av = avma;
    1658                 :          0 :   return gerepilecopy(av, ZV_gcdext_i(A));
    1659                 :            : }
    1660                 :            : /* as ZV_extgcd, transforming the gcd into a t_MAT, for mathnf0 */
    1661                 :            : static GEN
    1662                 :      10976 : ZV_hnfgcdext(GEN A)
    1663                 :            : {
    1664                 :      10976 :   pari_sp av = avma;
    1665                 :            :   GEN z;
    1666         [ +  + ]:      10976 :   if (lg(A) == 1) retmkvec2(cgetg(1,t_MAT),cgetg(1,t_MAT));
    1667                 :      10969 :   z = ZV_gcdext_i(A);
    1668                 :      10969 :   gel(z,1) = mkmat(mkcol(gel(z,1)));
    1669                 :      10976 :   return gerepilecopy(av, z);
    1670                 :            : }
    1671                 :            : 
    1672                 :            : /* HNF with permutation. */
    1673                 :            : GEN
    1674                 :        126 : ZM_hnfperm(GEN A, GEN *ptU, GEN *ptperm)
    1675                 :            : {
    1676                 :            :   GEN U, c, l, perm, d, p, q, b;
    1677                 :        126 :   pari_sp av = avma, av1;
    1678                 :            :   long r, t, i, j, j1, k, m, n;
    1679                 :            : 
    1680                 :        126 :   n = lg(A)-1;
    1681         [ +  + ]:        126 :   if (!n)
    1682                 :            :   {
    1683         [ +  - ]:          7 :     if (ptU) *ptU = cgetg(1,t_MAT);
    1684         [ +  - ]:          7 :     if (ptperm) *ptperm = cgetg(1,t_VEC);
    1685                 :          7 :     return cgetg(1, t_MAT);
    1686                 :            :   }
    1687                 :        119 :   m = nbrows(A);
    1688                 :        119 :   c = zero_zv(m);
    1689                 :        119 :   l = zero_zv(n);
    1690                 :        119 :   perm = cgetg(m+1, t_VECSMALL);
    1691                 :        119 :   av1 = avma;
    1692                 :        119 :   A = RgM_shallowcopy(A);
    1693         [ +  + ]:        119 :   U = ptU? matid(n): NULL;
    1694                 :            :   /* U base change matrix : A0*U = A all along */
    1695         [ +  + ]:        560 :   for (r=0, k=1; k <= n; k++)
    1696                 :            :   {
    1697         [ +  + ]:       1435 :     for (j=1; j<k; j++)
    1698                 :            :     {
    1699         [ +  + ]:        994 :       if (!l[j]) continue;
    1700                 :        882 :       t=l[j]; b=gcoeff(A,t,k);
    1701         [ +  + ]:        882 :       if (!signe(b)) continue;
    1702                 :            : 
    1703                 :        364 :       ZC_elem(b,gcoeff(A,t,j), A,U,k,j);
    1704                 :        364 :       d = gcoeff(A,t,j);
    1705         [ -  + ]:        364 :       if (signe(d) < 0)
    1706                 :            :       {
    1707                 :          0 :         ZV_neg_inplace(gel(A,j));
    1708         [ #  # ]:          0 :         if (U) ZV_togglesign(gel(U,j));
    1709                 :          0 :         d = gcoeff(A,t,j);
    1710                 :            :       }
    1711         [ +  + ]:        966 :       for (j1=1; j1<j; j1++)
    1712                 :            :       {
    1713         [ +  + ]:        602 :         if (!l[j1]) continue;
    1714                 :        595 :         q = truedivii(gcoeff(A,t,j1),d);
    1715         [ +  + ]:        595 :         if (!signe(q)) continue;
    1716                 :            : 
    1717                 :        329 :         togglesign(q);
    1718                 :        329 :         ZC_lincomb1_inplace(gel(A,j1), gel(A,j), q);
    1719         [ +  - ]:        329 :         if (U) ZC_lincomb1_inplace(gel(U,j1), gel(U,j), q);
    1720                 :            :       }
    1721                 :            :     }
    1722 [ +  + ][ +  + ]:       1288 :     t = m; while (t && (c[t] || !signe(gcoeff(A,t,k)))) t--;
                 [ +  + ]
    1723         [ +  + ]:        441 :     if (t)
    1724                 :            :     {
    1725                 :        357 :       p = gcoeff(A,t,k);
    1726         [ +  + ]:       6230 :       for (i=t-1; i; i--)
    1727                 :            :       {
    1728                 :       5873 :         q = gcoeff(A,i,k);
    1729 [ +  + ][ +  + ]:       5873 :         if (signe(q) && absi_cmp(p,q) > 0) { p = q; t = i; }
    1730                 :            :       }
    1731                 :        357 :       perm[++r] = l[k] = t; c[t] = k;
    1732         [ +  + ]:        357 :       if (signe(p) < 0)
    1733                 :            :       {
    1734                 :         68 :         ZV_neg_inplace(gel(A,k));
    1735         [ -  + ]:         68 :         if (U) ZV_togglesign(gel(U,k));
    1736                 :         68 :         p = gcoeff(A,t,k);
    1737                 :            :       }
    1738                 :            :       /* p > 0 */
    1739         [ +  + ]:       1029 :       for (j=1; j<k; j++)
    1740                 :            :       {
    1741         [ +  + ]:        672 :         if (!l[j]) continue;
    1742                 :        665 :         q = truedivii(gcoeff(A,t,j),p);
    1743         [ +  + ]:        665 :         if (!signe(q)) continue;
    1744                 :            : 
    1745                 :        167 :         togglesign(q);
    1746                 :        167 :         ZC_lincomb1_inplace(gel(A,j), gel(A,k), q);
    1747         [ +  + ]:        167 :         if (U) ZC_lincomb1_inplace(gel(U,j), gel(U,k), q);
    1748                 :            :       }
    1749                 :            :     }
    1750         [ -  + ]:        441 :     if (gc_needed(av1,1))
    1751                 :            :     {
    1752         [ #  # ]:          0 :       if (DEBUGMEM>1) pari_warn(warnmem,"hnfperm");
    1753         [ #  # ]:          0 :       gerepileall(av1, U? 2: 1, &A, &U);
    1754                 :            :     }
    1755                 :            :   }
    1756         [ +  + ]:        119 :   if (r < m)
    1757                 :            :   {
    1758         [ +  + ]:       1351 :     for (i=1,k=r; i<=m; i++)
    1759         [ +  + ]:       1281 :       if (!c[i]) perm[++k] = i;
    1760                 :            :   }
    1761                 :            : 
    1762                 :            : /* We have A0*U=A, U in Gl(n,Z)
    1763                 :            :  * basis for Im(A):  columns of A s.t l[j]>0 (r   cols)
    1764                 :            :  * basis for Ker(A): columns of U s.t l[j]=0 (n-r cols) */
    1765                 :        119 :   p = cgetg(r+1,t_MAT);
    1766         [ +  + ]:        784 :   for (i=1; i<=m/2; i++) lswap(perm[i], perm[m+1-i]);
    1767         [ +  + ]:        119 :   if (U)
    1768                 :            :   {
    1769                 :         56 :     GEN u = cgetg(n+1,t_MAT);
    1770         [ +  + ]:        266 :     for (t=1,k=r,j=1; j<=n; j++)
    1771         [ +  + ]:        210 :       if (l[j])
    1772                 :            :       {
    1773                 :        126 :         u[k + n-r] = U[j];
    1774                 :        126 :         gel(p,k--) = vecpermute(gel(A,j), perm);
    1775                 :            :       }
    1776                 :            :       else
    1777                 :         84 :         u[t++] = U[j];
    1778                 :         56 :     *ptU = u;
    1779         [ +  - ]:         56 :     if (ptperm) *ptperm = perm;
    1780         [ +  - ]:         56 :     gerepileall(av, ptperm? 3: 2, &p, ptU, ptperm);
    1781                 :            :   }
    1782                 :            :   else
    1783                 :            :   {
    1784         [ +  + ]:        294 :     for (k=r,j=1; j<=n; j++)
    1785         [ +  - ]:        231 :       if (l[j]) gel(p,k--) = vecpermute(gel(A,j), perm);
    1786         [ +  + ]:         63 :     if (ptperm) *ptperm = perm;
    1787         [ +  + ]:         63 :     gerepileall(av, ptperm? 2: 1, &p, ptperm);
    1788                 :            :   }
    1789                 :        126 :   return p;
    1790                 :            : }
    1791                 :            : 
    1792                 :            : GEN
    1793                 :         14 : hnfperm(GEN A)
    1794                 :            : {
    1795                 :         14 :   GEN y = cgetg(4, t_VEC);
    1796                 :         14 :   gel(y,1) = ZM_hnfperm(A, &gel(y,2), &gel(y,3));
    1797                 :         14 :   return y;
    1798                 :            : }
    1799                 :            : 
    1800                 :            : /* Hermite Normal Form, with base change matrix if ptB != NULL.
    1801                 :            :  * If 'remove' = 1, remove 0 columns (do NOT update *ptB accordingly)
    1802                 :            :  * If 'remove' = 2, remove 0 columns and update *ptB accordingly */
    1803                 :            : GEN
    1804                 :     103725 : ZM_hnfall(GEN A, GEN *ptB, long remove)
    1805                 :            : {
    1806                 :     103725 :   pari_sp av = avma, av1;
    1807                 :            :   long m, n, r, i, j, k, li;
    1808                 :            :   GEN B, c, h, a;
    1809                 :            : 
    1810                 :     103725 :   RgM_dimensions(A, &m,&n);
    1811         [ +  + ]:     103725 :   if (!n)
    1812                 :            :   {
    1813         [ +  - ]:          7 :     if (ptB) *ptB = cgetg(1,t_MAT);
    1814                 :          7 :     return cgetg(1,t_MAT);
    1815                 :            :   }
    1816                 :     103718 :   c = zero_zv(m);
    1817                 :     103718 :   h = const_vecsmall(n, m);
    1818                 :     103718 :   av1 = avma;
    1819                 :     103718 :   A = RgM_shallowcopy(A);
    1820         [ +  + ]:     103718 :   B = ptB? matid(n): NULL;
    1821                 :     103718 :   r = n+1;
    1822         [ +  + ]:     365718 :   for (li=m; li; li--)
    1823                 :            :   {
    1824         [ +  - ]:     835642 :     for (j=1; j<r; j++)
    1825                 :            :     {
    1826         [ +  + ]:    1087906 :       for (i=h[j]; i>li; i--)
    1827                 :            :       {
    1828                 :     252264 :         a = gcoeff(A,i,j);
    1829                 :     252264 :         k = c[i];
    1830                 :            :         /* zero a = Aij  using  Aik */
    1831         [ +  + ]:     252264 :         if (signe(a)) ZC_elem(a,gcoeff(A,i,k), A,B,j,k);
    1832                 :     252264 :         ZM_reduce(A,B, i,k); /* ensure reduced entries */
    1833         [ -  + ]:     252264 :         if (gc_needed(av1,1))
    1834                 :            :         {
    1835         [ #  # ]:          0 :           if (DEBUGMEM>1) pari_warn(warnmem,"hnfall[1], li = %ld", li);
    1836         [ #  # ]:          0 :           gerepileall(av1, B? 2: 1, &A, &B);
    1837                 :            :         }
    1838                 :            :       }
    1839         [ +  + ]:     835642 :       if (signe( gcoeff(A,li,j) )) break;
    1840                 :     573642 :       h[j] = li-1;
    1841                 :            :     }
    1842         [ -  + ]:     262000 :     if (j == r) continue;
    1843                 :     262000 :     r--;
    1844         [ +  + ]:     262000 :     if (j < r) /* A[j] != 0 */
    1845                 :            :     {
    1846                 :     177882 :       swap(gel(A,j), gel(A,r));
    1847         [ +  + ]:     177882 :       if (B) swap(gel(B,j), gel(B,r));
    1848                 :     177882 :       h[j] = h[r]; h[r] = li; c[li] = r;
    1849                 :            :     }
    1850         [ +  + ]:     262000 :     if (signe(gcoeff(A,li,r)) < 0)
    1851                 :            :     {
    1852                 :      51820 :       ZV_neg_inplace(gel(A,r));
    1853         [ +  + ]:      51820 :       if (B) ZV_togglesign(gel(B,r));
    1854                 :            :     }
    1855                 :     262000 :     ZM_reduce(A,B, li,r);
    1856         [ -  + ]:     262000 :     if (gc_needed(av1,1))
    1857                 :            :     {
    1858         [ #  # ]:          0 :       if (DEBUGMEM>1) pari_warn(warnmem,"hnfall[2], li = %ld", li);
    1859         [ #  # ]:          0 :       gerepileall(av1, B? 2: 1, &A, &B);
    1860                 :            :     }
    1861                 :            :   }
    1862                 :            : 
    1863         [ -  + ]:     103718 :   if (DEBUGLEVEL>5) err_printf("\nhnfall, final phase: ");
    1864                 :     103718 :   r--; /* first r cols are in the image the n-r (independent) last ones */
    1865         [ +  + ]:     622759 :   for (j=1; j<=r; j++)
    1866         [ +  + ]:    1189001 :     for (i=h[j]; i; i--)
    1867                 :            :     {
    1868                 :     669960 :       a = gcoeff(A,i,j);
    1869                 :     669960 :       k = c[i];
    1870         [ +  + ]:     669960 :       if (signe(a)) ZC_elem(a,gcoeff(A,i,k), A,B, j,k);
    1871                 :     669960 :       ZM_reduce(A,B, i,k); /* ensure reduced entries, even if a = 0 */
    1872         [ -  + ]:     669960 :       if (gc_needed(av1,1))
    1873                 :            :       {
    1874         [ #  # ]:          0 :         if (DEBUGMEM>1) pari_warn(warnmem,"hnfall[3], j = %ld", j);
    1875         [ #  # ]:          0 :         gerepileall(av1, B? 2: 1, &A, &B);
    1876                 :            :       }
    1877                 :            :     }
    1878         [ -  + ]:     103718 :   if (DEBUGLEVEL>5) err_printf("\n");
    1879         [ +  + ]:     103718 :   if (remove) remove_0cols(r, &A, &B, remove);
    1880         [ +  + ]:     103718 :   gerepileall(av, B? 2: 1, &A, &B);
    1881         [ +  + ]:     103718 :   if (B) *ptB = B;
    1882                 :     103725 :   return A;
    1883                 :            : }
    1884                 :            : 
    1885                 :            : GEN
    1886                 :         14 : hnfall(GEN x)
    1887                 :            : {
    1888                 :         14 :   GEN z = cgetg(3, t_VEC);
    1889                 :         14 :   gel(z,1) = ZM_hnfall(x, (GEN*)(z+2), 1);
    1890                 :         14 :   return z;
    1891                 :            : }
    1892                 :            : GEN
    1893                 :          0 : hnf(GEN x) { return mathnf0(x,0); }
    1894                 :            : 
    1895                 :            : /* C = A^(-1)(tB) where A, B, C are integral, A is upper triangular, t t_INT */
    1896                 :            : GEN
    1897                 :      25508 : hnf_divscale(GEN A, GEN B, GEN t)
    1898                 :            : {
    1899                 :      25508 :   long n = lg(A)-1, i,j,k;
    1900                 :      25508 :   GEN m, c = cgetg(n+1,t_MAT);
    1901                 :            : 
    1902         [ -  + ]:      25508 :   if (!n) return c;
    1903         [ +  + ]:      96124 :   for (k=1; k<=n; k++)
    1904                 :            :   {
    1905                 :      70616 :     GEN u = cgetg(n+1, t_COL), b = gel(B,k);
    1906                 :      70616 :     pari_sp av = avma;
    1907                 :      70616 :     gel(c,k) = u; m = mulii(gel(b,n),t);
    1908                 :      70616 :     gel(u,n) = gerepileuptoint(av, diviiexact(m, gcoeff(A,n,n)));
    1909         [ +  + ]:     349230 :     for (i=n-1; i>0; i--)
    1910                 :            :     {
    1911                 :     278614 :       av = avma; m = mulii(gel(b,i),t);
    1912         [ +  + ]:    3179855 :       for (j=i+1; j<=n; j++) m = subii(m, mulii(gcoeff(A,i,j),gel(u,j)));
    1913                 :     278614 :       gel(u,i) = gerepileuptoint(av, diviiexact(m, gcoeff(A,i,i)));
    1914                 :            :     }
    1915                 :            :   }
    1916                 :      25508 :   return c;
    1917                 :            : }
    1918                 :            : 
    1919                 :            : /* A, B integral upper HNF. A^(-1) B integral ? */
    1920                 :            : int
    1921                 :         56 : hnfdivide(GEN A, GEN B)
    1922                 :            : {
    1923                 :         56 :   pari_sp av = avma;
    1924                 :         56 :   long n = lg(A)-1, i,j,k;
    1925                 :            :   GEN u, b, m, r;
    1926                 :            : 
    1927         [ -  + ]:         56 :   if (!n) return 1;
    1928         [ -  + ]:         56 :   if (lg(B)-1 != n) pari_err_DIM("hnfdivide");
    1929                 :         56 :   u = cgetg(n+1, t_COL);
    1930         [ +  + ]:        273 :   for (k=1; k<=n; k++)
    1931                 :            :   {
    1932                 :        217 :     b = gel(B,k);
    1933                 :        217 :     m = gel(b,k);
    1934                 :        217 :     gel(u,k) = dvmdii(m, gcoeff(A,k,k), &r);
    1935         [ -  + ]:        217 :     if (r != gen_0) { avma = av; return 0; }
    1936         [ +  + ]:        588 :     for (i=k-1; i>0; i--)
    1937                 :            :     {
    1938                 :        371 :       m = gel(b,i);
    1939         [ +  + ]:       1092 :       for (j=i+1; j<=k; j++) m = subii(m, mulii(gcoeff(A,i,j),gel(u,j)));
    1940                 :        371 :       m = dvmdii(m, gcoeff(A,i,i), &r);
    1941         [ -  + ]:        371 :       if (r != gen_0) { avma = av; return 0; }
    1942                 :        371 :       gel(u,i) = m;
    1943                 :            :     }
    1944                 :            :   }
    1945                 :         56 :   avma = av; return 1;
    1946                 :            : }
    1947                 :            : 
    1948                 :            : /* A upper HNF, b integral vector. Return A^(-1) b if integral,
    1949                 :            :  * NULL otherwise. Assume #A[,1] = #b. */
    1950                 :            : GEN
    1951                 :      88116 : hnf_invimage(GEN A, GEN b)
    1952                 :            : {
    1953                 :      88116 :   pari_sp av = avma;
    1954                 :      88116 :   long n = lg(A)-1, m, i, k;
    1955                 :            :   GEN u, r;
    1956                 :            : 
    1957         [ -  + ]:      88116 :   if (!n) return NULL;
    1958                 :      88116 :   m = nbrows(A); /* m >= n */
    1959                 :      88116 :   u = cgetg(n+1, t_COL);
    1960         [ +  - ]:     600250 :   for (i = n, k = m; k > 0; k--)
    1961                 :            :   {
    1962                 :     600250 :     pari_sp av2 = avma;
    1963                 :            :     long j;
    1964                 :     600250 :     GEN t = gel(b,k), Aki = gcoeff(A,k,i);
    1965         [ -  + ]:     600250 :     if (typ(t) != t_INT) pari_err_TYPE("hnf_invimage",t);
    1966         [ +  + ]:    2700145 :     for (j=i+1; j<=n; j++) t = subii(t, mulii(gcoeff(A,k,j),gel(u,j)));
    1967         [ -  + ]:     600250 :     if (!signe(Aki))
    1968                 :            :     {
    1969         [ #  # ]:          0 :       if (signe(t)) { avma = av;return NULL; }
    1970                 :          0 :       avma = av2; gel(u,i) = gen_0; continue;
    1971                 :            :     }
    1972                 :     600250 :     t = dvmdii(t, Aki, &r);
    1973         [ +  + ]:     600250 :     if (r != gen_0) { avma = av; return NULL; }
    1974                 :     593075 :     gel(u,i) = gerepileuptoint(av2, t);
    1975         [ +  + ]:     593075 :     if (--i == 0) break;
    1976                 :            :   }
    1977                 :            :   /* If there is a solution, it must be u. Check remaining equations */
    1978         [ +  + ]:     161882 :   for (; k > 0; k--)
    1979                 :            :   {
    1980                 :      80941 :     pari_sp av2 = avma;
    1981                 :            :     long j;
    1982                 :      80941 :     GEN t = gel(b,k);
    1983         [ -  + ]:      80941 :     if (typ(t) != t_INT) pari_err_TYPE("hnf_invimage",t);
    1984         [ +  + ]:     666659 :     for (j=1; j<=n; j++) t = subii(t, mulii(gcoeff(A,k,j),gel(u,j)));
    1985         [ -  + ]:      80941 :     if (signe(t)) { avma = av;return NULL; }
    1986                 :      80941 :     avma = av2;
    1987                 :            :   }
    1988                 :      88116 :   return u;
    1989                 :            : }
    1990                 :            : 
    1991                 :            : /* A upper HNF, B integral matrix or column. Return A^(-1) B if integral,
    1992                 :            :  * NULL otherwise */
    1993                 :            : GEN
    1994                 :      21707 : hnf_solve(GEN A, GEN B)
    1995                 :            : {
    1996                 :            :   pari_sp av;
    1997                 :            :   long i, l;
    1998                 :            :   GEN C;
    1999                 :            : 
    2000         [ +  + ]:      21707 :   if (typ(B) == t_COL) return hnf_invimage(A, B);
    2001                 :      18263 :   av = avma; C = cgetg_copy(B, &l);
    2002         [ +  + ]:      97076 :   for (i = 1; i < l; i++) {
    2003                 :      82649 :     GEN c = hnf_invimage(A, gel(B,i));
    2004         [ +  + ]:      82649 :     if (!c) { avma = av; return NULL; }
    2005                 :      78813 :     gel(C,i) = c;
    2006                 :            :   }
    2007                 :      21707 :   return C;
    2008                 :            : }
    2009                 :            : 
    2010                 :            : /***************************************************************/
    2011                 :            : /**                                                           **/
    2012                 :            : /**               SMITH NORMAL FORM REDUCTION                 **/
    2013                 :            : /**                                                           **/
    2014                 :            : /***************************************************************/
    2015                 :            : 
    2016                 :            : static GEN
    2017                 :          0 : trivsmith(long all)
    2018                 :            : {
    2019                 :            :   GEN z;
    2020         [ #  # ]:          0 :   if (!all) return cgetg(1,t_VEC);
    2021                 :          0 :   z=cgetg(4,t_VEC);
    2022                 :          0 :   gel(z,1) = cgetg(1,t_MAT);
    2023                 :          0 :   gel(z,2) = cgetg(1,t_MAT);
    2024                 :          0 :   gel(z,3) = cgetg(1,t_MAT); return z;
    2025                 :            : }
    2026                 :            : 
    2027                 :            : static void
    2028                 :      35791 : snf_pile(pari_sp av, GEN *x, GEN *U, GEN *V)
    2029                 :            : {
    2030                 :            :   GEN *gptr[3];
    2031                 :      35791 :   int c = 1; gptr[0]=x;
    2032         [ +  + ]:      35791 :   if (*U) gptr[c++] = U;
    2033         [ +  + ]:      35791 :   if (*V) gptr[c++] = V;
    2034                 :      35791 :   gerepilemany(av,gptr,c);
    2035                 :      35791 : }
    2036                 :            : 
    2037                 :            : static GEN
    2038                 :      73115 : bezout_step(GEN *pa, GEN *pb, GEN *pu, GEN *pv)
    2039                 :            : {
    2040                 :      73115 :   GEN a = *pa, b = *pb, d;
    2041         [ +  + ]:      73115 :   if (absi_equal(a,b))
    2042                 :            :   {
    2043                 :      29078 :     long sa = signe(a), sb = signe(b);
    2044                 :      29078 :     *pv = gen_0;
    2045         [ +  + ]:      29078 :     if (sb == sa) {
    2046                 :      28903 :       *pa = *pb = gen_1;
    2047         [ +  - ]:      28903 :       if (sa > 0) {*pu=gen_1; return a;} else {*pu=gen_m1; return absi(a);}
    2048                 :            :     }
    2049         [ +  - ]:        175 :     if (sa > 0) { *pa = *pu = gen_1; *pb = gen_m1; return a; }
    2050                 :          0 :     *pa = *pu = gen_m1; *pb = gen_1; return b;
    2051                 :            :   }
    2052                 :      44037 :   d = bezout(a,b, pu,pv);
    2053                 :      44037 :   *pa = diviiexact(a, d);
    2054                 :      73115 :   *pb = diviiexact(b, d); return d;
    2055                 :            : }
    2056                 :            : 
    2057                 :            : static int
    2058                 :     224994 : negcmpii(void *E, GEN x, GEN y) { (void)E; return -cmpii(x,y); }
    2059                 :            : 
    2060                 :            : /* does b = x[i,i] divide all entries in x[1..i-1,1..i-1] ? If so, return 0;
    2061                 :            :  * else return the index of a problematic row */
    2062                 :            : static long
    2063                 :     135556 : ZM_snf_no_divide(GEN x, long i)
    2064                 :            : {
    2065                 :     135556 :   GEN b = gcoeff(x,i,i);
    2066                 :            :   long j, k;
    2067                 :            : 
    2068         [ -  + ]:     135556 :   if (!signe(b))
    2069                 :            :   { /* impossible in the current implementation : x square of maximal rank */
    2070         [ #  # ]:          0 :     for (k = 1; k < i; k++)
    2071         [ #  # ]:          0 :       for (j = 1; j < i; j++)
    2072         [ #  # ]:          0 :         if (signe(gcoeff(x,k,j))) return k;
    2073                 :          0 :     return 0;
    2074                 :            :   }
    2075         [ +  + ]:     135556 :   if (is_pm1(b)) return 0;
    2076         [ +  + ]:     285447 :   for (k=1; k<i; k++)
    2077                 :            :   {
    2078         [ +  + ]:    1063588 :     for (j=1; j<i; j++)
    2079         [ +  + ]:     854533 :       if (signe(remii(gcoeff(x,k,j),b))) return k;
    2080                 :            :   }
    2081                 :     135556 :   return 0;
    2082                 :            : }
    2083                 :            : 
    2084                 :            : /* Return the SNF D of matrix X. If ptU/ptV non-NULL set them to U/V
    2085                 :            :  * to that D = UXV */
    2086                 :            : GEN
    2087                 :      36596 : ZM_snfall_i(GEN x, GEN *ptU, GEN *ptV, int return_vec)
    2088                 :            : {
    2089                 :      36596 :   pari_sp av0 = avma, av;
    2090                 :            :   long i, j, k, m0, m, n0, n;
    2091                 :      36596 :   GEN p1, u, v, U, V, V0, mdet, ys, perm = NULL;
    2092                 :            : 
    2093                 :      36596 :   n0 = n = lg(x)-1;
    2094         [ +  + ]:      36596 :   if (!n) {
    2095         [ +  + ]:        805 :     if (ptU) *ptU = cgetg(1,t_MAT);
    2096         [ +  + ]:        805 :     if (ptV) *ptV = cgetg(1,t_MAT);
    2097         [ +  + ]:        805 :     return cgetg(1, return_vec? t_VEC: t_MAT);
    2098                 :            :   }
    2099                 :      35791 :   av = avma;
    2100                 :      35791 :   m0 = m = nbrows(x);
    2101                 :            : 
    2102         [ +  + ]:      35791 :   U = ptU? gen_1: NULL; /* TRANSPOSE of row transform matrix [act on columns]*/
    2103         [ +  + ]:      35791 :   V = ptV? gen_1: NULL;
    2104                 :      35791 :   V0 = NULL;
    2105                 :      35791 :   x = RgM_shallowcopy(x);
    2106 [ +  + ][ +  + ]:      35791 :   if (m == n && ZM_ishnf(x))
    2107                 :            :   {
    2108                 :      29568 :     mdet = ZM_det_triangular(x);
    2109         [ +  + ]:      29568 :     if (V) *ptV = matid(n);
    2110                 :            :   }
    2111                 :            :   else
    2112                 :            :   {
    2113                 :       6223 :     mdet = ZM_detmult(x);
    2114         [ +  + ]:       6223 :     if (signe(mdet))
    2115                 :            :     {
    2116         [ +  + ]:       6209 :       if (!V)
    2117                 :         63 :         p1 = ZM_hnfmod(x,mdet);
    2118                 :            :       else
    2119                 :            :       {
    2120         [ +  + ]:       6146 :         if (m == n)
    2121                 :            :         {
    2122                 :       6132 :           p1 = ZM_hnfmod(x,mdet);
    2123                 :       6132 :           *ptV = RgM_solve(x,p1);
    2124                 :            :         }
    2125                 :            :         else
    2126         [ +  - ]:         14 :           p1 = ZM_hnfperm(x, ptV, ptU? &perm: NULL);
    2127                 :            :       }
    2128                 :       6209 :       mdet = ZM_det_triangular(p1);
    2129                 :            :     }
    2130                 :            :     else
    2131         [ +  + ]:         14 :       p1 = ZM_hnfperm(x, ptV, ptU? &perm: NULL);
    2132                 :       6223 :     x = p1;
    2133                 :            :   }
    2134                 :      35791 :   n = lg(x)-1;
    2135         [ +  + ]:      35791 :   if (V)
    2136                 :            :   {
    2137                 :      26285 :     V = *ptV;
    2138         [ +  + ]:      26285 :     if (n != n0)
    2139                 :            :     {
    2140                 :         14 :       V0 = vecslice(V, 1, n0 - n); /* kernel */
    2141                 :         14 :       V  = vecslice(V, n0-n+1, n0);
    2142                 :         14 :       av = avma;
    2143                 :            :     }
    2144                 :            :   }
    2145                 :            :   /* independent columns */
    2146         [ +  + ]:      35791 :   if (!signe(mdet))
    2147                 :            :   {
    2148         [ +  - ]:         14 :     if (n)
    2149                 :            :     {
    2150                 :         14 :       x = ZM_snfall_i(shallowtrans(x), ptV, ptU, return_vec); /* swap ptV,ptU */
    2151 [ +  + ][ +  - ]:         14 :       if (typ(x) == t_MAT && n != m) x = shallowtrans(x);
    2152         [ +  + ]:         14 :       if (V) V = ZM_mul(V, shallowtrans(*ptV));
    2153         [ +  + ]:         14 :       if (U) U = *ptU; /* TRANSPOSE */
    2154                 :            :     }
    2155                 :            :     else /* 0 matrix */
    2156                 :            :     {
    2157                 :          0 :       x = cgetg(1,t_MAT);
    2158         [ #  # ]:          0 :       if (V) V = cgetg(1, t_MAT);
    2159         [ #  # ]:          0 :       if (U) U = matid(m);
    2160                 :            :     }
    2161                 :         14 :     goto THEEND;
    2162                 :            :   }
    2163         [ +  + ]:      35777 :   if (U) U = matid(n);
    2164                 :            : 
    2165                 :            :   /* square, maximal rank n */
    2166                 :      35777 :   p1 = gen_indexsort(RgM_diagonal_shallow(x), NULL, &negcmpii);
    2167                 :      35777 :   ys = cgetg(n+1,t_MAT);
    2168         [ +  + ]:     204778 :   for (j=1; j<=n; j++) gel(ys,j) = vecpermute(gel(x, p1[j]), p1);
    2169                 :      35777 :   x = ys;
    2170         [ +  + ]:      35777 :   if (U) U = vecpermute(U, p1);
    2171         [ +  + ]:      35777 :   if (V) V = vecpermute(V, p1);
    2172                 :            : 
    2173                 :      35777 :   p1 = ZM_hnfmod(x, mdet);
    2174         [ +  + ]:      35777 :   if (V) V = ZM_mul(V, RgM_solve(x,p1));
    2175                 :      35777 :   x = p1;
    2176                 :            : 
    2177         [ -  + ]:      35777 :   if (DEBUGLEVEL>7) err_printf("starting SNF loop");
    2178         [ +  + ]:     169001 :   for (i=n; i>1; i--)
    2179                 :            :   {
    2180         [ -  + ]:     133224 :     if (DEBUGLEVEL>7) err_printf("\ni = %ld: ",i);
    2181                 :            :     for(;;)
    2182                 :            :     {
    2183                 :     170213 :       int c = 0;
    2184                 :            :       GEN a, b;
    2185         [ +  + ]:     906559 :       for (j=i-1; j>=1; j--)
    2186                 :            :       {
    2187         [ +  + ]:     736346 :         b = gcoeff(x,i,j); if (!signe(b)) continue;
    2188                 :       5216 :         a = gcoeff(x,i,i);
    2189                 :       5216 :         ZC_elem(b, a, x,V, j,i);
    2190         [ -  + ]:       5216 :         if (gc_needed(av,1))
    2191                 :            :         {
    2192         [ #  # ]:          0 :           if (DEBUGMEM>1) pari_warn(warnmem,"[1]: ZM_snfall i = %ld", i);
    2193                 :          0 :           snf_pile(av, &x,&U,&V);
    2194                 :            :         }
    2195                 :            :       }
    2196         [ -  + ]:     170213 :       if (DEBUGLEVEL>7) err_printf("; ");
    2197         [ +  + ]:     906559 :       for (j=i-1; j>=1; j--)
    2198                 :            :       {
    2199                 :            :         GEN d;
    2200         [ +  + ]:     736346 :         b = gcoeff(x,j,i); if (!signe(b)) continue;
    2201                 :      73115 :         a = gcoeff(x,i,i);
    2202                 :      73115 :         d = bezout_step(&a, &b, &u, &v);
    2203         [ +  + ]:     600831 :         for (k = 1; k < i; k++)
    2204                 :            :         {
    2205                 :     527716 :           GEN t = addii(mulii(u,gcoeff(x,i,k)),mulii(v,gcoeff(x,j,k)));
    2206                 :     527716 :           gcoeff(x,j,k) = subii(mulii(a,gcoeff(x,j,k)),
    2207                 :     527716 :                                 mulii(b,gcoeff(x,i,k)));
    2208                 :     527716 :           gcoeff(x,i,k) = t;
    2209                 :            :         }
    2210                 :      73115 :         gcoeff(x,j,i) = gen_0;
    2211                 :      73115 :         gcoeff(x,i,i) = d;
    2212         [ +  + ]:      73115 :         if (U) update(u,v,a,b,(GEN*)(U+i),(GEN*)(U+j));
    2213         [ -  + ]:      73115 :         if (gc_needed(av,1))
    2214                 :            :         {
    2215         [ #  # ]:          0 :           if (DEBUGMEM>1) pari_warn(warnmem,"[2]: ZM_snfall, i = %ld", i);
    2216                 :          0 :           snf_pile(av, &x,&U,&V);
    2217                 :            :         }
    2218                 :      73115 :         c = 1;
    2219                 :            :       }
    2220         [ +  + ]:     170213 :       if (!c)
    2221                 :            :       {
    2222                 :     135556 :         k = ZM_snf_no_divide(x, i);
    2223         [ +  + ]:     135556 :         if (!k) break;
    2224                 :            : 
    2225                 :            :         /* x[k,j] != 0 mod b */
    2226         [ +  + ]:       9776 :         for (j=1; j<=i; j++)
    2227                 :       7444 :           gcoeff(x,i,j) = addii(gcoeff(x,i,j),gcoeff(x,k,j));
    2228         [ +  + ]:       2332 :         if (U) gel(U,i) = gadd(gel(U,i),gel(U,k));
    2229                 :            :       }
    2230         [ -  + ]:      36989 :       if (gc_needed(av,1))
    2231                 :            :       {
    2232         [ #  # ]:          0 :         if (DEBUGMEM>1) pari_warn(warnmem,"[3]: ZM_snfall");
    2233                 :          0 :         snf_pile(av, &x,&U,&V);
    2234                 :            :       }
    2235                 :      36989 :     }
    2236                 :            :   }
    2237         [ -  + ]:      35777 :   if (DEBUGLEVEL>7) err_printf("\n");
    2238         [ +  + ]:     204778 :   for (k=1; k<=n; k++)
    2239         [ +  + ]:     169001 :     if (signe(gcoeff(x,k,k)) < 0)
    2240                 :            :     {
    2241         [ +  + ]:         63 :       if (V) ZV_togglesign(gel(V,k));
    2242                 :         63 :       togglesign(gcoeff(x,k,k));
    2243                 :            :     }
    2244                 :            : THEEND:
    2245         [ +  + ]:      35791 :   if (return_vec)
    2246                 :            :   {
    2247                 :      34601 :     long l = lg(x)-1;
    2248         [ +  + ]:      34601 :     if (typ(x) == t_MAT) x = RgM_diagonal_shallow(x);
    2249         [ +  + ]:      34601 :     if (m0 > l) x = shallowconcat(zerovec(m0-l), x);
    2250                 :            :   }
    2251                 :            : 
    2252         [ +  + ]:      35791 :   if (V0)
    2253                 :            :   {
    2254         [ +  - ]:         14 :     if (!return_vec) x = shallowconcat(zeromat(m,n0-n), x);
    2255         [ +  - ]:         14 :     if (V) V = shallowconcat(V0, V);
    2256                 :            :   }
    2257         [ +  + ]:      35791 :   if (U)
    2258                 :            :   {
    2259                 :      30247 :     U = shallowtrans(U);
    2260         [ +  + ]:      30247 :     if (perm) U = vecpermute(U, perm_inv(perm));
    2261                 :            :   }
    2262                 :      35791 :   snf_pile(av0, &x,&U,&V);
    2263         [ +  + ]:      35791 :   if (ptU) *ptU = U;
    2264         [ +  + ]:      35791 :   if (ptV) *ptV = V;
    2265                 :      36596 :   return x;
    2266                 :            : }
    2267                 :            : GEN
    2268                 :       1379 : ZM_snfall(GEN x, GEN *U, GEN *V) { return ZM_snfall_i(x, U, V, 0); }
    2269                 :            : GEN
    2270                 :        217 : ZM_snf(GEN x) { return ZM_snfall_i(x, NULL,NULL, 1); }
    2271                 :            : 
    2272                 :            : GEN
    2273                 :         28 : smith(GEN x) {
    2274         [ -  + ]:         28 :   if (typ(x)!=t_MAT) pari_err_TYPE("smith",x);
    2275                 :         28 :   RgM_check_ZM(x, "smith");
    2276                 :         28 :   return ZM_snfall_i(x, NULL,NULL, 1);
    2277                 :            : }
    2278                 :            : GEN
    2279                 :        301 : smithall(GEN x)
    2280                 :            : {
    2281                 :        301 :   GEN z = cgetg(4, t_VEC);
    2282         [ -  + ]:        301 :   if (typ(x)!=t_MAT) pari_err_TYPE("smithall",x);
    2283                 :        301 :   RgM_check_ZM(x, "smithall");
    2284                 :        301 :   gel(z,3) = ZM_snfall_i(x, (GEN*)(z+1),(GEN*)(z+2), 0);
    2285                 :        301 :   return z;
    2286                 :            : }
    2287                 :            : 
    2288                 :            : void
    2289                 :      31570 : ZM_snfclean(GEN d, GEN u, GEN v)
    2290                 :            : {
    2291                 :      31570 :   long i, c, l = lg(d);
    2292                 :            : 
    2293         [ +  - ]:      31570 :   if (typ(d) == t_VEC)
    2294 [ +  + ][ +  + ]:     135261 :     for (c=1; c<l; c++) { GEN t = gel(d,c); if (is_pm1(t)) break; }
    2295                 :            :   else
    2296                 :            :   {
    2297 [ #  # ][ #  # ]:          0 :     for (c=1; c<l; c++) { GEN t = gcoeff(d,c,c); if (is_pm1(t)) break; }
    2298 [ #  # ][ #  # ]:          0 :     if (c < l) for (i = 1; i < c; i++) setlg(gel(d,i), c);
    2299                 :            :   }
    2300                 :      31570 :   setlg(d, c);
    2301 [ +  + ][ +  + ]:     179389 :   if (u) for (i=1; i<l; i++) setlg(gel(u,i), c);
    2302         [ +  + ]:      31570 :   if (v) setlg(v, c);
    2303                 :      31570 : }
    2304                 :            : 
    2305                 :            : /* Assume z was computed by [g]smithall(). Remove the 1s on the diagonal */
    2306                 :            : GEN
    2307                 :        525 : smithclean(GEN z)
    2308                 :            : {
    2309                 :            :   long i, j, h, l, c, d;
    2310                 :            :   GEN U, V, y, D, t;
    2311                 :            : 
    2312         [ -  + ]:        525 :   if (typ(z) != t_VEC) pari_err_TYPE("smithclean",z);
    2313         [ +  + ]:        525 :   l = lg(z); if (l == 1) return cgetg(1,t_VEC);
    2314                 :        518 :   U = gel(z,1);
    2315 [ +  + ][ +  + ]:        518 :   if (l != 4 || typ(U) != t_MAT)
    2316                 :            :   { /* assume z = vector of elementary divisors */
    2317         [ +  + ]:        630 :     for (c=1; c<l; c++)
    2318         [ +  + ]:        581 :       if (gequal1(gel(z,c))) break;
    2319                 :        217 :     return gcopy_lg(z, c);
    2320                 :            :   }
    2321                 :        301 :   V = gel(z,2);
    2322                 :        301 :   D = gel(z,3);
    2323                 :        301 :   l = lg(D);
    2324         [ -  + ]:        301 :   if (l == 1) return gcopy(z);
    2325                 :        301 :   h = lgcols(D);
    2326         [ +  + ]:        301 :   if (h > l)
    2327                 :            :   { /* D = vconcat(zero matrix, diagonal matrix) */
    2328         [ +  - ]:         14 :     for (c=1+h-l, d=1; c<h; c++,d++)
    2329         [ +  - ]:          7 :       if (gequal1(gcoeff(D,c,d))) break;
    2330                 :            :   }
    2331         [ +  + ]:        294 :   else if (h < l)
    2332                 :            :   { /* D = concat(zero matrix, diagonal matrix) */
    2333         [ +  - ]:         14 :     for (c=1, d=1+l-h; d<l; c++,d++)
    2334         [ +  - ]:          7 :       if (gequal1(gcoeff(D,c,d))) break;
    2335                 :            :   }
    2336                 :            :   else
    2337                 :            :   { /* D diagonal */
    2338         [ +  + ]:        644 :     for (c=1; c<l; c++)
    2339         [ +  + ]:        392 :       if (gequal1(gcoeff(D,c,c))) break;
    2340                 :        287 :     d = c;
    2341                 :            :   }
    2342                 :            :   /* U was (h-1)x(h-1), V was (l-1)x(l-1), D was (h-1)x(l-1) */
    2343                 :        301 :   y = cgetg(4,t_VEC);
    2344                 :            :   /* truncate U to (c-1) x (h-1) */
    2345                 :        301 :   gel(y,1) = t = cgetg(h,t_MAT);
    2346         [ +  + ]:        721 :   for (j=1; j<h; j++) gel(t,j) = gcopy_lg(gel(U,j), c);
    2347                 :            :   /* truncate V to (l-1) x (d-1) */
    2348                 :        301 :   gel(y,2) = gcopy_lg(V, d);
    2349                 :        301 :   gel(y,3) = t = zeromatcopy(c-1, d-1);
    2350                 :            :   /* truncate D to a (c-1) x (d-1) matrix */
    2351         [ +  + ]:        301 :   if (d > 1)
    2352                 :            :   {
    2353         [ -  + ]:        294 :     if (h > l)
    2354                 :            :     {
    2355         [ #  # ]:          0 :       for (i=1+h-l, j=1; i<h; i++,j++)
    2356                 :          0 :         gcoeff(t,i,j) = gcopy(gcoeff(D,i,j));
    2357                 :            :     }
    2358         [ +  + ]:        294 :     else if (h < l)
    2359                 :            :     {
    2360         [ -  + ]:          7 :       for (i=1, j=1+l-h; j<d; i++,j++)
    2361                 :          0 :         gcoeff(t,i,j) = gcopy(gcoeff(D,i,j));
    2362                 :            :     }
    2363                 :            :     else
    2364                 :            :     {
    2365         [ +  + ]:        644 :       for (j=1; j<d; j++)
    2366                 :        357 :         gcoeff(t,j,j) = gcopy(gcoeff(D,j,j));
    2367                 :            :     }
    2368                 :            :   }
    2369                 :        525 :   return y;
    2370                 :            : }
    2371                 :            : 
    2372                 :            : /* does b = x[i,i] divide all entries in x[1..i-1,1..i-1] ? If so, return 0;
    2373                 :            :  * else return the index of a problematic row */
    2374                 :            : static long
    2375                 :        112 : gsnf_no_divide(GEN x, long i, long vx)
    2376                 :            : {
    2377                 :        112 :   GEN b = gcoeff(x,i,i);
    2378                 :            :   long j, k;
    2379                 :            : 
    2380         [ +  + ]:        112 :   if (gequal0(b))
    2381                 :            :   {
    2382         [ +  - ]:         14 :     for (k = 1; k < i; k++)
    2383         [ +  - ]:         14 :       for (j = 1; j < i; j++)
    2384         [ +  - ]:         14 :         if (!gequal0(gcoeff(x,k,j))) return k;
    2385                 :          0 :     return 0;
    2386                 :            :   }
    2387                 :            : 
    2388 [ +  - ][ +  + ]:         98 :   if (!is_RgX(b,vx) || degpol(b)<=0) return 0;
    2389         [ +  + ]:        133 :   for (k = 1; k < i; k++)
    2390         [ +  + ]:        231 :     for (j = 1; j < i; j++)
    2391                 :            :     {
    2392                 :        154 :       GEN z = gcoeff(x,k,j), r;
    2393         [ +  + ]:        154 :       if (!is_RgX(z,vx)) z = scalarpol(z, vx);
    2394                 :        154 :       r = RgX_rem(z, b);
    2395         [ +  + ]:        154 :       if (signe(r) && (! isinexactreal(r) ||
           [ -  +  #  # ]
    2396                 :          0 :              gexpo(r) > 16 + gexpo(b) - prec2nbits(gprecision(r)))
    2397                 :          7 :          ) return k;
    2398                 :            :     }
    2399                 :        112 :   return 0;
    2400                 :            : }
    2401                 :            : 
    2402                 :            : /* Hermite Normal Form, with base change matrix if ptB != NULL.
    2403                 :            :  * If 'remove' = 1, remove 0 columns (do NOT update *ptB accordingly)
    2404                 :            :  * If 'remove' = 2, remove 0 columns and update *ptB accordingly */
    2405                 :            : GEN
    2406                 :         21 : RgM_hnfall(GEN A, GEN *pB, long remove)
    2407                 :            : {
    2408                 :            :   pari_sp av;
    2409                 :            :   long li, j, k, m, n, def, ldef;
    2410                 :            :   GEN B;
    2411                 :         21 :   long vx = gvar(A);
    2412                 :            : 
    2413                 :         21 :   n = lg(A)-1;
    2414 [ +  - ][ -  + ]:         21 :   if (vx==NO_VARIABLE || !n)
    2415                 :            :   {
    2416                 :          0 :     RgM_check_ZM(A, "mathnf0");
    2417                 :          0 :     return ZM_hnfall(A, pB, remove);
    2418                 :            :   }
    2419                 :         21 :   m = nbrows(A);
    2420                 :         21 :   av = avma;
    2421                 :         21 :   A = RgM_shallowcopy(A);
    2422         [ +  + ]:         21 :   B = pB? matid(n): NULL;
    2423         [ -  + ]:         21 :   def = n; ldef = (m>n)? m-n: 0;
    2424         [ +  + ]:         49 :   for (li=m; li>ldef; li--)
    2425                 :            :   {
    2426                 :            :     GEN T;
    2427         [ +  + ]:         49 :     for (j=def-1; j; j--)
    2428                 :            :     {
    2429                 :         21 :       GEN a = gcoeff(A,li,j);
    2430         [ +  + ]:         21 :       if (gequal0(a)) continue;
    2431                 :            : 
    2432         [ -  + ]:          7 :       k = (j==1)? def: j-1;
    2433                 :          7 :       RgC_elem(a,gcoeff(A,li,k), A,B, j,k, li, vx);
    2434                 :            :     }
    2435                 :         28 :     T = gcoeff(A,li,def);
    2436         [ -  + ]:         28 :     if (gequal0(T))
    2437         [ #  # ]:          0 :     { if (ldef) ldef--; }
    2438                 :            :     else
    2439                 :            :     {
    2440                 :            :       GEN d;
    2441                 :         28 :       gcoeff(A,li,def) = RgX_normalize_all(T, &d);
    2442 [ +  + ][ -  + ]:         28 :       if (B && !gequal1(d)) gel(B, def) = RgC_Rg_div(gel(B, def), d);
    2443                 :         28 :       RgM_reduce(A, B, li, def, vx);
    2444                 :         28 :       def--;
    2445                 :            :     }
    2446         [ -  + ]:         28 :     if (gc_needed(av,1))
    2447                 :            :     {
    2448         [ #  # ]:          0 :       if (DEBUGMEM>1) pari_warn(warnmem,"ghnfall");
    2449         [ #  # ]:          0 :       gerepileall(av, B? 2: 1, &A, &B);
    2450                 :            :     }
    2451                 :            :   }
    2452                 :            :   /* rank A = n - def */
    2453         [ +  - ]:         21 :   if (remove) remove_0cols(def, &A, &B, remove);
    2454         [ +  + ]:         21 :   gerepileall(av, B? 2: 1, &A, &B);
    2455         [ +  + ]:         21 :   if (B) *pB = B;
    2456                 :         21 :   return A;
    2457                 :            : }
    2458                 :            : 
    2459                 :            : static GEN
    2460                 :         42 : gsmithall_i(GEN x,long all)
    2461                 :            : {
    2462                 :            :   pari_sp av;
    2463                 :            :   long i, j, k, n;
    2464                 :            :   GEN z, u, v, U, V;
    2465                 :         42 :   long vx = gvar(x);
    2466         [ -  + ]:         42 :   if (typ(x)!=t_MAT) pari_err_TYPE("gsmithall",x);
    2467 [ +  + ][ -  + ]:         42 :   if (vx==NO_VARIABLE) return all? smithall(x): smith(x);
    2468                 :         35 :   n = lg(x)-1;
    2469         [ -  + ]:         35 :   if (!n) return trivsmith(all);
    2470         [ -  + ]:         35 :   if (lgcols(x) != n+1) pari_err_DIM("gsmithall");
    2471                 :         35 :   av = avma;
    2472                 :         35 :   x = RgM_shallowcopy(x);
    2473         [ -  + ]:         35 :   if (all) { U = matid(n); V = matid(n); }
    2474         [ +  + ]:        126 :   for (i=n; i>=2; i--)
    2475                 :            :   {
    2476                 :            :     for(;;)
    2477                 :            :     {
    2478                 :            :       GEN a, b, d;
    2479                 :        189 :       int c = 0;
    2480         [ +  + ]:        546 :       for (j=i-1; j>=1; j--)
    2481                 :            :       {
    2482         [ +  + ]:        357 :         b = gcoeff(x,i,j); if (gequal0(b)) continue;
    2483                 :        140 :         a = gcoeff(x,i,i);
    2484                 :        140 :         d = gbezout_step(&b, &a, &v, &u, vx);
    2485         [ +  + ]:        469 :         for (k = 1; k < i; k++)
    2486                 :            :         {
    2487                 :        329 :           GEN t = gadd(gmul(u,gcoeff(x,k,i)),gmul(v,gcoeff(x,k,j)));
    2488                 :        329 :           gcoeff(x,k,j) = gsub(gmul(a,gcoeff(x,k,j)),gmul(b,gcoeff(x,k,i)));
    2489                 :        329 :           gcoeff(x,k,i) = t;
    2490                 :            :         }
    2491                 :        140 :         gcoeff(x,i,j) = gen_0;
    2492                 :        140 :         gcoeff(x,i,i) = d;
    2493         [ -  + ]:        140 :         if (all) update(u,v,a,b,(GEN*)(V+i),(GEN*)(V+j));
    2494                 :            :       }
    2495         [ +  + ]:        546 :       for (j=i-1; j>=1; j--)
    2496                 :            :       {
    2497         [ +  + ]:        357 :         b = gcoeff(x,j,i); if (gequal0(b)) continue;
    2498                 :        119 :         a = gcoeff(x,i,i);
    2499                 :        119 :         d = gbezout_step(&b, &a, &v, &u, vx);
    2500         [ +  + ]:        406 :         for (k = 1; k < i; k++)
    2501                 :            :         {
    2502                 :        287 :           GEN t = gadd(gmul(u,gcoeff(x,i,k)),gmul(v,gcoeff(x,j,k)));
    2503                 :        287 :           gcoeff(x,j,k) = gsub(gmul(a,gcoeff(x,j,k)),gmul(b,gcoeff(x,i,k)));
    2504                 :        287 :           gcoeff(x,i,k) = t;
    2505                 :            :         }
    2506                 :        119 :         gcoeff(x,j,i) = gen_0;
    2507                 :        119 :         gcoeff(x,i,i) = d;
    2508         [ -  + ]:        119 :         if (all) update(u,v,a,b,(GEN*)(U+i),(GEN*)(U+j));
    2509                 :        119 :         c = 1;
    2510                 :            :       }
    2511         [ +  + ]:        189 :       if (!c)
    2512                 :            :       {
    2513                 :        112 :         k = gsnf_no_divide(x, i, vx);
    2514         [ +  + ]:        112 :         if (!k) break;
    2515                 :            : 
    2516         [ +  + ]:         70 :         for (j=1; j<=i; j++)
    2517                 :         49 :           gcoeff(x,i,j) = gadd(gcoeff(x,i,j),gcoeff(x,k,j));
    2518         [ -  + ]:         21 :         if (all) gel(U,i) = gadd(gel(U,i),gel(U,k));
    2519                 :            :       }
    2520         [ -  + ]:         98 :       if (gc_needed(av,1))
    2521                 :            :       {
    2522         [ #  # ]:          0 :         if (DEBUGMEM>1) pari_warn(warnmem,"gsmithall");
    2523         [ #  # ]:          0 :         gerepileall(av, all? 3: 1, &x, &U, &V);
    2524                 :            :       }
    2525                 :         98 :     }
    2526                 :            :   }
    2527         [ +  + ]:        161 :   for (k=1; k<=n; k++)
    2528                 :            :   {
    2529                 :        126 :     GEN d, T = gcoeff(x,k,k);
    2530         [ +  + ]:        126 :     if (!signe(T)) continue;
    2531         [ +  - ]:        112 :     if (is_RgX(T,vx)) T = RgX_normalize_all(T, &d); else { d = T; T = gen_1; }
    2532 [ -  + ][ #  # ]:        112 :     if (all && !gequal1(d)) gel(V,k) = RgC_Rg_div(gel(V,k), d);
    2533                 :        112 :     gcoeff(x,k,k) = T;
    2534                 :            :   }
    2535         [ -  + ]:         35 :   z = all? mkvec3(shallowtrans(U), V, x): RgM_diagonal_shallow(x);
    2536                 :         42 :   return gerepilecopy(av, z);
    2537                 :            : }
    2538                 :            : 
    2539                 :            : GEN
    2540                 :        364 : matsnf0(GEN x,long flag)
    2541                 :            : {
    2542                 :        364 :   pari_sp av = avma;
    2543         [ -  + ]:        364 :   if (flag > 7) pari_err_FLAG("matsnf");
    2544 [ -  + ][ #  # ]:        364 :   if (typ(x) == t_VEC && flag & 4) return smithclean(x);
    2545 [ +  + ][ -  + ]:        364 :   if (flag & 2) x = flag&1 ? gsmithall(x): gsmith(x);
    2546         [ +  + ]:        322 :   else          x = flag&1 ?  smithall(x):  smith(x);
    2547         [ +  + ]:        364 :   if (flag & 4) x = gerepileupto(av, smithclean(x));
    2548                 :        364 :   return x;
    2549                 :            : }
    2550                 :            : 
    2551                 :            : GEN
    2552                 :         42 : gsmith(GEN x) { return gsmithall_i(x,0); }
    2553                 :            : 
    2554                 :            : GEN
    2555                 :          0 : gsmithall(GEN x) { return gsmithall_i(x,1); }
    2556                 :            : 
    2557                 :            : /* H relation matrix among row of generators g in HNF.  Let URV = D its SNF,
    2558                 :            :  * newU R newV = newD its clean SNF (no 1 in Dnew). Return the diagonal of
    2559                 :            :  * newD, newU and newUi such that  1/U = (newUi, ?).
    2560                 :            :  * Rationale: let (G,0) = g Ui be the new generators then
    2561                 :            :  * 0 = G U R --> G D = 0,  g = G newU  and  G = g newUi */
    2562                 :            : GEN
    2563                 :      31570 : ZM_snf_group(GEN H, GEN *newU, GEN *newUi)
    2564                 :            : {
    2565                 :      31570 :   GEN D = ZM_snfall_i(H, newU, newUi, 1);
    2566                 :            :   long i, j, l;
    2567                 :            : 
    2568 [ +  + ][ +  + ]:      31570 :   ZM_snfclean(D, newU? *newU: NULL, newUi? *newUi: NULL);
    2569                 :      31570 :   l = lg(D);
    2570         [ +  + ]:      31570 :   if (newU) {
    2571                 :      25970 :     GEN U = *newU;
    2572         [ +  + ]:     118790 :     for (i = 1; i < l; i++)
    2573                 :            :     {
    2574                 :      92820 :       GEN d = gel(D,i), d2 = shifti(d, 1);
    2575         [ +  + ]:     810803 :       for (j = 1; j < lg(U); j++)
    2576                 :     717983 :         gcoeff(U,i,j) = centermodii(gcoeff(U,i,j), d, d2);
    2577                 :            :     }
    2578                 :      25970 :     *newU = U;
    2579                 :            :   }
    2580         [ +  + ]:      31570 :   if (newUi) { /* UHV=D -> U^-1 = HVD^-1 -> U^-1 = H(VD^-1 mod 1) mod H */
    2581         [ +  + ]:      25396 :     if (l > 1)
    2582                 :            :     { /* Ui = ZM_inv(U, gen_1); setlg(Ui, l); */
    2583                 :      25074 :       GEN V = FpM_red(*newUi, gel(D,1));
    2584                 :      25074 :       GEN Ui = ZM_mul(H, V);
    2585         [ +  + ]:      96796 :       for (i = 1; i < l; i++) gel(Ui,i) = ZC_Z_divexact(gel(Ui,i), gel(D,i));
    2586                 :      25074 :       *newUi = ZM_hnfrem(Ui, H);
    2587                 :            :     }
    2588                 :            :   }
    2589                 :      31570 :   return D;
    2590                 :            : }
    2591                 :            : 

Generated by: LCOV version 1.9