Code coverage tests

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

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

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

LCOV - code coverage report
Current view: top level - basemath - hnf_snf.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.8.0 lcov report (development 19226-b907b8d) Lines: 1386 1537 90.2 %
Date: 2016-07-29 07:10:27 Functions: 75 81 92.6 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2000  The PARI group.
       2             : 
       3             : This file is part of the PARI/GP package.
       4             : 
       5             : PARI/GP is free software; you can redistribute it and/or modify it under the
       6             : terms of the GNU General Public License as published by the Free Software
       7             : Foundation. It is distributed in the hope that it will be useful, but WITHOUT
       8             : ANY WARRANTY WHATSOEVER.
       9             : 
      10             : Check the License for details. You should have received a copy of it, along
      11             : with the package; see the file 'COPYING'. If not, write to the Free Software
      12             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      13             : 
      14             : #include "pari.h"
      15             : #include "paripriv.h"
      16             : /**************************************************************/
      17             : /**                                                          **/
      18             : /**                HERMITE NORMAL FORM REDUCTION             **/
      19             : /**                                                          **/
      20             : /**************************************************************/
      21             : static GEN ZV_hnfgcdext(GEN A);
      22             : static GEN
      23           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         140 : mathnf0(GEN x, long flag)
      31             : {
      32         140 :   switch(typ(x))
      33             :   {
      34             :     case t_VEC:
      35          70 :       if (RgV_is_ZV(x))
      36          42 :         switch (flag)
      37             :         {
      38             :           case 0:
      39          14 :             if (lg(x) == 1) return cgetg(1, t_MAT);
      40           7 :             retmkmat(mkcol(ZV_content(x)));
      41             :           case 1:
      42             :           case 4:
      43          21 :             return ZV_hnfgcdext(x);
      44             :         }
      45          35 :       x = gtomat(x); break;
      46          70 :     case t_MAT: break;
      47           0 :     default: pari_err_TYPE("mathnf0",x);
      48             :   }
      49             : 
      50         105 :   switch(flag)
      51             :   {
      52          63 :     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           0 :   return NULL; /* not reached */
      61             : }
      62             : 
      63             : /*******************************************************************/
      64             : /*                                                                 */
      65             : /*                SPECIAL HNF (FOR INTERNAL USE !!!)               */
      66             : /*                                                                 */
      67             : /*******************************************************************/
      68             : static int
      69     3054006 : count(GEN mat, long row, long len, long *firstnonzero)
      70             : {
      71     3054006 :   long j, n = 0;
      72             : 
      73   194449577 :   for (j=1; j<=len; j++)
      74             :   {
      75   192114314 :     long p = mael(mat,j,row);
      76   192114314 :     if (p)
      77             :     {
      78     7941579 :       if (labs(p)!=1) return -1;
      79     7222836 :       n++; *firstnonzero=j;
      80             :     }
      81             :   }
      82     2335263 :   return n;
      83             : }
      84             : 
      85             : static int
      86      122419 : count2(GEN mat, long row, long len)
      87             : {
      88             :   long j;
      89     1018086 :   for (j=len; j; j--)
      90      984217 :     if (labs(mael(mat,j,row)) == 1) return j;
      91       33869 :   return 0;
      92             : }
      93             : 
      94             : static GEN
      95       59368 : hnffinal(GEN matgen,GEN perm,GEN* ptdep,GEN* ptB,GEN* ptC)
      96             : {
      97             :   GEN p1,p2,U,H,Hnew,Bnew,Cnew,diagH1;
      98       59368 :   GEN B = *ptB, C = *ptC, dep = *ptdep, depnew;
      99             :   pari_sp av;
     100             :   long i,j,k,s,i1,j1,zc;
     101       59368 :   long co = lg(C);
     102       59368 :   long col = lg(matgen)-1;
     103             :   long lnz, nlze, lig;
     104             : 
     105       59368 :   if (col == 0) return matgen;
     106       59214 :   lnz = nbrows(matgen); /* was called lnz-1 - nr in hnfspec */
     107       59214 :   nlze = nbrows(dep);
     108       59214 :   lig = nlze + lnz;
     109             :   /* H: lnz x lnz [disregarding initial 0 cols], U: col x col */
     110       59214 :   H = ZM_hnflll(matgen, &U, 0);
     111       59214 :   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       59214 :   zc = col - lnz; /* # of 0 columns, correspond to units */
     115       59214 :   if (nlze) { dep = ZM_mul(dep,U); dep += zc; }
     116             : 
     117       59214 :   diagH1 = new_chunk(lnz+1); /* diagH1[i] = 0 iff H[i,i] != 1 (set later) */
     118             : 
     119       59214 :   av = avma;
     120       59214 :   Cnew = cgetg(co, typ(C));
     121       59214 :   setlg(C, col+1); p1 = gmul(C,U);
     122       59214 :   for (j=1; j<=col; j++) gel(Cnew,j) = gel(p1,j);
     123       59214 :   for (   ; j<co ; j++)  gel(Cnew,j) = gel(C,j);
     124             : 
     125             :   /* Clean up B using new H */
     126      254950 :   for (s=0,i=lnz; i; i--)
     127             :   {
     128      195736 :     GEN Di = gel(dep,i), Hi = gel(H,i);
     129      195736 :     GEN h = gel(Hi,i); /* H[i,i] */
     130      195736 :     if ( (diagH1[i] = is_pm1(h)) ) { h = NULL; s++; }
     131     7442536 :     for (j=col+1; j<co; j++)
     132             :     {
     133     7246800 :       GEN z = gel(B,j-col);
     134     7246800 :       p1 = gel(z,i+nlze);
     135     7246800 :       if (h) p1 = truedivii(p1,h);
     136     7246800 :       if (!signe(p1)) continue;
     137     4448578 :       for (k=1; k<=nlze; k++) gel(z,k) = subii(gel(z,k), mulii(p1, gel(Di,k)));
     138     4448578 :       for (   ; k<=lig;  k++) gel(z,k) = subii(gel(z,k), mulii(p1, gel(Hi,k-nlze)));
     139     4448578 :       gel(Cnew,j) = gsub(gel(Cnew,j), gmul(p1, gel(Cnew,i+zc)));
     140             :     }
     141      195736 :     if (gc_needed(av,3))
     142             :     {
     143          71 :       if(DEBUGMEM>1) pari_warn(warnmem,"hnffinal, i = %ld",i);
     144          71 :       gerepileall(av, 2, &Cnew, &B);
     145             :     }
     146             :   }
     147       59214 :   p1 = cgetg(lnz+1,t_VEC); p2 = perm + nlze;
     148      254950 :   for (i1=0, j1=lnz-s, i=1; i<=lnz; i++) /* push the 1 rows down */
     149      195736 :     if (diagH1[i])
     150       98137 :       gel(p1,++j1) = gel(p2,i);
     151             :     else
     152       97599 :       gel(p2,++i1) = gel(p2,i);
     153       59214 :   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       59214 :   lig -= s; col -= s; lnz -= s;
     164       59214 :   Hnew = cgetg(lnz+1,t_MAT);
     165       59214 :   depnew = cgetg(lnz+1,t_MAT); /* only used if nlze > 0 */
     166       59214 :   Bnew = cgetg(co-col,t_MAT);
     167       59214 :   C = shallowcopy(Cnew);
     168      254950 :   for (j=1,i1=j1=0; j<=lnz+s; j++)
     169             :   {
     170      195736 :     GEN z = gel(H,j);
     171      195736 :     if (diagH1[j])
     172             :     { /* hit exactly s times */
     173       98137 :       i1++; C[i1+col] = Cnew[j+zc];
     174       98137 :       p1 = cgetg(lig+1,t_COL); gel(Bnew,i1) = p1;
     175       98137 :       for (i=1; i<=nlze; i++) gel(p1,i) = gcoeff(dep,i,j);
     176       98137 :       p1 += nlze;
     177             :     }
     178             :     else
     179             :     {
     180       97599 :       j1++; C[j1+zc] = Cnew[j+zc];
     181       97599 :       p1 = cgetg(lnz+1,t_COL); gel(Hnew,j1) = p1;
     182       97599 :       depnew[j1] = dep[j];
     183             :     }
     184      842868 :     for (i=k=1; k<=lnz; i++)
     185      647132 :       if (!diagH1[i]) p1[k++] = z[i];
     186             :   }
     187     1264694 :   for (j=s+1; j<co-col; j++)
     188             :   {
     189     1205480 :     GEN z = gel(B,j-s);
     190     1205480 :     p1 = cgetg(lig+1,t_COL); gel(Bnew,j) = p1;
     191     1205480 :     for (i=1; i<=nlze; i++) gel(p1,i) = gel(z,i);
     192     1205480 :     z += nlze; p1 += nlze;
     193     5523948 :     for (i=k=1; k<=lnz; i++)
     194     4318468 :       if (!diagH1[i]) gel(p1,k++) = gel(z,i);
     195             :   }
     196       59214 :   *ptdep = depnew;
     197       59214 :   *ptC = C;
     198       59214 :   *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      816452 : col_dup(long l, GEN col)
     215             : {
     216      816452 :   GEN c = new_chunk(l);
     217      816452 :   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       43033 : 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       43033 :   const long li = lg(perm); /* = lgcols(mat0) */
     241       43033 :   const long CO = lg(mat0);
     242             : 
     243       43033 :   n = 0; /* -Wall */
     244             : 
     245       43033 :   C = *ptC; co = CO;
     246       43033 :   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       43033 :   if (DEBUGLEVEL>5)
     253             :   {
     254           0 :     err_printf("Entering hnfspec\n");
     255           0 :     p_mat(mat0,perm,0);
     256             :   }
     257       43033 :   matt = cgetg(co, t_MAT); /* dense part of mat (top) */
     258       43033 :   mat = cgetg(co, t_MAT);
     259      859485 :   for (j = 1; j < co; j++)
     260             :   {
     261      816452 :     GEN matj = col_dup(li, gel(mat0,j));
     262      816452 :     p1 = cgetg(k0+1,t_COL); gel(matt,j) = p1; gel(mat,j) = matj;
     263      816452 :     for (i=1; i<=k0; i++) gel(p1,i) = stoi(matj[perm[i]]);
     264             :   }
     265       43033 :   av = avma;
     266             : 
     267       43033 :   i = lig = li-1; col = co-1; lk0 = k0;
     268       43033 :   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     2385911 :   while (i > lk0 && col)
     271     2299845 :     switch( count(mat,perm[i],col,&n) )
     272             :     {
     273             :       case 0: /* move zero lines between k0+1 and lk0 */
     274       17138 :         lk0++; lswap(perm[i], perm[lk0]);
     275       17138 :         i = lig; continue;
     276             : 
     277             :       case 1: /* move trivial generator between lig+1 and li */
     278      142974 :         lswap(perm[i], perm[lig]);
     279      142974 :         if (T) swap(gel(T,n), gel(T,col));
     280      142974 :         swap(gel(mat,n), gel(mat,col)); p = gel(mat,col);
     281      142974 :         if (p[perm[lig]] < 0) /* = -1 */
     282             :         { /* convert relation -g = 0 to g = 0 */
     283      127725 :           for (i=lk0+1; i<lig; i++) p[perm[i]] = -p[perm[i]];
     284      127725 :           if (T)
     285             :           {
     286      127725 :             p1 = gel(T,col);
     287     5962673 :             for (i=1; ; i++) /* T is a permuted identity: single non-0 entry */
     288     5962673 :               if (signe(gel(p1,i))) { togglesign_safe(&gel(p1,i)); break; }
     289     5834948 :           }
     290             :         }
     291      142974 :         lig--; col--; i = lig; continue;
     292             : 
     293     2139733 :       default: i--;
     294             :     }
     295       43033 :   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       43033 :   s = 0;
     302      314341 :   while (lig > lk0 && col && s < (long)(HIGHBIT>>1))
     303             :   {
     304      788768 :     for (i=lig; i>lk0; i--)
     305      754161 :       if (count(mat,perm[i],col,&n) > 0) break;
     306      262882 :     if (i == lk0) break;
     307             : 
     308             :     /* only 0, +/- 1 entries, at least 2 of them non-zero */
     309      228275 :     lswap(perm[i], perm[lig]);
     310      228275 :     swap(gel(mat,n), gel(mat,col)); p = gel(mat,col);
     311      228275 :     if (T) swap(gel(T,n), gel(T,col));
     312      228275 :     if (p[perm[lig]] < 0)
     313             :     {
     314      171737 :       for (i=lk0+1; i<=lig; i++) p[perm[i]] = -p[perm[i]];
     315      171737 :       if (T) ZV_togglesign(gel(T,col));
     316             :     }
     317     9352035 :     for (j=1; j<col; j++)
     318             :     {
     319     9123760 :       GEN matj = gel(mat,j);
     320             :       long t;
     321     9123760 :       if (! (t = matj[perm[lig]]) ) continue;
     322      592673 :       if (t == 1) {
     323      267557 :         for (i=lk0+1; i<=lig; i++) absmax(s, matj[perm[i]] -= p[perm[i]]);
     324             :       }
     325             :       else { /* t = -1 */
     326      325116 :         for (i=lk0+1; i<=lig; i++) absmax(s, matj[perm[i]] += p[perm[i]]);
     327             :       }
     328      592673 :       if (T) ZC_lincomb1_inplace(gel(T,j), gel(T,col), stoi(-t));
     329             :     }
     330      228275 :     lig--; col--;
     331      228275 :     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       43033 :   vmax = cgetg(co,t_VECSMALL);
     340      488236 :   for (j=1; j<=col; j++)
     341             :   {
     342      445203 :     GEN matj = gel(mat,j);
     343      445203 :     for (s=0, i=lk0+1; i<=lig; i++) absmax(s, matj[i]);
     344      445203 :     vmax[j] = s;
     345             :   }
     346      174616 :   while (lig > lk0 && col)
     347             :   {
     348      129067 :     for (i=lig; i>lk0; i--)
     349      122419 :       if ( (n = count2(mat,perm[i],col)) ) break;
     350       95198 :     if (i == lk0) break;
     351             : 
     352       88550 :     lswap(vmax[n], vmax[col]);
     353       88550 :     lswap(perm[i], perm[lig]);
     354       88550 :     swap(gel(mat,n), gel(mat,col)); p = gel(mat,col);
     355       88550 :     if (T) swap(gel(T,n), gel(T,col));
     356       88550 :     if (p[perm[lig]] < 0)
     357             :     {
     358       54085 :       for (i=lk0+1; i<=lig; i++) p[perm[i]] = -p[perm[i]];
     359       54085 :       if (T) ZV_togglesign(gel(T,col));
     360             :     }
     361     1319630 :     for (j=1; j<col; j++)
     362             :     {
     363     1231080 :       GEN matj = gel(mat,j);
     364             :       long t;
     365     1231080 :       if (! (t = matj[perm[lig]]) ) continue;
     366      621186 :       if (vmax[col] && (ulong)labs(t) >= (HIGHBIT-vmax[j]) / vmax[col])
     367           0 :         goto END2;
     368             : 
     369      621186 :       for (s=0, i=lk0+1; i<=lig; i++) absmax(s, matj[perm[i]] -= t*p[perm[i]]);
     370      621186 :       vmax[j] = s;
     371      621186 :       if (T) ZC_lincomb1_inplace(gel(T,j), gel(T,col), stoi(-t));
     372             :     }
     373       88550 :     lig--; col--;
     374       88550 :     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       43033 :   matb = cgetg(co,t_MAT); /* bottom part (complement of matt) */
     384      859485 :   for (j=1; j<co; j++)
     385             :   {
     386      816452 :     GEN matj = gel(mat,j);
     387      816452 :     p1 = cgetg(li-k0,t_COL); gel(matb,j) = p1;
     388      816452 :     p1 -= k0;
     389      816452 :     for (i=k0+1; i<li; i++) gel(p1,i) = stoi(matj[perm[i]]);
     390             :   }
     391       43033 :   if (DEBUGLEVEL>5)
     392             :   {
     393           0 :     err_printf("    after phase2:\n");
     394           0 :     p_mat(mat,perm,lk0);
     395             :   }
     396      460169 :   for (i=li-2; i>lig; i--)
     397             :   {
     398      417136 :     long h, i0 = i - k0, k = i + co-li;
     399      417136 :     GEN Bk = gel(matb,k);
     400    13065087 :     for (j=k+1; j<co; j++)
     401             :     {
     402    12647951 :       GEN Bj = gel(matb,j), v = gel(Bj,i0);
     403    12647951 :       s = signe(v); if (!s) continue;
     404             : 
     405     2595435 :       gel(Bj,i0) = gen_0;
     406     2595435 :       if (is_pm1(v))
     407             :       {
     408     1475796 :         if (s > 0) /* v = 1 */
     409      725112 :         { for (h=1; h<i0; h++) gel(Bj,h) = subii(gel(Bj,h), gel(Bk,h)); }
     410             :         else /* v = -1 */
     411      750684 :         { for (h=1; h<i0; h++) gel(Bj,h) = addii(gel(Bj,h), gel(Bk,h)); }
     412             :       }
     413             :       else {
     414     1119639 :         for (h=1; h<i0; h++) gel(Bj,h) = subii(gel(Bj,h), mulii(v,gel(Bk,h)));
     415             :       }
     416     2595435 :       if (T) ZC_lincomb1_inplace(gel(T,j), gel(T,k), negi(v));
     417     2595435 :       if (gc_needed(av,3))
     418             :       {
     419          26 :         if(DEBUGMEM>1) pari_warn(warnmem,"hnfspec[3], (i,j) = %ld,%ld", i,j);
     420          26 :         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       43033 :   for (j=1; j<co; j++) setlg(matb[j], lig-k0+1); /* bottom can be forgotten */
     427       43033 :   gerepileall(av, T? 2: 1, &matb, &T);
     428       43033 :   if (DEBUGLEVEL>5) err_printf("    matb cleaned up (using Id block)\n");
     429             : 
     430       43033 :   nlze = lk0 - k0;  /* # of 0 rows */
     431       43033 :   lnz = lig-nlze+1; /* 1 + # of non-0 rows (!= 0...0 1 0 ... 0) */
     432       43033 :   if (T) matt = ZM_mul(matt,T); /* update top rows */
     433       43033 :   extramat = cgetg(col+1,t_MAT); /* = new C minus the 0 rows */
     434      399686 :   for (j=1; j<=col; j++)
     435             :   {
     436      356653 :     GEN z = gel(matt,j);
     437      356653 :     GEN t = (gel(matb,j)) + nlze - k0;
     438      356653 :     p2=cgetg(lnz,t_COL); gel(extramat,j) = p2;
     439      356653 :     for (i=1; i<=k0; i++) gel(p2,i) = gel(z,i); /* top k0 rows */
     440      356653 :     for (   ; i<lnz; i++) gel(p2,i) = gel(t,i); /* other non-0 rows */
     441             :   }
     442       43033 :   if (!col) {
     443         154 :     permpro = identity_perm(lnz);
     444         154 :     nr = lnz;
     445             :   }
     446             :   else
     447       42879 :     permpro = ZM_imagecomplspec(extramat, &nr);
     448             :   /* lnz = lg(permpro) */
     449       43033 :   if (nlze)
     450             :   { /* put the nlze 0 rows (trivial generators) at the top */
     451        3979 :     p1 = new_chunk(lk0+1);
     452        3979 :     for (i=1; i<=nlze; i++) p1[i] = perm[i + k0];
     453        3979 :     for (   ; i<=lk0; i++)  p1[i] = perm[i - nlze];
     454        3979 :     for (i=1; i<=lk0; i++)  perm[i] = p1[i];
     455             :   }
     456             :   /* sort other rows according to permpro (nr redundant generators first) */
     457       43033 :   p1 = new_chunk(lnz); p2 = perm + nlze;
     458       43033 :   for (i=1; i<lnz; i++) p1[i] = p2[permpro[i]];
     459       43033 :   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       43033 :   matbnew = cgetg(col+1,t_MAT); /* dense+toobig, maximal rank. For hnffinal */
     475       43033 :   dep    = cgetg(col+1,t_MAT); /* rows dependent from the ones in matbnew */
     476      399686 :   for (j=1; j<=col; j++)
     477             :   {
     478      356653 :     GEN z = gel(extramat,j);
     479      356653 :     p1 = cgetg(nlze+nr+1,t_COL); gel(dep,j) = p1;
     480      356653 :     p2 = cgetg(lnz-nr,t_COL); gel(matbnew,j) = p2;
     481      356653 :     for (i=1; i<=nlze; i++) gel(p1,i) = gen_0;
     482      356653 :     p1 += nlze; for (i=1; i<=nr; i++) p1[i] = z[permpro[i]];
     483      356653 :     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       43033 :   B = cgetg(co-col,t_MAT);
     489      502832 :   for (j=col+1; j<co; j++)
     490             :   {
     491      459799 :     GEN y = gel(matt,j);
     492      459799 :     GEN z = gel(matb,j);
     493      459799 :     p1=cgetg(lig+1,t_COL); gel(B,j-col) = p1;
     494      459799 :     for (i=1; i<=nlze; i++) gel(p1,i) = gel(z,i);
     495      459799 :     p1 += nlze; z += nlze-k0;
     496     3718017 :     for (k=1; k<lnz; k++)
     497             :     {
     498     3258218 :       i = permpro[k];
     499     3258218 :       gel(p1,k) = (i <= k0)? gel(y,i): gel(z,i);
     500             :     }
     501             :   }
     502       43033 :   if (T) C = typ(C)==t_MAT? RgM_mul(C,T): RgV_RgM_mul(C,T);
     503       43033 :   gerepileall(av, 4, &matbnew, &B, &dep, &C);
     504       43033 :   *ptdep = dep;
     505       43033 :   *ptB = B;
     506       43033 :   H = hnffinal(matbnew, perm, ptdep, ptB, &C);
     507       43033 :   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       43033 :   *ptC = C; return H;
     534             : }
     535             : 
     536             : GEN
     537         147 : hnfspec(GEN mat, GEN perm, GEN* ptdep, GEN* ptB, GEN* ptC, long k0)
     538             : {
     539         147 :   pari_sp av = avma;
     540         147 :   GEN H = hnfspec_i(mat, perm, ptdep, ptB, ptC, k0);
     541         147 :   gerepileall(av, 4, ptC, ptdep, ptB, &H); return H;
     542             : }
     543             : 
     544             : /* HNF reduce x, apply same transforms to C */
     545             : GEN
     546         147 : mathnfspec(GEN x, GEN *ptperm, GEN *ptdep, GEN *ptB, GEN *ptC)
     547             : {
     548         147 :   long i,j,k,ly,lx = lg(x);
     549             :   GEN z, perm;
     550         147 :   if (lx == 1) return cgetg(1, t_MAT);
     551         147 :   ly = lgcols(x);
     552         147 :   *ptperm = perm = identity_perm(ly-1);
     553         147 :   z = cgetg(lx,t_MAT);
     554        1967 :   for (i=1; i<lx; i++)
     555             :   {
     556        1820 :     GEN C = cgetg(ly,t_COL), D = gel(x,i);
     557        1820 :     gel(z,i) = C;
     558       85652 :     for (j=1; j<ly; j++)
     559             :     {
     560       83832 :       GEN d = gel(D,j);
     561       83832 :       if (is_bigint(d)) goto TOOLARGE;
     562       83832 :       C[j] = itos(d);
     563             :     }
     564             :   }
     565             :   /*  [ dep |     ]
     566             :    *  [-----|  B  ]
     567             :    *  [  H  |     ]
     568             :    *  [-----|-----]
     569             :    *  [  0  | Id  ] */
     570         147 :   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           0 :   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       16335 : hnfadd_i(GEN H, GEN perm, GEN* ptdep, GEN* ptB, GEN* ptC, /* cf hnfspec */
     595             :        GEN extramat,GEN extraC)
     596             : {
     597       16335 :   GEN matb, extratop, Cnew, permpro, B = *ptB, C = *ptC, dep = *ptdep;
     598             :   long i, lH, lB, li, lig, co, col, nlze;
     599             : 
     600       16335 :   if (lg(extramat) == 1) return H;
     601       16335 :   co = lg(C)-1;
     602       16335 :   lH = lg(H)-1;
     603       16335 :   lB = lg(B)-1;
     604       16335 :   li = lg(perm)-1;
     605       16335 :   lig = li - lB;
     606       16335 :   col = co - lB;
     607       16335 :   nlze = lig - lH;
     608             : 
     609             :  /*               col    co
     610             :   *       [ 0 |dep |     ]
     611             :   *  nlze [--------|  B  ]
     612             :   *       [ 0 | H  |     ]
     613             :   *       [--------|-----] lig
     614             :   *       [   0    | Id  ]
     615             :   *       [        |     ] li */
     616       16335 :   extratop = zm_to_ZM( rowslicepermute(extramat, perm, 1, lig) );
     617       16335 :   if (li != lig)
     618             :   { /* zero out bottom part, using the Id block */
     619       16335 :     GEN A = vecslice(C, col+1, co);
     620       16335 :     GEN c = rowslicepermute(extramat, perm, lig+1, li);
     621       16335 :     extraC   = gsub(extraC, typ(A)==t_MAT? RgM_zm_mul(A, c): RgV_zm_mul(A,c));
     622       16335 :     extratop = ZM_sub(extratop, ZM_zm_mul(B, c));
     623             :   }
     624             : 
     625       16335 :   extramat = shallowconcat(extratop, vconcat(dep, H));
     626       16335 :   Cnew     = shallowconcat(extraC, vecslice(C, col-lH+1, co));
     627       16335 :   if (DEBUGLEVEL>5) err_printf("    1st phase done\n");
     628       16335 :   permpro = ZM_imagecomplspec(extramat, &nlze);
     629       16335 :   extramat = rowpermute(extramat, permpro);
     630       16335 :   *ptB     = rowpermute(B,        permpro);
     631       16335 :   permpro = vecsmallpermute(perm, permpro);
     632       16335 :   for (i=1; i<=lig; i++) perm[i] = permpro[i]; /* perm o= permpro */
     633             : 
     634       16335 :   *ptdep  = rowslice(extramat, 1, nlze);
     635       16335 :   matb    = rowslice(extramat, nlze+1, lig);
     636       16335 :   if (DEBUGLEVEL>5) err_printf("    2nd phase done\n");
     637       16335 :   H = hnffinal(matb,perm,ptdep,ptB,&Cnew);
     638       16335 :   *ptC = shallowconcat(vecslice(C, 1, col-lH), Cnew);
     639       16335 :   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    27239868 : ZC_elem(GEN aj, GEN ak, GEN A, GEN U, long j, long k)
     655             : {
     656             :   GEN p1,u,v,d;
     657             : 
     658    27239868 :   if (!signe(ak)) {
     659     5099516 :     swap(gel(A,j), gel(A,k));
     660     5099516 :     if (U) swap(gel(U,j), gel(U,k));
     661    29666612 :     return;
     662             :   }
     663    22140352 :   d = bezout(aj,ak,&u,&v);
     664             :   /* frequent special case (u,v) = (1,0) or (0,1) */
     665    22140352 :   if (!signe(u))
     666             :   { /* ak | aj */
     667    15059238 :     p1 = diviiexact(aj,ak); togglesign(p1);
     668    15059238 :     ZC_lincomb1_inplace(gel(A,j), gel(A,k), p1);
     669    15059238 :     if (U)
     670     1091976 :       ZC_lincomb1_inplace(gel(U,j), gel(U,k), p1);
     671    15059238 :     return;
     672             :   }
     673     7081114 :   if (!signe(v))
     674             :   { /* aj | ak */
     675     4408342 :     p1 = diviiexact(ak,aj); togglesign(p1);
     676     4408342 :     ZC_lincomb1_inplace(gel(A,k), gel(A,j), p1);
     677     4408342 :     swap(gel(A,j), gel(A,k));
     678     4408342 :     if (U) {
     679      101376 :       ZC_lincomb1_inplace(gel(U,k), gel(U,j), p1);
     680      101376 :       swap(gel(U,j), gel(U,k));
     681             :     }
     682     4408342 :     return;
     683             :   }
     684             : 
     685     2672772 :   if (!is_pm1(d)) { aj = diviiexact(aj, d); ak = diviiexact(ak, d); }
     686     2672772 :   p1 = gel(A,k); aj = negi(aj); /* NOT togglesign */
     687     2672772 :   gel(A,k) = ZC_lincomb(u,v, gel(A,j),p1);
     688     2672772 :   gel(A,j) = ZC_lincomb(aj,ak, p1,gel(A,j));
     689     2672772 :   if (U)
     690             :   {
     691      277039 :     p1 = gel(U,k);
     692      277039 :     gel(U,k) = ZC_lincomb(u,v, gel(U,j),p1);
     693      277039 :     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_coeff(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      268856 : col_mul(GEN x, GEN c)
     736             : {
     737      268856 :   if (typ(x) == t_INT)
     738             :   {
     739      268856 :     long s = signe(x);
     740      268856 :     if (!s) return NULL;
     741      201761 :     if (is_pm1(x)) return (s > 0)? c: RgC_neg(c);
     742             :   }
     743       38997 :   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       67214 : update(GEN u, GEN v, GEN a, GEN b, GEN *c1, GEN *c2)
     755             : {
     756             :   GEN p1,p2;
     757             : 
     758       67214 :   u = col_mul(u,*c1);
     759       67214 :   v = col_mul(v,*c2);
     760       67214 :   if (u) p1 = v? gadd(u,v): u;
     761         161 :   else   p1 = v? v: NULL;
     762             : 
     763       67214 :   a = col_mul(a,*c2);
     764       67214 :   b = col_mul(gneg_i(b),*c1);
     765       67214 :   if (a) p2 = b? RgC_add(a,b): a;
     766           0 :   else   p2 = b? b: NULL;
     767             : 
     768       67214 :   if (!p1) do_zero(*c1); else *c1 = p1;
     769       67214 :   if (!p2) do_zero(*c2); else *c2 = p2;
     770       67214 : }
     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     1230975 : ZM_reduce(GEN A, GEN U, long i, long j0)
     794             : {
     795     1230975 :   long j, lA = lg(A);
     796     1230975 :   GEN d = gcoeff(A,i,j0);
     797     1230975 :   if (signe(d) < 0)
     798             :   {
     799        3914 :     ZV_neg_inplace(gel(A,j0));
     800        3914 :     if (U) ZV_togglesign(gel(U,j0));
     801        3914 :     d = gcoeff(A,i,j0);
     802             :   }
     803     3342916 :   for (j=j0+1; j<lA; j++)
     804             :   {
     805     2111941 :     GEN q = truedivii(gcoeff(A,i,j), d);
     806     2111941 :     if (!signe(q)) continue;
     807             : 
     808      107177 :     togglesign(q);
     809      107177 :     ZC_lincomb1_inplace(gel(A,j), gel(A,j0), q);
     810      107177 :     if (U) ZC_lincomb1_inplace(gel(U,j), gel(U,j0), q);
     811             :   }
     812     1230975 : }
     813             : 
     814             : static GEN
     815         168 : RgX_normalize_all(GEN T, GEN *pd)
     816             : {
     817         168 :   GEN d = leading_coeff(T);
     818         336 :   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_coeff(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      319189 : hnfmerge_get_1(GEN A, GEN B)
     862             : {
     863      319189 :   pari_sp av = avma;
     864      319189 :   long j, k, c, l = lg(A), lb;
     865      319189 :   GEN b, t, U = cgetg(l + 1, t_MAT), C = cgetg(l + 1, t_VEC);
     866             : 
     867      319189 :   t = NULL; /* -Wall */
     868      319189 :   b = gcoeff(B,1,1); lb = lgefint(b);
     869      319189 :   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      755691 :   for (j = 1; j < l; j++)
     874             :   {
     875      755684 :     c = j+1;
     876      755684 :     gel(U,j) = col_ei(l-1, j);
     877      755684 :     gel(U,c) = zerocol(l-1); /* dummy */
     878      755684 :     gel(C,j) = vecslice(gel(A,j), 1,j);
     879      755684 :     gel(C,c) = vecslice(gel(B,j), 1,j);
     880     2155964 :     for (k = j; k > 0; k--)
     881             :     {
     882     1400280 :       t = gcoeff(C,k,c);
     883     1400280 :       if (gequal0(t)) continue;
     884     1285820 :       setlg(C[c], k+1);
     885     1285820 :       ZC_elem(t, gcoeff(C,k,k), C, U, c, k);
     886     1285820 :       if (lgefint(gcoeff(C,k,k)) > lb) gel(C,k) = FpC_red(gel(C,k), b);
     887     1285820 :       if (j > 4)
     888             :       {
     889       57358 :         GEN u = gel(U,k);
     890             :         long h;
     891      695058 :         for (h=1; h<l; h++)
     892      637700 :           if (lgefint(gel(u,h)) > lb) gel(u,h) = remii(gel(u,h), b);
     893             :       }
     894             :     }
     895      755684 :     if (j == 1)
     896      319189 :       t = gcoeff(C,1,1);
     897             :     else
     898             :     {
     899             :       GEN u;
     900      436495 :       t = bezout(gcoeff(C,1,1), b, &u, NULL); /* >= 0 */
     901      436495 :       if (signe(u) && !equali1(u)) gel(U,1) = ZC_Z_mul(gel(U,1), u);
     902      436495 :       gcoeff(C,1,1) = t;
     903             :     }
     904      755684 :     if (equali1(t)) break;
     905             :   }
     906      319189 :   if (j >= l) return NULL;
     907      319182 :   return gerepileupto(av, ZM_ZC_mul(A,gel(U,1)));
     908             : }
     909             : 
     910             : /* remove the first r columns */
     911             : static void
     912       19354 : remove_0cols(long r, GEN *pA, GEN *pB, long remove)
     913             : {
     914       19354 :   GEN A = *pA, B = *pB;
     915       19354 :   long l = lg(A);
     916       19354 :   A += r; A[0] = evaltyp(t_MAT) | evallg(l-r);
     917       19354 :   if (B && remove == 2) { B += r; B[0] = A[0]; }
     918       19354 :   *pA = A; *pB = B;
     919       19354 : }
     920             : 
     921             : /* Inefficient compared to hnfall. 'remove' = throw away lin.dep columns */
     922             : static GEN
     923        8362 : hnf_i(GEN A, int remove)
     924             : {
     925        8362 :   pari_sp av0 = avma, av;
     926             :   long s, n, m, j, k, li, def, ldef;
     927             : 
     928        8362 :   RgM_dimensions(A, &m, &n);
     929        8362 :   if (!n) return cgetg(1,t_MAT);
     930        8285 :   av = avma;
     931        8285 :   A = RgM_shallowcopy(A);
     932        8285 :   def = n; ldef = (m>n)? m-n: 0;
     933       29282 :   for (li=m; li>ldef; li--)
     934             :   {
     935       52070 :     for (j=def-1; j; j--)
     936             :     {
     937       31073 :       GEN a = gcoeff(A,li,j);
     938       31073 :       if (!signe(a)) continue;
     939             : 
     940             :       /* zero a = Aij  using  b = Aik */
     941       14200 :       k = (j==1)? def: j-1;
     942       14200 :       ZC_elem(a,gcoeff(A,li,k), A,NULL, j,k);
     943       14200 :       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       20997 :     s = signe(gcoeff(A,li,def));
     950       20997 :     if (s)
     951             :     {
     952       20969 :       if (s < 0) ZV_neg_inplace(gel(A,def));
     953       20969 :       ZM_reduce(A, NULL, li,def);
     954       20969 :       def--;
     955             :     }
     956             :     else
     957          28 :       if (ldef) ldef--;
     958       20997 :     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        8285 :   if (remove) { GEN B = NULL; remove_0cols(def, &A, &B, remove); }
     966        8285 :   return gerepileupto(av0, ZM_copy(A));
     967             : }
     968             : 
     969             : GEN
     970       10165 : 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     1065352 : FpV_Fp_mul_part_ip(GEN z, GEN u, GEN p, long k)
     975             : {
     976             :   long i;
     977     1065352 :   if (is_pm1(u)) {
     978      554951 :     if (signe(u) > 0) {
     979     1756768 :       for (i = 1; i <= k; i++)
     980     1246663 :         if (signe(gel(z,i))) gel(z,i) = modii(gel(z,i), p);
     981             :     } else {
     982      133217 :       for (i = 1; i <= k; i++)
     983       88371 :         if (signe(gel(z,i))) gel(z,i) = modii(negi(gel(z,i)), p);
     984             :     }
     985             :   }
     986             :   else {
     987     2619391 :     for (i = 1; i <= k; i++)
     988     2108990 :       if (signe(gel(z,i))) gel(z,i) = Fp_mul(u,gel(z,i), p);
     989             :   }
     990     1065352 : }
     991             : static void
     992    19438986 : FpV_red_part_ipvec(GEN z, GEN p, long k)
     993             : {
     994             :   long i;
     995    19438986 :   for (i = 1; i <= k; i++) gel(z,i) = modii(gel(z,i), gel(p,i));
     996    19438986 : }
     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       45364 : ZpM_echelon(GEN x, long early_abort, GEN p, GEN pm)
    1002             : {
    1003       45364 :   pari_sp av0 = avma, av;
    1004             :   long m, li, co, i, j, k, def, ldef;
    1005             : 
    1006       45364 :   co = lg(x); if (co == 1) return cgetg(1,t_MAT);
    1007       45364 :   li = lgcols(x);
    1008       45364 :   av = avma;
    1009       45364 :   x = RgM_shallowcopy(x);
    1010       45364 :   m = Z_pval(pm, p);
    1011             : 
    1012       45364 :   ldef = (li > co)? li - co: 0;
    1013      413287 :   for (def = co-1,i = li-1; i > ldef; i--)
    1014             :   {
    1015      368347 :     long vmin = LONG_MAX, kmin = 0;
    1016      368347 :     GEN umin = gen_0, pvmin, q;
    1017     3403508 :     for (k = 1; k <= def; k++)
    1018             :     {
    1019     3104688 :       GEN u = gcoeff(x,i,k);
    1020             :       long v;
    1021     3104688 :       if (!signe(u)) continue;
    1022     1491855 :       v = Z_pvalrem(u, p, &u);
    1023     1491855 :       if (v >= m) gcoeff(x,i,k) = gen_0;
    1024      964042 :       else if (v < vmin) {
    1025      357139 :         vmin = v; kmin = k; umin = u;
    1026      357139 :         if (!vmin) break;
    1027             :       }
    1028             :     }
    1029      368347 :     if (!kmin)
    1030             :     {
    1031      151339 :       if (early_abort) return NULL;
    1032      150915 :       gcoeff(x,i,def) = gen_0;
    1033      150915 :       ldef--;
    1034      150915 :       if (ldef < 0) ldef = 0;
    1035      150915 :       continue;
    1036             :     }
    1037      217008 :     if (kmin != def) swap(gel(x,def), gel(x,kmin));
    1038      217008 :     q = vmin? powiu(p, m-vmin): pm;
    1039             :     /* pivot has valuation vmin */
    1040      217008 :     umin = modii(umin, q);
    1041      217008 :     if (!equali1(umin))
    1042      172155 :       FpV_Fp_mul_part_ip(gel(x,def), Fp_inv(umin,q), pm, i-1);
    1043      217008 :     gcoeff(x, i, def) = pvmin = powiu(p, vmin);
    1044     2560376 :     for (j = def-1; j; j--)
    1045             :     { /* zero x[i, 1..def-1] using x[i,def] = pvmin */
    1046     2343368 :       GEN t, a = gcoeff(x,i,j) = modii(gcoeff(x,i,j), pm);
    1047     2343368 :       if (!signe(a)) continue;
    1048             : 
    1049     1143085 :       t = diviiexact(a, pvmin); togglesign(t);
    1050     1143085 :       ZC_lincomb1_inplace(gel(x,j), gel(x,def), t);
    1051     1143085 :       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      217008 :     def--;
    1058             :   }
    1059       44940 :   if (co > li)
    1060             :   {
    1061           0 :     x += co - li;
    1062           0 :     x[0] = evaltyp(t_MAT) | evallg(li);
    1063             :   }
    1064       44940 :   return gerepilecopy(av0, x);
    1065             : }
    1066             : GEN
    1067      186049 : zlm_echelon(GEN x, long early_abort, ulong p, ulong pm)
    1068             : {
    1069      186049 :   pari_sp av0 = avma;
    1070             :   long li, co, i, j, k, def, ldef;
    1071             :   ulong m;
    1072             : 
    1073      186049 :   co = lg(x); if (co == 1) return cgetg(1,t_MAT);
    1074      186049 :   li = lgcols(x);
    1075      186049 :   x = Flm_copy(x);
    1076      186049 :   m = u_lval(pm, p);
    1077             : 
    1078      186049 :   ldef = (li > co)? li - co: 0;
    1079     1584492 :   for (def = co-1,i = li-1; i > ldef; i--)
    1080             :   {
    1081     1424759 :     long vmin = LONG_MAX, kmin = 0;
    1082     1424759 :     ulong umin = 0, pvmin, q;
    1083     8083625 :     for (k = 1; k <= def; k++)
    1084             :     {
    1085     7347082 :       ulong u = ucoeff(x,i,k);
    1086             :       long v;
    1087     7347082 :       if (!u) continue;
    1088     3342834 :       v = u_lvalrem(u, p, &u);
    1089     3342834 :       if (v >= (long) m) ucoeff(x,i,k) = 0;
    1090     3342834 :       else if (v < vmin) {
    1091     1456001 :         vmin = v; kmin = k; umin = u;
    1092     1456001 :         if (!vmin) break;
    1093             :       }
    1094             :     }
    1095     1424759 :     if (!kmin)
    1096             :     {
    1097      203120 :       if (early_abort) return NULL;
    1098      176804 :       ucoeff(x,i,def) = 0;
    1099      176804 :       ldef--;
    1100      176804 :       if (ldef < 0) ldef = 0;
    1101      176804 :       continue;
    1102             :     }
    1103     1221639 :     if (kmin != def) swap(gel(x,def), gel(x,kmin));
    1104     1221639 :     q = vmin? upowuu(p, m-vmin): pm;
    1105             :     /* pivot has valuation vmin */
    1106     1221639 :     umin %= q;
    1107     1221639 :     if (umin != 1)
    1108      993870 :       Flv_Fl_mul_part_inplace(gel(x,def), Fl_inv(umin,q), pm, i-1);
    1109     1221639 :     ucoeff(x, i, def) = pvmin = upowuu(p, vmin);
    1110    13390118 :     for (j = def-1; j; j--)
    1111             :     { /* zero x[i, 1..def-1] using x[i,def] = pvmin */
    1112    12168479 :       ulong t, a = ucoeff(x,i,j);
    1113    12168479 :       if (!a) continue;
    1114             : 
    1115     7183793 :       t = Fl_neg(a / pvmin, q);
    1116     7183793 :       Flc_lincomb1_inplace(gel(x,j), gel(x,def), t, pm);
    1117             :     }
    1118     1221639 :     def--;
    1119             :   }
    1120      159733 :   if (co > li)
    1121             :   {
    1122           0 :     x += co - li;
    1123           0 :     x[0] = evaltyp(t_MAT) | evallg(li);
    1124             :   }
    1125      159733 :   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     3422211 : ZM_hnfmodall_i(GEN x, GEN dm, long flag)
    1135             : {
    1136             :   pari_sp av;
    1137     3422211 :   const long center = (flag & hnf_CENTER);
    1138     3422211 :   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     3422211 :   co = lg(x); if (co == 1) return cgetg(1,t_MAT);
    1143     3422211 :   li = lgcols(x); if (li == 1) return cgetg(1,t_MAT);
    1144     3422148 :   if (typ(dm) == t_INT)
    1145             :   {
    1146     3409107 :     long ldm = lgefint(dm);
    1147     3409107 :     dm2 = shifti(dm, -1);
    1148     3409107 :     dm = const_vec(li-1,dm);
    1149     3409107 :     dm2= const_vec(li-1,dm2);
    1150     3409107 :     LDM= const_vecsmall(li-1,ldm);
    1151             :   }
    1152             :   else
    1153             :   {
    1154       13041 :     if (lg(dm) != li) pari_err_DIM("ZM_hnfmod");
    1155       13041 :     moddiag = 1;
    1156       13041 :     dm2 = cgetg(li, t_VEC);
    1157       13041 :     LDM = cgetg(li, t_VECSMALL);
    1158       39837 :     for (i=1; i<li; i++)
    1159             :     {
    1160       26796 :       gel(dm2,i) = shifti(gel(dm,i),-1);
    1161       26796 :       LDM[i] = lgefint(gel(dm,i));
    1162             :     }
    1163             :   }
    1164     3422148 :   av = avma;
    1165     3422148 :   x = RgM_shallowcopy(x);
    1166             : 
    1167     3422148 :   ldef = 0;
    1168     3422148 :   if (li > co)
    1169             :   {
    1170        5061 :     ldef = li - co;
    1171        5061 :     if (!moddiag)
    1172           7 :       pari_err_DOMAIN("ZM_hnfmod","nb lines",">", strtoGENstr("nb columns"), x);
    1173             :   }
    1174    13979977 :   for (def = co-1,i = li-1; i > ldef; i--,def--)
    1175             :   {
    1176    10557836 :     GEN d = gel(dm,i), d2 = gel(dm2,i);
    1177    10557836 :     gcoeff(x,i,def) = centermodii(gcoeff(x,i,def), d,d2);
    1178    40328455 :     for (j = def-1; j; j--)
    1179             :     {
    1180    29770619 :       gcoeff(x,i,j) = centermodii(gcoeff(x,i,j), d,d2);
    1181    29770619 :       a = gcoeff(x,i,j);
    1182    29770619 :       if (!signe(a)) continue;
    1183             : 
    1184    15886651 :       k = (j==1)? def: j-1;
    1185    15886651 :       gcoeff(x,i,k) = centermodii(gcoeff(x,i,k), d,d2);
    1186    15886651 :       ZC_elem(a,gcoeff(x,i,k), x,NULL, j,k);
    1187    15886651 :       p1 = gel(x,j);
    1188    15886651 :       p2 = gel(x,k);
    1189             :       /* prevent coeffs explosion: reduce mod dm when lg() > ldm */
    1190    96999770 :       for (k = 1; k < i; k++)
    1191             :       {
    1192    81113119 :         if (lgefint(gel(p1,k)) > LDM[k])
    1193     7210192 :           gel(p1,k) = centermodii(gel(p1,k), gel(dm,k),gel(dm2,k));
    1194    81113119 :         if (lgefint(gel(p2,k)) > LDM[k])
    1195     1237288 :           gel(p2,k) = centermodii(gel(p2,k), gel(dm,k),gel(dm2,k));
    1196             :       }
    1197             :     }
    1198    10557836 :     if (gc_needed(av,1))
    1199             :     {
    1200          59 :       if (DEBUGMEM>1) pari_warn(warnmem,"ZM_hnfmod[1]. i=%ld",i);
    1201          59 :       x = gerepilecopy(av, x);
    1202             :     }
    1203    10557836 :     if (moddiag && !signe(gcoeff(x,i,def)))
    1204             :     { /* missing pivot on line i, insert column */
    1205     3845487 :       GEN a = cgetg(co + 1, t_MAT);
    1206     3845487 :       for (k = 1; k <= def; k++) gel(a,k) = gel(x,k);
    1207     3845487 :       gel(a,k++) = Rg_col_ei(gel(dm,i), li-1, i);
    1208     3845487 :       for (     ; k <= co;  k++) gel(a,k) = gel(x,k-1);
    1209     3845487 :       ldef--; if (ldef < 0) ldef = 0;
    1210     3845487 :       co++; def++; x = a;
    1211             :     }
    1212             :   }
    1213     3422141 :   if (co < li)
    1214             :   { /* implies moddiag, add missing diag(dm) components */
    1215        4235 :     GEN a = cgetg(li+1, t_MAT);
    1216        4235 :     for (k = 1; k <= li-co; k++) gel(a,k) = Rg_col_ei(gel(dm,k), li-1, k);
    1217        4235 :     for (i = 1; i < co; i++) gel(a,k-1+i) = gel(x,i);
    1218        4235 :     gel(a,li) = zerocol(li-1); x = a;
    1219             :   }
    1220             :   else
    1221             :   {
    1222     3417906 :     x += co - li;
    1223     3417906 :     x[0] = evaltyp(t_MAT) | evallg(li); /* kill 0 columns */
    1224             :   }
    1225     3422141 :   if (moddiag)
    1226             :   { /* one column extra: an accumulator, discarded at the end */
    1227     3125412 :     if (lg(x) == li) x = shallowconcat(x, zerocol(li-1));
    1228             :     /* add up missing diag(dm) components */
    1229    12805290 :     for (i = li-1; i > 0; i--)
    1230             :     {
    1231     9679878 :       gcoeff(x, i, li) = gel(dm,i);
    1232    36065092 :       for (j = i; j > 0; j--)
    1233             :       {
    1234    26385214 :         GEN a = gcoeff(x, j, li);
    1235    26385214 :         if (!signe(a)) continue;
    1236     9719493 :         ZC_elem(a, gcoeff(x,j,j), x, NULL, li,j);
    1237     9719493 :         FpV_red_part_ipvec(gel(x,li), dm, j-1);
    1238     9719493 :         FpV_red_part_ipvec(gel(x,j),  dm, j-1);
    1239             :       }
    1240     9679878 :       if (gc_needed(av,1))
    1241             :       {
    1242           1 :         if (DEBUGMEM>1) pari_warn(warnmem,"ZM_hnfmod[2]. i=%ld", i);
    1243           1 :         x = gerepilecopy(av, x);
    1244             :       }
    1245             :     }
    1246             :   }
    1247             :   else
    1248             :   {
    1249      296729 :     b = gel(dm,1);
    1250     1189926 :     for (i = li-1; i > 0; i--)
    1251             :     {
    1252      893197 :       GEN d = bezout(gcoeff(x,i,i),b, &u,NULL);
    1253      893197 :       gcoeff(x,i,i) = d;
    1254      893197 :       FpV_Fp_mul_part_ip(gel(x,i), u, b, i-1);
    1255      893197 :       if (i > 1) b = diviiexact(b,d);
    1256             :     }
    1257             :   }
    1258     3422141 :   x[0] = evaltyp(t_MAT) | evallg(li); /* kill 0 columns / discard accumulator */
    1259     3422141 :   if (flag & hnf_PART) return x;
    1260             : 
    1261     3422141 :   if (!moddiag)
    1262             :   { /* compute optimal value for dm */
    1263      296729 :     b = cgetg(li, t_VEC); gel(b,1) = gcoeff(x,1,1);
    1264      296729 :     for (i = 2; i < li; i++) gel(b,i) = mulii(gel(b,i-1), gcoeff(x,i,i));
    1265      296729 :     dm = b;
    1266             :   }
    1267             : 
    1268    13995216 :   for (i = li-1; i > 0; i--)
    1269             :   {
    1270    10573075 :     GEN diag = gcoeff(x,i,i);
    1271    10573075 :     if (signe(diag) < 0) { gel(x,i) = ZC_neg(gel(x,i)); diag = gcoeff(x,i,i); }
    1272    28711105 :     for (j = i+1; j < li; j++)
    1273             :     {
    1274    18138030 :       b = gcoeff(x,i,j);
    1275    18138030 :       b = center? diviiround(b,diag): truedivii(b, diag);
    1276    18138030 :       if (!signe(b)) continue;
    1277     7376903 :       togglesign(b);
    1278     7376903 :       ZC_lincomb1_inplace(gel(x,j), gel(x,i),b);
    1279     7376903 :       p1 = gel(x,j);
    1280    29602307 :       for (k=1; k<i; k++)
    1281    22225404 :         if (lgefint(gel(p1,k)) > LDM[k]) gel(p1,k) = remii(gel(p1,k), gel(dm,i));
    1282             :     }
    1283    10573075 :     if (gc_needed(av,1))
    1284             :     {
    1285          32 :       if (DEBUGMEM>1) pari_warn(warnmem,"ZM_hnfmod[3]. i=%ld", i);
    1286          32 :       gerepileall(av, 2, &x, &dm); diag = gcoeff(x,i,i);
    1287             :     }
    1288             :   }
    1289     3422141 :   return x;
    1290             : }
    1291             : GEN
    1292     3418396 : ZM_hnfmodall(GEN x, GEN dm, long flag)
    1293             : {
    1294     3418396 :   pari_sp av = avma;
    1295     3418396 :   return gerepilecopy(av, ZM_hnfmodall_i(x, dm, flag));
    1296             : }
    1297             : GEN
    1298      296722 : ZM_hnfmod(GEN x, GEN d) { return ZM_hnfmodall(x,d,0); }
    1299             : GEN
    1300     3118594 : ZM_hnfmodid(GEN x, GEN d) { return ZM_hnfmodall(x,d,hnf_MODID); }
    1301             : 
    1302             : static GEN
    1303        3059 : allhnfmod(GEN x, GEN dm, int flag)
    1304             : {
    1305        3059 :   if (typ(x)!=t_MAT) pari_err_TYPE("allhnfmod",x);
    1306        3059 :   RgM_check_ZM(x, "allhnfmod");
    1307        3059 :   if (isintzero(dm)) return ZM_hnf(x);
    1308        3059 :   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        3045 : hnfmodid(GEN x, GEN d)
    1318             : {
    1319        3045 :   switch(typ(d))
    1320             :   {
    1321        3031 :     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        3045 :   return allhnfmod(x, d, hnf_MODID);
    1327             : }
    1328             : 
    1329             : /* M a ZM in HNF. Normalize with *centered* residues */
    1330             : GEN
    1331         959 : ZM_hnfcenter(GEN M)
    1332             : {
    1333         959 :   long i, j, k, N = lg(M)-1;
    1334         959 :   pari_sp av = avma;
    1335             : 
    1336        3640 :   for (j=N-1; j>0; j--) /* skip last line */
    1337             :   {
    1338        2681 :     GEN Mj = gel(M,j), a = gel(Mj,j);
    1339       15421 :     for (k = j+1; k <= N; k++)
    1340             :     {
    1341       12740 :       GEN Mk = gel(M,k), q = diviiround(gel(Mk,j), a);
    1342       12740 :       long s = signe(q);
    1343       12740 :       if (!s) continue;
    1344       10129 :       if (is_pm1(q))
    1345             :       {
    1346        1785 :         if (s < 0)
    1347         280 :           for (i = 1; i <= j; i++) gel(Mk,i) = addii(gel(Mk,i), gel(Mj,i));
    1348             :         else
    1349        1505 :           for (i = 1; i <= j; i++) gel(Mk,i) = subii(gel(Mk,i), gel(Mj,i));
    1350             :       }
    1351             :       else
    1352        8344 :         for (i = 1; i <= j; i++) gel(Mk,i) = subii(gel(Mk,i), mulii(q,gel(Mj,i)));
    1353       10129 :       if (gc_needed(av,1))
    1354             :       {
    1355           0 :         if (DEBUGMEM>1) pari_warn(warnmem,"ZM_hnfcenter, j = %ld",j);
    1356           0 :         M = gerepilecopy(av, M);
    1357             :       }
    1358             :     }
    1359             :   }
    1360         959 :   return M;
    1361             : }
    1362             : 
    1363             : /***********************************************************************/
    1364             : /*                                                                     */
    1365             : /*                 HNFLLL (Havas, Majewski, Mathews)                   */
    1366             : /*                                                                     */
    1367             : /***********************************************************************/
    1368             : 
    1369             : static void
    1370      712056 : Minus(long j, GEN lambda)
    1371             : {
    1372      712056 :   long k, n = lg(lambda);
    1373             : 
    1374      712056 :   for (k=1  ; k<j; k++) togglesign_safe(&gcoeff(lambda,k,j));
    1375      712056 :   for (k=j+1; k<n; k++) togglesign_safe(&gcoeff(lambda,j,k));
    1376      712056 : }
    1377             : 
    1378             : /* index of first non-zero entry */
    1379             : static long
    1380    34463394 : findi(GEN M)
    1381             : {
    1382    34463394 :   long i, n = lg(M);
    1383   210083709 :   for (i=1; i<n; i++)
    1384   191775895 :     if (signe(gel(M,i))) return i;
    1385    18307814 :   return 0;
    1386             : }
    1387             : 
    1388             : static long
    1389    34463394 : findi_normalize(GEN Aj, GEN B, long j, GEN lambda)
    1390             : {
    1391    34463394 :   long r = findi(Aj);
    1392    34463394 :   if (r && signe(gel(Aj,r)) < 0)
    1393             :   {
    1394      712056 :     ZV_togglesign(Aj); if (B) ZV_togglesign(gel(B,j));
    1395      712056 :     Minus(j,lambda);
    1396             :   }
    1397    34463394 :   return r;
    1398             : }
    1399             : 
    1400             : static void
    1401    17231540 : 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    17231540 :   *row0 = findi_normalize(gel(A,j), B,j,lambda);
    1407    17231540 :   *row1 = findi_normalize(gel(A,k), B,k,lambda);
    1408    17231540 :   if (*row0)
    1409     5717698 :     q = truedivii(gcoeff(A,*row0,k), gcoeff(A,*row0,j));
    1410    11513842 :   else if (absi_cmp(shifti(gcoeff(lambda,j,k), 1), gel(D,j)) > 0)
    1411     3247350 :     q = diviiround(gcoeff(lambda,j,k), gel(D,j));
    1412             :   else
    1413    25498032 :     return;
    1414             : 
    1415     8965048 :   if (signe(q))
    1416             :   {
    1417     5750530 :     GEN Lk = gel(lambda,k), Lj = gel(lambda,j);
    1418     5750530 :     togglesign_safe(&q);
    1419     5750530 :     if (*row0) ZC_lincomb1_inplace(gel(A,k),gel(A,j),q);
    1420     5750530 :     if (B) ZC_lincomb1_inplace(gel(B,k),gel(B,j),q);
    1421     5750530 :     gel(Lk,j) = addii(gel(Lk,j), mulii(q,gel(D,j)));
    1422     5750530 :     if (is_pm1(q))
    1423             :     {
    1424     1973085 :       if (signe(q) > 0)
    1425             :       {
    1426     5437933 :         for (i=1; i<j; i++)
    1427     4657753 :           if (signe(gel(Lj,i))) gel(Lk,i) = addii(gel(Lk,i), gel(Lj,i));
    1428             :       }
    1429             :       else
    1430             :       {
    1431     6727060 :         for (i=1; i<j; i++)
    1432     5534155 :           if (signe(gel(Lj,i))) gel(Lk,i) = subii(gel(Lk,i), gel(Lj,i));
    1433             :       }
    1434             :     }
    1435             :     else
    1436             :     {
    1437    19377757 :       for (i=1; i<j; i++)
    1438    15600312 :         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     3018622 : hnfswap(GEN A, GEN B, long k, GEN lambda, GEN D)
    1445             : {
    1446     3018622 :   GEN t, p1, p2, Lk = gel(lambda,k);
    1447     3018622 :   long i,j,n = lg(A);
    1448             : 
    1449     3018622 :   swap(gel(A,k), gel(A,k-1));
    1450     3018622 :   if (B) swap(gel(B,k), gel(B,k-1));
    1451     3018622 :   for (j=k-2; j; j--) swap(gcoeff(lambda,j,k-1), gel(Lk,j));
    1452    31996435 :   for (i=k+1; i<n; i++)
    1453             :   {
    1454    28977813 :     GEN Li = gel(lambda,i);
    1455    28977813 :     p1 = mulii(gel(Li,k-1), gel(D,k));
    1456    28977813 :     p2 = mulii(gel(Li,k), gel(Lk,k-1));
    1457    28977813 :     t = subii(p1,p2);
    1458             : 
    1459    28977813 :     p1 = mulii(gel(Li,k), gel(D,k-2));
    1460    28977813 :     p2 = mulii(gel(Li,k-1), gel(Lk,k-1));
    1461    28977813 :     gel(Li,k-1) = diviiexact(addii(p1,p2), gel(D,k-1));
    1462    28977813 :     gel(Li,k) = diviiexact(t, gel(D,k-1));
    1463             :   }
    1464     3018622 :   p1 = mulii(gel(D,k-2), gel(D,k));
    1465     3018622 :   p2 = sqri(gel(Lk,k-1));
    1466     3018622 :   gel(D,k-1) = diviiexact(addii(p1,p2), gel(D,k-1));
    1467     3018622 : }
    1468             : 
    1469             : /* reverse row order in matrix A, IN PLACE */
    1470             : static GEN
    1471      121980 : reverse_rows(GEN A)
    1472             : {
    1473      121980 :   long i, j, h, n = lg(A);
    1474      121980 :   if (n == 1) return A;
    1475      121966 :   h = lgcols(A);
    1476     1097448 :   for (j=1; j<n; j++)
    1477             :   {
    1478      975482 :     GEN c = gel(A,j);
    1479             :     /* start at (h-1) >>1 : if h = 2i even, no need to swap c[i] and itself */
    1480      975482 :     for (i=(h-1)>>1; i; i--) swap(gel(c,i), gel(c,h-i));
    1481             :   }
    1482      121966 :   return A;
    1483             : }
    1484             : 
    1485             : GEN
    1486       60990 : ZM_hnflll(GEN A, GEN *ptB, int remove)
    1487             : {
    1488       60990 :   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       60990 :   n = lg(A);
    1496       60990 :   A = reverse_rows(ZM_copy(A)); /* ZM_copy for in place findi_normalize() */
    1497       60990 :   B = ptB? matid(n-1): NULL;
    1498       60990 :   D = const_vec(n, gen_1) + 1;
    1499       60990 :   lambda = zeromatcopy(n-1,n-1);
    1500       60990 :   k = kmax = 2;
    1501     5451061 :   while (k < n)
    1502             :   {
    1503             :     long row0, row1;
    1504             :     int do_swap;
    1505     5329081 :     reduce2(A,B,k,k-1,&row0,&row1,lambda,D);
    1506     5329081 :     if (row0)
    1507     3178361 :       do_swap = (!row1 || row0 <= row1);
    1508     2150720 :     else if (!row1)
    1509             :     { /* row0 == row1 == 0 */
    1510     1842063 :       pari_sp av1 = avma;
    1511     1842063 :       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     1842063 :       do_swap = (cmpii(z, sqri(gel(D,k-1))) < 0);
    1516             : #endif
    1517     1842063 :       avma = av1;
    1518             :     }
    1519             :     else
    1520      308657 :       do_swap = 0;
    1521     5329081 :     if (do_swap)
    1522             :     {
    1523     2998492 :       hnfswap(A,B,k,lambda,D);
    1524     2998492 :       if (k > 2) k--;
    1525             :     }
    1526             :     else
    1527             :     {
    1528             :       long i;
    1529    14233048 :       for (i=k-2; i; i--)
    1530             :       {
    1531             :         long row0, row1;
    1532    11902459 :         reduce2(A,B,k,i,&row0,&row1,lambda,D);
    1533    11902459 :         if (gc_needed(av,3))
    1534             :         {
    1535          31 :           GEN b = D-1;
    1536          31 :           if (DEBUGMEM>1) pari_warn(warnmem,"hnflll (reducing), kmax = %ld",kmax);
    1537          31 :           gerepileall(av, B? 4: 3, &A, &lambda, &b, &B);
    1538          31 :           D = b+1;
    1539             :         }
    1540             :       }
    1541     2330589 :       if (++k > kmax) kmax = k;
    1542             :     }
    1543     5329081 :     if (gc_needed(av,3))
    1544             :     {
    1545          55 :       GEN b = D-1;
    1546          55 :       if (DEBUGMEM>1) pari_warn(warnmem,"hnflll, kmax = %ld / %ld",kmax,n-1);
    1547          55 :       gerepileall(av, B? 4: 3, &A, &lambda, &b, &B);
    1548          55 :       D = b+1;
    1549             :     }
    1550             :   }
    1551             :   /* handle trivial case: return negative diag coefficient otherwise */
    1552       60990 :   if (n == 2) (void)findi_normalize(gel(A,1), B,1,lambda);
    1553       60990 :   A = reverse_rows(A);
    1554       60990 :   if (remove)
    1555             :   {
    1556             :     long i;
    1557       21264 :     for (i = 1; i < n; i++)
    1558       20928 :       if (!ZV_equal0(gel(A,i))) break;
    1559        1776 :     remove_0cols(i-1, &A, &B, remove);
    1560             :   }
    1561       60990 :   gerepileall(av, B? 2: 1, &A, &B);
    1562       60990 :   if (B) *ptB = B;
    1563       60990 :   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      260799 : reduce1(GEN A, GEN B, long k, long j, GEN lambda, GEN D)
    1578             : {
    1579             :   GEN q;
    1580             :   long i;
    1581             : 
    1582      260799 :   if (signe(gel(A,j)))
    1583       16613 :     q = diviiround(gel(A,k),gel(A,j));
    1584      244186 :   else if (absi_cmp(shifti(gcoeff(lambda,j,k), 1), gel(D,j)) > 0)
    1585        2375 :     q = diviiround(gcoeff(lambda,j,k), gel(D,j));
    1586             :   else
    1587      502610 :     return;
    1588             : 
    1589       18988 :   if (signe(q))
    1590             :   {
    1591        8693 :     GEN Lk = gel(lambda,k), Lj = gel(lambda,j);
    1592        8693 :     togglesign_safe(&q);
    1593        8693 :     gel(A,k) = addii(gel(A,k), mulii(q,gel(A,j)));
    1594        8693 :     ZC_lincomb1_inplace(gel(B,k),gel(B,j),q);
    1595        8693 :     gel(Lk,j) = addii(gel(Lk,j), mulii(q,gel(D,j)));
    1596       22923 :     for (i=1; i<j; i++)
    1597       14230 :       if (signe(gel(Lj,i))) gel(Lk,i) = addii(gel(Lk,i), mulii(q,gel(Lj,i)));
    1598             :   }
    1599             : }
    1600             : 
    1601             : static GEN
    1602       14021 : 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       14021 :   long k, n = lg(A);
    1608             :   GEN B, lambda, D;
    1609             : 
    1610       14021 :   if (n == 1) retmkvec2(gen_1, cgetg(1,t_MAT));
    1611       14021 :   A = leafcopy(A);
    1612       14021 :   B = matid(n-1);
    1613       14021 :   lambda = zeromatcopy(n-1,n-1);
    1614       14021 :   D = const_vec(n, gen_1) + 1;
    1615       14021 :   k = 2;
    1616      120830 :   while (k < n)
    1617             :   {
    1618             :     int do_swap;
    1619             : 
    1620       92788 :     reduce1(A,B,k,k-1,lambda,D);
    1621       92788 :     if (signe(gel(A,k-1))) do_swap = 1;
    1622       76175 :     else if (!signe(gel(A,k)))
    1623             :     {
    1624       48901 :       pari_sp av1 = avma;
    1625       48901 :       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       48901 :       do_swap = (cmpii(z, sqri(gel(D,k-1))) < 0);
    1630             : #endif
    1631       48901 :       avma = av1;
    1632             :     }
    1633       27274 :     else do_swap = 0;
    1634             : 
    1635       92788 :     if (do_swap)
    1636             :     {
    1637       20130 :       hnfswap(A,B,k,lambda,D);
    1638       20130 :       if (k > 2) k--;
    1639             :     }
    1640             :     else
    1641             :     {
    1642             :       long i;
    1643       72658 :       for (i=k-2; i; i--) reduce1(A,B,k,i,lambda,D);
    1644       72658 :       k++;
    1645             :     }
    1646             :   }
    1647       14021 :   if (signe(gel(A,n-1)) < 0)
    1648             :   {
    1649        5703 :     gel(A,n-1) = negi(gel(A,n-1));
    1650        5703 :     ZV_togglesign(gel(B,n-1));
    1651             :   }
    1652       14021 :   return mkvec2(gel(A,n-1), B);
    1653             : }
    1654             : GEN
    1655       14007 : ZV_extgcd(GEN A)
    1656             : {
    1657       14007 :   pari_sp av = avma;
    1658       14007 :   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          21 : ZV_hnfgcdext(GEN A)
    1663             : {
    1664          21 :   pari_sp av = avma;
    1665             :   GEN z;
    1666          21 :   if (lg(A) == 1) retmkvec2(cgetg(1,t_MAT),cgetg(1,t_MAT));
    1667          14 :   z = ZV_gcdext_i(A);
    1668          14 :   gel(z,1) = mkmat(mkcol(gel(z,1)));
    1669          14 :   return gerepilecopy(av, z);
    1670             : }
    1671             : 
    1672             : /* HNF with permutation. */
    1673             : GEN
    1674         259 : ZM_hnfperm(GEN A, GEN *ptU, GEN *ptperm)
    1675             : {
    1676             :   GEN U, c, l, perm, d, p, q, b;
    1677         259 :   pari_sp av = avma, av1;
    1678             :   long r, t, i, j, j1, k, m, n;
    1679             : 
    1680         259 :   n = lg(A)-1;
    1681         259 :   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         252 :   m = nbrows(A);
    1688         252 :   c = zero_zv(m);
    1689         252 :   l = zero_zv(n);
    1690         252 :   perm = cgetg(m+1, t_VECSMALL);
    1691         252 :   av1 = avma;
    1692         252 :   A = RgM_shallowcopy(A);
    1693         252 :   U = ptU? matid(n): NULL;
    1694             :   /* U base change matrix : A0*U = A all along */
    1695        1246 :   for (r=0, k=1; k <= n; k++)
    1696             :   {
    1697        3038 :     for (j=1; j<k; j++)
    1698             :     {
    1699        2044 :       if (!l[j]) continue;
    1700        1911 :       t=l[j]; b=gcoeff(A,t,k);
    1701        1911 :       if (!signe(b)) continue;
    1702             : 
    1703         553 :       ZC_elem(b,gcoeff(A,t,j), A,U,k,j);
    1704         553 :       d = gcoeff(A,t,j);
    1705         553 :       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        1309 :       for (j1=1; j1<j; j1++)
    1712             :       {
    1713         756 :         if (!l[j1]) continue;
    1714         742 :         q = truedivii(gcoeff(A,t,j1),d);
    1715         742 :         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         994 :     t = m; while (t && (c[t] || !signe(gcoeff(A,t,k)))) t--;
    1723         994 :     if (t)
    1724             :     {
    1725         896 :       p = gcoeff(A,t,k);
    1726       15365 :       for (i=t-1; i; i--)
    1727             :       {
    1728       14469 :         q = gcoeff(A,i,k);
    1729       14469 :         if (signe(q) && absi_cmp(p,q) > 0) { p = q; t = i; }
    1730             :       }
    1731         896 :       perm[++r] = l[k] = t; c[t] = k;
    1732         896 :       if (signe(p) < 0)
    1733             :       {
    1734         153 :         ZV_neg_inplace(gel(A,k));
    1735         153 :         if (U) ZV_togglesign(gel(U,k));
    1736         153 :         p = gcoeff(A,t,k);
    1737             :       }
    1738             :       /* p > 0 */
    1739        2597 :       for (j=1; j<k; j++)
    1740             :       {
    1741        1701 :         if (!l[j]) continue;
    1742        1680 :         q = truedivii(gcoeff(A,t,j),p);
    1743        1680 :         if (!signe(q)) continue;
    1744             : 
    1745         271 :         togglesign(q);
    1746         271 :         ZC_lincomb1_inplace(gel(A,j), gel(A,k), q);
    1747         271 :         if (U) ZC_lincomb1_inplace(gel(U,j), gel(U,k), q);
    1748             :       }
    1749             :     }
    1750         994 :     if (gc_needed(av1,1))
    1751             :     {
    1752           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"hnfperm, k=%ld",k);
    1753           0 :       gerepileall(av1, U? 2: 1, &A, &U);
    1754             :     }
    1755             :   }
    1756         252 :   if (r < m)
    1757             :   {
    1758        3640 :     for (i=1,k=r; i<=m; i++)
    1759        3444 :       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         252 :   p = cgetg(r+1,t_MAT);
    1766         252 :   for (i=1; i<=m/2; i++) lswap(perm[i], perm[m+1-i]);
    1767         252 :   if (U)
    1768             :   {
    1769          70 :     GEN u = cgetg(n+1,t_MAT);
    1770         322 :     for (t=1,k=r,j=1; j<=n; j++)
    1771         252 :       if (l[j])
    1772             :       {
    1773         154 :         u[k + n-r] = U[j];
    1774         154 :         gel(p,k--) = vecpermute(gel(A,j), perm);
    1775             :       }
    1776             :       else
    1777          98 :         u[t++] = U[j];
    1778          70 :     *ptU = u;
    1779          70 :     if (ptperm) *ptperm = perm;
    1780          70 :     gerepileall(av, ptperm? 3: 2, &p, ptU, ptperm);
    1781             :   }
    1782             :   else
    1783             :   {
    1784         924 :     for (k=r,j=1; j<=n; j++)
    1785         742 :       if (l[j]) gel(p,k--) = vecpermute(gel(A,j), perm);
    1786         182 :     if (ptperm) *ptperm = perm;
    1787         182 :     gerepileall(av, ptperm? 2: 1, &p, ptperm);
    1788             :   }
    1789         252 :   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      104669 : ZM_hnfall(GEN A, GEN *ptB, long remove)
    1805             : {
    1806      104669 :   pari_sp av = avma, av1;
    1807             :   long m, n, r, i, j, k, li;
    1808             :   GEN B, c, h, a;
    1809             : 
    1810      104669 :   RgM_dimensions(A, &m,&n);
    1811      104669 :   if (!n)
    1812             :   {
    1813           7 :     if (ptB) *ptB = cgetg(1,t_MAT);
    1814           7 :     return cgetg(1,t_MAT);
    1815             :   }
    1816      104662 :   c = zero_zv(m);
    1817      104662 :   h = const_vecsmall(n, m);
    1818      104662 :   av1 = avma;
    1819      104662 :   A = RgM_shallowcopy(A);
    1820      104662 :   B = ptB? matid(n): NULL;
    1821      104662 :   r = n+1;
    1822      371455 :   for (li=m; li; li--)
    1823             :   {
    1824      852996 :     for (j=1; j<r; j++)
    1825             :     {
    1826     1108252 :       for (i=h[j]; i>li; i--)
    1827             :       {
    1828      255418 :         a = gcoeff(A,i,j);
    1829      255418 :         k = c[i];
    1830             :         /* zero a = Aij  using  Aik */
    1831      255418 :         if (signe(a)) ZC_elem(a,gcoeff(A,i,k), A,B,j,k);
    1832      255418 :         ZM_reduce(A,B, i,k); /* ensure reduced entries */
    1833      255418 :         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      852834 :       if (signe( gcoeff(A,li,j) )) break;
    1840      586203 :       h[j] = li-1;
    1841             :     }
    1842      266793 :     if (j == r) continue;
    1843      266631 :     r--;
    1844      266631 :     if (j < r) /* A[j] != 0 */
    1845             :     {
    1846      179004 :       swap(gel(A,j), gel(A,r));
    1847      179004 :       if (B) swap(gel(B,j), gel(B,r));
    1848      179004 :       h[j] = h[r]; h[r] = li; c[li] = r;
    1849             :     }
    1850      266631 :     if (signe(gcoeff(A,li,r)) < 0)
    1851             :     {
    1852       52618 :       ZV_neg_inplace(gel(A,r));
    1853       52618 :       if (B) ZV_togglesign(gel(B,r));
    1854             :     }
    1855      266631 :     ZM_reduce(A,B, li,r);
    1856      266631 :     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      104662 :   if (DEBUGLEVEL>5) err_printf("\nhnfall, final phase: ");
    1864      104662 :   r--; /* first r cols are in the image the n-r (independent) last ones */
    1865      628327 :   for (j=1; j<=r; j++)
    1866     1211622 :     for (i=h[j]; i; i--)
    1867             :     {
    1868      687957 :       a = gcoeff(A,i,j);
    1869      687957 :       k = c[i];
    1870      687957 :       if (signe(a)) ZC_elem(a,gcoeff(A,i,k), A,B, j,k);
    1871      687957 :       ZM_reduce(A,B, i,k); /* ensure reduced entries, even if a = 0 */
    1872      687957 :       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      104662 :   if (DEBUGLEVEL>5) err_printf("\n");
    1879      104662 :   if (remove) remove_0cols(r, &A, &B, remove);
    1880      104662 :   gerepileall(av, B? 2: 1, &A, &B);
    1881      104662 :   if (B) *ptB = B;
    1882      104662 :   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       28929 : hnf_divscale(GEN A, GEN B, GEN t)
    1898             : {
    1899       28929 :   long n = lg(A)-1, i,j,k;
    1900       28929 :   GEN m, c = cgetg(n+1,t_MAT);
    1901             : 
    1902       28929 :   if (!n) return c;
    1903      108926 :   for (k=1; k<=n; k++)
    1904             :   {
    1905       79997 :     GEN u = cgetg(n+1, t_COL), b = gel(B,k);
    1906       79997 :     pari_sp av = avma;
    1907       79997 :     gel(c,k) = u; m = mulii(gel(b,n),t);
    1908       79997 :     gel(u,n) = gerepileuptoint(av, diviiexact(m, gcoeff(A,n,n)));
    1909      387593 :     for (i=n-1; i>0; i--)
    1910             :     {
    1911      307596 :       av = avma; m = mulii(gel(b,i),t);
    1912      307596 :       for (j=i+1; j<=n; j++) m = subii(m, mulii(gcoeff(A,i,j),gel(u,j)));
    1913      307596 :       gel(u,i) = gerepileuptoint(av, diviiexact(m, gcoeff(A,i,i)));
    1914             :     }
    1915             :   }
    1916       28929 :   return c;
    1917             : }
    1918             : 
    1919             : /* A, B integral upper HNF. A^(-1) B integral ? */
    1920             : int
    1921          91 : hnfdivide(GEN A, GEN B)
    1922             : {
    1923          91 :   pari_sp av = avma;
    1924          91 :   long n = lg(A)-1, i,j,k;
    1925             :   GEN u, b, m, r;
    1926             : 
    1927          91 :   if (!n) return 1;
    1928          91 :   if (lg(B)-1 != n) pari_err_DIM("hnfdivide");
    1929          91 :   u = cgetg(n+1, t_COL);
    1930         413 :   for (k=1; k<=n; k++)
    1931             :   {
    1932         322 :     b = gel(B,k);
    1933         322 :     m = gel(b,k);
    1934         322 :     gel(u,k) = dvmdii(m, gcoeff(A,k,k), &r);
    1935         322 :     if (r != gen_0) { avma = av; return 0; }
    1936         798 :     for (i=k-1; i>0; i--)
    1937             :     {
    1938         476 :       m = gel(b,i);
    1939         476 :       for (j=i+1; j<=k; j++) m = subii(m, mulii(gcoeff(A,i,j),gel(u,j)));
    1940         476 :       m = dvmdii(m, gcoeff(A,i,i), &r);
    1941         476 :       if (r != gen_0) { avma = av; return 0; }
    1942         476 :       gel(u,i) = m;
    1943             :     }
    1944             :   }
    1945          91 :   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       93758 : hnf_invimage(GEN A, GEN b)
    1952             : {
    1953       93758 :   pari_sp av = avma;
    1954       93758 :   long n = lg(A)-1, m, i, k;
    1955             :   GEN u, r;
    1956             : 
    1957       93758 :   if (!n) return lg(b)==1? cgetg(1,t_COL):NULL;
    1958       93716 :   m = nbrows(A); /* m >= n */
    1959       93716 :   u = cgetg(n+1, t_COL);
    1960      617512 :   for (i = n, k = m; k > 0; k--)
    1961             :   {
    1962      617512 :     pari_sp av2 = avma;
    1963             :     long j;
    1964      617512 :     GEN t = gel(b,k), Aki = gcoeff(A,k,i);
    1965      617512 :     if (typ(t) != t_INT) pari_err_TYPE("hnf_invimage",t);
    1966      617512 :     for (j=i+1; j<=n; j++) t = subii(t, mulii(gcoeff(A,k,j),gel(u,j)));
    1967      617512 :     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      617512 :     t = dvmdii(t, Aki, &r);
    1973      617512 :     if (r != gen_0) { avma = av; return NULL; }
    1974      608818 :     gel(u,i) = gerepileuptoint(av2, t);
    1975      608818 :     if (--i == 0) break;
    1976             :   }
    1977             :   /* If there is a solution, it must be u. Check remaining equations */
    1978      170044 :   for (; k > 0; k--)
    1979             :   {
    1980       85022 :     pari_sp av2 = avma;
    1981             :     long j;
    1982       85022 :     GEN t = gel(b,k);
    1983       85022 :     if (typ(t) != t_INT) pari_err_TYPE("hnf_invimage",t);
    1984       85022 :     for (j=1; j<=n; j++) t = subii(t, mulii(gcoeff(A,k,j),gel(u,j)));
    1985       85022 :     if (signe(t)) { avma = av;return NULL; }
    1986       85022 :     avma = av2;
    1987             :   }
    1988       85022 :   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       24458 : hnf_solve(GEN A, GEN B)
    1995             : {
    1996             :   pari_sp av;
    1997             :   long i, l;
    1998             :   GEN C;
    1999             : 
    2000       24458 :   if (typ(B) == t_COL) return hnf_invimage(A, B);
    2001       20727 :   av = avma; C = cgetg_copy(B, &l);
    2002      103026 :   for (i = 1; i < l; i++) {
    2003       86457 :     GEN c = hnf_invimage(A, gel(B,i));
    2004       86457 :     if (!c) { avma = av; return NULL; }
    2005       82299 :     gel(C,i) = c;
    2006             :   }
    2007       16569 :   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       68267 : snf_pile(pari_sp av, GEN *x, GEN *U, GEN *V)
    2029             : {
    2030             :   GEN *gptr[3];
    2031       68267 :   int c = 1; gptr[0]=x;
    2032       68267 :   if (*U) gptr[c++] = U;
    2033       68267 :   if (*V) gptr[c++] = V;
    2034       68267 :   gerepilemany(av,gptr,c);
    2035       68267 : }
    2036             : 
    2037             : static GEN
    2038       86997 : bezout_step(GEN *pa, GEN *pb, GEN *pu, GEN *pv)
    2039             : {
    2040       86997 :   GEN a = *pa, b = *pb, d;
    2041       86997 :   if (absi_equal(a,b))
    2042             :   {
    2043       30950 :     long sa = signe(a), sb = signe(b);
    2044       30950 :     *pv = gen_0;
    2045       30950 :     if (sb == sa) {
    2046       30754 :       *pa = *pb = gen_1;
    2047       30754 :       if (sa > 0) {*pu=gen_1; return a;} else {*pu=gen_m1; return absi(a);}
    2048             :     }
    2049         196 :     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       56047 :   d = bezout(a,b, pu,pv);
    2053       56047 :   *pa = diviiexact(a, d);
    2054       56047 :   *pb = diviiexact(b, d); return d;
    2055             : }
    2056             : 
    2057             : static int
    2058      266690 : 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      152231 : ZM_snf_no_divide(GEN x, long i)
    2064             : {
    2065      152231 :   GEN b = gcoeff(x,i,i);
    2066             :   long j, k;
    2067             : 
    2068      152231 :   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      152231 :   if (is_pm1(b)) return 0;
    2076      313437 :   for (k=1; k<i; k++)
    2077             :   {
    2078     1113082 :     for (j=1; j<i; j++)
    2079      888053 :       if (signe(remii(gcoeff(x,k,j),b))) return k;
    2080             :   }
    2081       84784 :   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       75933 : ZM_snfall_i(GEN x, GEN *ptU, GEN *ptV, int return_vec)
    2088             : {
    2089       75933 :   pari_sp av0 = avma, av;
    2090             :   long i, j, k, m0, m, n0, n;
    2091       75933 :   GEN p1, u, v, U, V, V0, mdet, ys, perm = NULL;
    2092             : 
    2093       75933 :   n0 = n = lg(x)-1;
    2094       75933 :   if (!n) {
    2095        8506 :     if (ptU) *ptU = cgetg(1,t_MAT);
    2096        8506 :     if (ptV) *ptV = cgetg(1,t_MAT);
    2097        8506 :     return cgetg(1, return_vec? t_VEC: t_MAT);
    2098             :   }
    2099       67427 :   av = avma;
    2100       67427 :   m0 = m = nbrows(x);
    2101             : 
    2102       67427 :   U = ptU? gen_1: NULL; /* TRANSPOSE of row transform matrix [act on columns]*/
    2103       67427 :   V = ptV? gen_1: NULL;
    2104       67427 :   V0 = NULL;
    2105       67427 :   x = RgM_shallowcopy(x);
    2106       67427 :   if (m == n && ZM_ishnf(x))
    2107             :   {
    2108       60924 :     mdet = ZM_det_triangular(x);
    2109       60924 :     if (V) *ptV = matid(n);
    2110             :   }
    2111             :   else
    2112             :   {
    2113        6503 :     mdet = ZM_detmult(x);
    2114        6503 :     if (signe(mdet))
    2115             :     {
    2116        6482 :       if (!V)
    2117          77 :         p1 = ZM_hnfmod(x,mdet);
    2118             :       else
    2119             :       {
    2120        6405 :         if (m == n)
    2121             :         {
    2122        6384 :           p1 = ZM_hnfmod(x,mdet);
    2123        6384 :           *ptV = RgM_solve(x,p1);
    2124             :         }
    2125             :         else
    2126          21 :           p1 = ZM_hnfperm(x, ptV, ptU? &perm: NULL);
    2127             :       }
    2128        6482 :       mdet = ZM_det_triangular(p1);
    2129             :     }
    2130             :     else
    2131          21 :       p1 = ZM_hnfperm(x, ptV, ptU? &perm: NULL);
    2132        6503 :     x = p1;
    2133             :   }
    2134       67427 :   n = lg(x)-1;
    2135       67427 :   if (V)
    2136             :   {
    2137       58040 :     V = *ptV;
    2138       58040 :     if (n != n0)
    2139             :     {
    2140          21 :       V0 = vecslice(V, 1, n0 - n); /* kernel */
    2141          21 :       V  = vecslice(V, n0-n+1, n0);
    2142          21 :       av = avma;
    2143             :     }
    2144             :   }
    2145             :   /* independent columns */
    2146       67427 :   if (!signe(mdet))
    2147             :   {
    2148          21 :     if (n)
    2149             :     {
    2150          21 :       x = ZM_snfall_i(shallowtrans(x), ptV, ptU, return_vec); /* swap ptV,ptU */
    2151          21 :       if (typ(x) == t_MAT && n != m) x = shallowtrans(x);
    2152          21 :       if (V) V = ZM_mul(V, shallowtrans(*ptV));
    2153          21 :       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          21 :     goto THEEND;
    2162             :   }
    2163       67406 :   if (U) U = matid(n);
    2164             : 
    2165             :   /* square, maximal rank n */
    2166       67406 :   p1 = gen_indexsort(RgM_diagonal_shallow(x), NULL, &negcmpii);
    2167       67406 :   ys = cgetg(n+1,t_MAT);
    2168       67406 :   for (j=1; j<=n; j++) gel(ys,j) = vecpermute(gel(x, p1[j]), p1);
    2169       67406 :   x = ys;
    2170       67406 :   if (U) U = vecpermute(U, p1);
    2171       67406 :   if (V) V = vecpermute(V, p1);
    2172             : 
    2173       67406 :   p1 = ZM_hnfmod(x, mdet);
    2174       67406 :   if (V) V = ZM_mul(V, RgM_solve(x,p1));
    2175       67406 :   x = p1;
    2176             : 
    2177       67406 :   if (DEBUGLEVEL>7) err_printf("starting SNF loop");
    2178      432026 :   for (i=n; i>1; i--)
    2179             :   {
    2180      148607 :     if (DEBUGLEVEL>7) err_printf("\ni = %ld: ",i);
    2181             :     for(;;)
    2182             :     {
    2183      199112 :       int c = 0;
    2184             :       GEN a, b;
    2185      984188 :       for (j=i-1; j>=1; j--)
    2186             :       {
    2187      785076 :         b = gcoeff(x,i,j); if (!signe(b)) continue;
    2188       10245 :         a = gcoeff(x,i,i);
    2189       10245 :         ZC_elem(b, a, x,V, j,i);
    2190       10245 :         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      199112 :       if (DEBUGLEVEL>7) err_printf("; ");
    2197      984188 :       for (j=i-1; j>=1; j--)
    2198             :       {
    2199             :         GEN d;
    2200      785076 :         b = gcoeff(x,j,i); if (!signe(b)) continue;
    2201       86997 :         a = gcoeff(x,i,i);
    2202       86997 :         d = bezout_step(&a, &b, &u, &v);
    2203      641201 :         for (k = 1; k < i; k++)
    2204             :         {
    2205      554204 :           GEN t = addii(mulii(u,gcoeff(x,i,k)),mulii(v,gcoeff(x,j,k)));
    2206     1108408 :           gcoeff(x,j,k) = subii(mulii(a,gcoeff(x,j,k)),
    2207      554204 :                                 mulii(b,gcoeff(x,i,k)));
    2208      554204 :           gcoeff(x,i,k) = t;
    2209             :         }
    2210       86997 :         gcoeff(x,j,i) = gen_0;
    2211       86997 :         gcoeff(x,i,i) = d;
    2212       86997 :         if (U) update(u,v,a,b,(GEN*)(U+i),(GEN*)(U+j));
    2213       86997 :         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       86997 :         c = 1;
    2219             :       }
    2220      199112 :       if (!c)
    2221             :       {
    2222      152231 :         k = ZM_snf_no_divide(x, i);
    2223      152231 :         if (!k) break;
    2224             : 
    2225             :         /* x[k,j] != 0 mod b */
    2226       14757 :         for (j=1; j<=i; j++)
    2227       11133 :           gcoeff(x,i,j) = addii(gcoeff(x,i,j),gcoeff(x,k,j));
    2228        3624 :         if (U) gel(U,i) = gadd(gel(U,i),gel(U,k));
    2229             :       }
    2230       50505 :       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       50505 :     }
    2236             :   }
    2237       67406 :   if (DEBUGLEVEL>7) err_printf("\n");
    2238      283419 :   for (k=1; k<=n; k++)
    2239      216013 :     if (signe(gcoeff(x,k,k)) < 0)
    2240             :     {
    2241          98 :       if (V) ZV_togglesign(gel(V,k));
    2242          98 :       togglesign(gcoeff(x,k,k));
    2243             :     }
    2244             : THEEND:
    2245       67427 :   if (return_vec)
    2246             :   {
    2247       66386 :     long l = lg(x)-1;
    2248       66386 :     if (typ(x) == t_MAT) x = RgM_diagonal_shallow(x);
    2249       66386 :     if (m0 > l) x = shallowconcat(zerovec(m0-l), x);
    2250             :   }
    2251             : 
    2252       67427 :   if (V0)
    2253             :   {
    2254          21 :     if (!return_vec) x = shallowconcat(zeromat(m,n0-n), x);
    2255          21 :     if (V) V = shallowconcat(V0, V);
    2256             :   }
    2257       67427 :   if (U)
    2258             :   {
    2259       32807 :     U = shallowtrans(U);
    2260       32807 :     if (perm) U = vecpermute(U, perm_inv(perm));
    2261             :   }
    2262       67427 :   snf_pile(av0, &x,&U,&V);
    2263       67427 :   if (ptU) *ptU = U;
    2264       67427 :   if (ptV) *ptV = V;
    2265       67427 :   return x;
    2266             : }
    2267             : GEN
    2268        1727 : 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          21 : smithall(GEN x)
    2280             : {
    2281          21 :   GEN z = cgetg(4, t_VEC);
    2282          21 :   if (typ(x)!=t_MAT) pari_err_TYPE("smithall",x);
    2283          21 :   RgM_check_ZM(x, "smithall");
    2284          21 :   gel(z,3) = ZM_snfall_i(x, (GEN*)(z+1),(GEN*)(z+2), 0);
    2285          21 :   return z;
    2286             : }
    2287             : 
    2288             : void
    2289       71924 : ZM_snfclean(GEN d, GEN u, GEN v)
    2290             : {
    2291       71924 :   long i, c, l = lg(d);
    2292             : 
    2293       71924 :   if (typ(d) == t_VEC)
    2294       71924 :     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       71924 :   setlg(d, c);
    2301       71924 :   if (u) for (i=1; i<l; i++) setlg(gel(u,i), c);
    2302       71924 :   if (v) setlg(v, c);
    2303       71924 : }
    2304             : 
    2305             : /* Assume z was computed by [g]smithall(). Remove the 1s on the diagonal */
    2306             : GEN
    2307         245 : smithclean(GEN z)
    2308             : {
    2309             :   long i, j, h, l, c, d;
    2310             :   GEN U, V, y, D, t;
    2311             : 
    2312         245 :   if (typ(z) != t_VEC) pari_err_TYPE("smithclean",z);
    2313         245 :   l = lg(z); if (l == 1) return cgetg(1,t_VEC);
    2314         238 :   U = gel(z,1);
    2315         238 :   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          21 :   V = gel(z,2);
    2322          21 :   D = gel(z,3);
    2323          21 :   l = lg(D);
    2324          21 :   if (l == 1) return gcopy(z);
    2325          21 :   h = lgcols(D);
    2326          21 :   if (h > l)
    2327             :   { /* D = vconcat(zero matrix, diagonal matrix) */
    2328          21 :     for (c=1+h-l, d=1; c<h; c++,d++)
    2329          21 :       if (gequal1(gcoeff(D,c,d))) break;
    2330             :   }
    2331           7 :   else if (h < l)
    2332             :   { /* D = concat(zero matrix, diagonal matrix) */
    2333           7 :     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           0 :     for (c=1; c<l; c++)
    2339           0 :       if (gequal1(gcoeff(D,c,c))) break;
    2340           0 :     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          21 :   y = cgetg(4,t_VEC);
    2344             :   /* truncate U to (c-1) x (h-1) */
    2345          21 :   gel(y,1) = t = cgetg(h,t_MAT);
    2346          21 :   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          21 :   gel(y,2) = gcopy_lg(V, d);
    2349          21 :   gel(y,3) = t = zeromatcopy(c-1, d-1);
    2350             :   /* truncate D to a (c-1) x (d-1) matrix */
    2351          21 :   if (d > 1)
    2352             :   {
    2353          14 :     if (h > l)
    2354             :     {
    2355          14 :       for (i=1+h-l, j=1; i<c; i++,j++)
    2356           7 :         gcoeff(t,i,j) = gcopy(gcoeff(D,i,j));
    2357             :     }
    2358           7 :     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           0 :       for (j=1; j<d; j++)
    2366           0 :         gcoeff(t,j,j) = gcopy(gcoeff(D,j,j));
    2367             :     }
    2368             :   }
    2369          21 :   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          49 :   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         252 :   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          35 :   return gerepilecopy(av, z);
    2537             : }
    2538             : 
    2539             : GEN
    2540          84 : matsnf0(GEN x,long flag)
    2541             : {
    2542          84 :   pari_sp av = avma;
    2543          84 :   if (flag > 7) pari_err_FLAG("matsnf");
    2544          84 :   if (typ(x) == t_VEC && flag & 4) return smithclean(x);
    2545          84 :   if (flag & 2) x = flag&1 ? gsmithall(x): gsmith(x);
    2546          42 :   else          x = flag&1 ?  smithall(x):  smith(x);
    2547          84 :   if (flag & 4) x = gerepileupto(av, smithclean(x));
    2548          84 :   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 is a relation matrix, either in HNF or a t_VEC (diagonal HNF) */
    2558             : static GEN
    2559       71637 : snf_group(GEN H, GEN D, GEN *newU, GEN *newUi)
    2560             : {
    2561             :   long i, j, l;
    2562             : 
    2563       71637 :   ZM_snfclean(D, newU? *newU: NULL, newUi? *newUi: NULL);
    2564       71637 :   l = lg(D);
    2565       71637 :   if (newU) {
    2566       29484 :     GEN U = *newU;
    2567      127729 :     for (i = 1; i < l; i++)
    2568             :     {
    2569       98245 :       GEN d = gel(D,i), d2 = shifti(d, 1);
    2570      833840 :       for (j = 1; j < lg(U); j++)
    2571      735595 :         gcoeff(U,i,j) = centermodii(gcoeff(U,i,j), d, d2);
    2572             :     }
    2573       29484 :     *newU = U;
    2574             :   }
    2575       71637 :   if (newUi) { /* UHV=D -> U^-1 = (HV)D^-1 -> U^-1 = H(VD^-1 mod 1) mod H */
    2576       65330 :     if (l > 1)
    2577             :     { /* Ui = ZM_inv(U, gen_1); setlg(Ui, l); */
    2578       57510 :       GEN V = *newUi, Ui;
    2579       57510 :       for (i = 1; i < l; i++) gel(V,i) = FpC_red(gel(V,i), gel(D,i));
    2580       57510 :       Ui = typ(H) == t_VEC? ZM_diag_mul(H, V): ZM_mul(H, V);
    2581       57510 :       for (i = 1; i < l; i++) gel(Ui,i) = ZC_Z_divexact(gel(Ui,i), gel(D,i));
    2582       57510 :       if (typ(H) == t_VEC)
    2583         553 :       { for (i = 1; i < l; i++) gel(Ui,i) = vecmodii(gel(Ui,i), H); }
    2584             :       else
    2585       56957 :         Ui = ZM_hnfrem(Ui, H);
    2586       57510 :       *newUi = Ui;
    2587             :     }
    2588             :   }
    2589       71637 :   return D;
    2590             : }
    2591             : /* H relation matrix among row of generators g in HNF.  Let URV = D its SNF,
    2592             :  * newU R newV = newD its clean SNF (no 1 in Dnew). Return the diagonal of
    2593             :  * newD, newU and newUi such that  1/U = (newUi, ?).
    2594             :  * Rationale: let (G,0) = g Ui be the new generators then
    2595             :  * 0 = G U R --> G D = 0,  g = G newU  and  G = g newUi */
    2596             : GEN
    2597       71084 : ZM_snf_group(GEN H, GEN *newU, GEN *newUi)
    2598             : {
    2599       71084 :   GEN D = ZM_snfall_i(H, newU, newUi, 1);
    2600       71084 :   return snf_group(H, D, newU, newUi);
    2601             : }
    2602             : 
    2603             : /* D a ZV: SNF for matdiagonal(D). Faster because we only ensure elementary
    2604             :  * divisors condition: d[n] | ... | d[1] and need not convert D to matrix form*/
    2605             : GEN
    2606         840 : ZV_snfall(GEN D, GEN *pU, GEN *pV)
    2607             : {
    2608         840 :   pari_sp av = avma;
    2609         840 :   long j, n = lg(D)-1;
    2610         840 :   GEN U = pU? matid(n): NULL;
    2611         840 :   GEN V = pV? matid(n): NULL;
    2612             :   GEN p;
    2613             : 
    2614         840 :   D = leafcopy(D);
    2615        2345 :   for (j = n; j > 0; j--)
    2616             :   {
    2617        1505 :     GEN b = gel(D,j);
    2618        1505 :     if (signe(b) < 0)
    2619             :     {
    2620           0 :       gel(D,j) = absi(b);
    2621           0 :       if (V) ZV_togglesign(gel(V,j));
    2622             :     }
    2623             :   }
    2624             :   /* entries are non-negative integers */
    2625         840 :   p = gen_indexsort(D, NULL, &negcmpii);
    2626         840 :   D = vecpermute(D, p);
    2627         840 :   if (U) U = vecpermute(U, p);
    2628         840 :   if (V) V = vecpermute(V, p);
    2629             :   /* entries are sorted by decreasing value */
    2630        2345 :   for (j = n; j > 0; j--)
    2631             :   {
    2632        1505 :     GEN b = gel(D,j);
    2633             :     long i;
    2634        2275 :     for (i = j-1; i > 0; i--)
    2635             :     { /* invariant: a >= b. If au+bv = d is a Bezout relation, A=a/d and B=b/d
    2636             :        * we have [B,-A;u,v]*diag(a,b)*[1-u*A,1; -u*A,1]] = diag(Ab, d) */
    2637         847 :       GEN a = gel(D,i), u,v, d = bezout(a,b, &u,&v), A, Wi, Wj;
    2638         847 :       if (equalii(d,b)) continue;
    2639         161 :       A = diviiexact(a,d);
    2640         161 :       if (V)
    2641             :       {
    2642         126 :         GEN t = mulii(u,A);
    2643         126 :         Wi = ZC_lincomb(subui(1,t), negi(t), gel(V,i), gel(V,j));
    2644         126 :         Wj = ZC_add(gel(V,i), gel(V,j));
    2645         126 :         gel(V,i) = Wi;
    2646         126 :         gel(V,j) = Wj;
    2647             :       }
    2648         161 :       if (U)
    2649             :       {
    2650         161 :         GEN B = diviiexact(b,d);
    2651         161 :         Wi = ZC_lincomb(B, negi(A), gel(U,i), gel(U,j));
    2652         161 :         Wj = ZC_lincomb(u, v, gel(U,i), gel(U,j));
    2653         161 :         gel(U,i) = Wi;
    2654         161 :         gel(U,j) = Wj;
    2655             :       }
    2656         161 :       gel(D,i) = mulii(A,b); /* lcm(a,b) */
    2657         161 :       gel(D,j) = d; /* gcd(a,b) */
    2658         161 :       b = gel(D,j); if (equali1(b)) break;
    2659             :     }
    2660             :   }
    2661         840 :   snf_pile(av, &D,&U,&V);
    2662         840 :   if (U) *pU = shallowtrans(U);
    2663         840 :   if (V) *pV = V;
    2664         840 :   return D;
    2665             : }
    2666             : GEN
    2667         553 : ZV_snf_group(GEN d, GEN *newU, GEN *newUi)
    2668             : {
    2669         553 :   GEN D = ZV_snfall(d, newU, newUi);
    2670         553 :   return snf_group(d, D, newU, newUi);
    2671             : }

Generated by: LCOV version 1.11